Velocity estimation from imagery using symmetric displaced frame difference equation

Abstract

A method and apparatus for processing an image sequence described herein that provides a symmetric displaced frame difference equation. Next, an input image sequence can be received that includes a pair of image frames individually including multidimensional image data corresponding to a plurality of pixel locations at different times. Finally, using at least one processor, the symmetric displaced frame difference equation can be solved using an iteration equation and the image data of the pair of image frames to determine a displacement field describing displacement vectors at half of the total displacement vector.

Claims

1. A method for processing an image sequence, comprising; providing a symmetric displaced frame difference equation; receiving an input image sequence comprising a pair of image frames individually including multidimensional image data corresponding to a plurality of pixel locations at different times; using at least one processor, solving the symmetric displaced frame difference equation using: an iteration equation and the image data of the pair of image frames to determine a displacement field describing displacement vectors at half of the total displacement vector, and at least one bilinear polynomial function expressing a multidimensional displacement field; and determining, using the at least one processor, displacement by motion for the input image sequence based on the solved symmetric displaced frame difference equation.

2. The method of claim 1, wherein the multidimensional image data of the pair of frames is two-dimensional, and wherein the symmetric displaced frame difference equation is
I(i+u.sub.ijt.sub.h,j+v.sub.ijt.sub.h,t.sub.2)I(iu.sub.ijt.sub.h,jv.sub.ijt.sub.h,t.sub.1)=0, wherein: I is an intensity of the image data, i is a pixel index in a horizontal direction x, j is a pixel index in a vertical direction y orthogonal to the horizontal direction x, u.sub.ij is a component of the velocity vector on a pixel point, v.sub.ij is another component of the velocity vector on a pixel point, and t.sub.h is equal to t/2, wherein t is a time difference between the pair of image frames (t.sub.1 and t.sub.2).

3. The method of claim 1, wherein solving the symmetric displaced frame difference equation comprises using at least one iteration equation derived from the symmetric displaced frame difference equation by conversion to a fully or over-constrained system using a nonlinear least squares model of the displacement field.

4. A motion estimator apparatus, comprising: at least one processor; and a memory storing a symmetric displaced frame difference equation, wherein the at least one processor is operative to: receive an input image sequence comprising a pair of image frames individually including multidimensional image data corresponding to a plurality of pixel locations at different times, solve the symmetric displaced frame difference equation using; an iteration equation and the image data of the pair of image frames to determine a displacement field describing displacement vectors at half of the total displacement vector, and at least one bilinear polynomial function expressing a multidimensional displacement field, and determine displacement by motion for the input image sequence based on the solved symmetric displaced frame difference equation.

5. The motion estimator apparatus of claim 4, wherein the multidimensional image data of the pair of frames is two-dimensional, and wherein the symmetric displaced frame difference equation is
I(i+u.sub.ijt.sub.h,j+v.sub.ijt.sub.h,t.sub.2)I(iu.sub.ijt.sub.h,jv.sub.ijt.sub.h,t.sub.1)=0, wherein: I is an intensity of the image data, i is a pixel index in a horizontal direction x, j is a pixel index in a vertical direction y orthogonal to the horizontal direction x, u.sub.ij is a component of the velocity vector on a pixel point, v.sub.ij is another component of the velocity vector on a pixel point, and t.sub.h is equal to t/2, wherein t is a time difference between the pair of image frames (t.sub.1 and t.sub.2).

6. The motion estimator apparatus of claim 4, wherein the at least one processor is operative to solve the symmetric displaced frame difference equation using at least one iteration equation derived from the symmetric displaced frame difference equation by conversion to a fully or over-constrained system using a nonlinear least squares model of the displacement field.

7. A non-transitory computer readable medium with computer executable instructions for: providing a symmetric displaced frame difference equation; receiving an input image sequence comprising a pair of image frames individually including multidimensional image data corresponding to a plurality of pixel locations at different times; solving the symmetric displaced frame difference equation using: an iteration equation and the image data of the pair of image frames to determine a displacement field describing displacement vectors at half of the total displacement vector, and at least one bilinear polynomial function expressing a multidimensional displacement field; and determining displacement by motion for the input image sequence, based on the solved symmetric displaced frame difference equation.

8. The non-transitory computer readable medium of claim 7: wherein the multidimensional image data of the pair of frames is two-dimensional, and wherein the symmetric displaced frame difference equation is
I(i+u.sub.ijt.sub.h,j+v.sub.ijt.sub.h,t.sub.2)I(iu.sub.ijt.sub.h,jv.sub.ijt.sub.h,t.sub.1)=0, wherein: I is an intensity of the image data, i is a pixel index in a horizontal direction x, j is a pixel index in a vertical direction y orthogonal to the horizontal direction x, u.sub.ij is a component of the velocity vector on a pixel point, v.sub.ij is another component of the velocity vector on a pixel point, and t.sub.h is equal to t/2, wherein t is a time difference between the pair of image frames (t.sub.1 and t.sub.2).

9. The non-transitory computer readable medium of claim 7, comprising computer-executable instructions for solving the symmetric displaced frame difference equation using at least one iteration equation derived from the symmetric displaced frame difference equation by conversion to a fully or over-constrained system using a nonlinear least squares model of the displacement field.

10. The method of claim 1, wherein receiving the input image sequence comprises receiving only two image frames.

11. The motion estimator apparatus of claim 4, wherein the input image sequence comprises only two image frames.

12. The non-transitory computer readable medium of claim 7, wherein receiving the input image sequence comprises receiving only two image frames.

13. The method of claim 1, further comprising: determining, using the at least one processor, displacement by motion for the input image sequence based on the solved symmetric displaced frame difference equation.

14. The method of claim 1, wherein solving the symmetric displaced frame difference equation using the iteration equation and the image data of the pair of image frames to determine a displacement field describing displacement vectors comprises: determining, using the at least one processor, a displacement vector corresponding to a moving particle in the image sequence.

15. The method of claim 14, further comprising: determining, using the at least one processor and the determined displacement vector corresponding to the moving particle, a plot of a motion particle trajectory of the moving particle.

16. The motion estimator apparatus of claim 4, wherein the at least one processor is further operative to: determine displacement by motion for the input image sequence based on the solved symmetric displaced frame difference equation.

17. The non-transitory computer readable medium of claim 7, further comprising computer-executable instructions for: determining, using the at least one processor, displacement by motion for the input image sequence based on the solved symmetric displaced frame difference equation.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The following description and drawings set forth certain illustrative implementations of the disclosure in detail, which are indicative of several exemplary ways in which the various principles of the disclosure may be carried out. The illustrated examples, however, are not exhaustive of the many possible embodiments of the disclosure. Other objects, advantages and novel features of the disclosure will be set forth in the following detailed description of the disclosure when considered in conjunction with the drawings, in which:

(2) FIGS. 1 and 2 are schematic diagrams illustrating an input video image sequence which includes first and second image frames and corresponding pixel data obtained at different times;

(3) FIG. 3 is a schematic diagram illustrating an exemplary processor-implemented motion estimator apparatus using a fully constrained set of nonlinear equations to determine a displacement vector field from the input image sequence of FIGS. 1 and 2;

(4) FIG. 4 is a schematic diagram illustrating an exemplary fully constrained set of nonlinear equations including forward and backward displaced frame difference equations and two forward displacement vector invariant equations in the motion estimator apparatus of FIG. 3;

(5) FIG. 5 is a schematic diagram illustrating another exemplary fully constrained set of nonlinear equations including forward and backward displaced frame difference equations and two backward displacement vector invariant equations in the motion estimator apparatus of FIG. 3;

(6) FIGS. 6A and 6B are flow diagrams illustrating an exemplary method for processing an image sequence having two image frames to estimate a motion vector field in accordance with the present disclosure;

(7) FIG. 7A is a vector diagram illustrating an exemplary forward displacement vector;

(8) FIG. 7B is a vector diagram illustrating a corresponding backward displacement vector;

(9) FIG. 7C is a vector diagram illustrating exemplary displacement vectors corresponding to a pair of moving particles in an image sequence;

(10) FIGS. 8 and 9 are schematic diagrams illustrating a portion of an exemplary image frame partitioned into blocks;

(11) FIG. 10 is a flow diagram illustrating an exemplary adaptive framework for solving a system of equations in accordance with further aspects of the present disclosure; and

(12) FIG. 11 is a detailed flow diagram illustrating an adaptive process or framework for solving a system of equations for motion estimation in accordance with further aspects of the present disclosure.

(13) FIG. 12 is a plot of a motion particle trajectory, instant velocities at time t.sub.1 and t.sub.2, and a displacement vector between initial and final positions at time t.sub.1 and t.sub.2.

(14) FIG. 13 is a coarser node point in grid.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

(15) One or more embodiments or implementations are hereinafter described in conjunction with the drawings, where like reference numerals refer to like elements throughout, and where the various features are not necessarily drawn to scale.

(16) FIG. 1 illustrates an input image sequence 10 with a plurality of image frames 12 including a first image frame 12.sub.1, a final frame 12.sub.n, and one or more intervening frames 12 wherein the input sequence 10 can include any number n frames 12, where n is an integer greater than 1. As seen in FIG. 2, each individual image frame 12.sub.i includes an integer number of pixel data values, each corresponding to a two-dimensional pixel location in an X-Y plane, where each pixel location is indicated with an integer subscript index i indicating an X direction pixel index (1iI) and an integer subscript index j indicating a Y direction pixel index (1jJ). In the illustrated example, moreover, the image data represents optical (e.g., visible) data, such as video frames 12, where the input sequence 10 shows a vehicle artifact 14 moving towards a lower right corner of the image frame as the sequence proceeds from an initial (e.g., first) time t.sub.1 (frame 12.sub.1) to a final (e.g., second) time t.sub.2 (frame 12.sub.n). In practice, as is known, the individual image data for a given pixel location i,j can be a single value representing light intensity (luminance) and/or an individual data point may include luminance as well as color or chrominance data values (e.g., L*a*b*, RGB, CMYK data, etc.). In other embodiments, the image sequence 10 can include multidimensional image data representing thermal images, for instance, where the data values represent temperature, or other image sequences can be used in which the image intensity is recorded from physical sensors, for tracer, heat, or optical flow in space and time.

(17) FIG. 3 illustrates an exemplary motion estimator apparatus 100 that can perform various image processing tasks including determining or estimating a motion or displacement field 140 from an input image sequence 10 having a first image frame 12.sub.t1 including multidimensional image data corresponding to a plurality of pixel locations i, j at a first time t.sub.1, and a second image frame 12.sub.t2 including multidimensional image data corresponding to a plurality of pixel locations at a second time t.sub.2. The apparatus 100 can be used to determine a displacement field 140 according to the disclosed techniques for any two image frames 12, which may but need not be the initial and final frames 12 of an input image sequence 10. The apparatus 100 includes one or more processors 110 operatively coupled with a memory 120, and the apparatus 100 can be implemented in a single processing device or in distributed fashion across two or more processing elements 110. In the illustrated embodiment of FIG. 3, the processor 110 implements a PROC equation solver 110a used for solving iteration equations 132 with a Progressive Relaxation of the Over-Constraint (PROC) algorithm, as well as a Motion Vector Solver 110b is designed for solving the motion vectors at time t=t.sub.2 with a bilinear motion vector function 134 during an iteration. In addition, the motion estimator apparatus 100 may form part of a larger system, such as a video processor coupled to a transmitter, a video camera, a frame rate up-converter, etc.

(18) The estimator apparatus 100 receives the input image sequence 10 and generates a displacement vector field 140 which can be stored in the internal memory 120 and/or maybe outputted by the apparatus 100 alone or as part of a processed image sequence 200. In addition, the estimator provides an equation system 130, which may be stored in the electronic memory 120. The illustrated estimator 100 further includes at least one iteration equation 132 and the bilinear displacement (or motion) vector function 134, which can be stored in the memory 120 or otherwise be accessible for use by the processor 110 in performing the displacement field estimation function set forth herein. In particular, the iteration equations 132 in certain embodiments are derived from the equation set 130 by conversion to a fully or over-constrained system using a nonlinear least squares model of the displacement field, as discussed further below. In addition, the bilinear motion vector function 134 in certain embodiments expresses a multidimensional displacement field. The bilinear displacement vector function 134 in certain embodiments can be represented by compact form shown in equation (10) below.

(19) The equation system 130 is a fully constrained nonlinear equation set of equations, where four exemplary equations are illustrated in the example of FIG. 3. In this regard, the derived or estimated displacement field 140 describes displacement vectors r.sub.ij(t.sub.1) or r.sub.ij(t.sub.2) at pixel locations i,j at one of the first and second times t.sub.1 or t.sub.2, for example a first displacement vector r.sub.ij(t.sub.1)=(x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)) at t.sub.1, or a second displacement vector r.sub.ij(t.sub.2)=(x.sub.ij(t.sub.2), y.sub.ij(t.sub.2)) at t.sub.2. In the illustrated two dimensional case, the equation set 130 includes a forward displaced frame difference equation (FDFD.sub.ij), a backward displaced frame difference equation (BDFD.sub.ij), and two displacement vector invariant equations (DVI.sub.ij), which are shown as equations (7) and (8) below. In one embodiment shown in FIG. 4, the third and fourth equations of the set 130 are forward displacement vector invariant equations (FDVI), and in another embodiment shown in FIG. 5, the third and fourth equations are backward DVI equations (BDVI). Using the image data from the two received input frames 12.sub.t1 and 12.sub.t2 corresponding to times t.sub.1 and t.sub.2, the processor 110 is programmed to solve the equation set 130 to determine a displacement field (140) describing displacement vectors r.sub.ij(t.sub.1) or r.sub.ij(t.sub.2) at pixel locations i, j at t.sub.1 or t.sub.2 using initial conditions 136 including presetting displacement field values at the other time. In this regard, the set of four equations 130 has a total of four unknown variables for each pixel location i, j, namely x.sub.ij(t.sub.1), y.sub.ij(t.sub.1), x.sub.ij(t.sub.2), and y.sub.ij(t.sub.2), and therefore the equation set 130 is fully constrained without approximation and without imposition of any additional constraints or assumptions.

(20) In certain embodiments, the PROC equation Solver 110a is programmed to solve the equation set 130 using an iterative numerical and PROC techniques to determine the displacement field 140, and may employ any suitable initial conditions and loop termination logic, including without limitation a maximum number of iterations per pixel location i, j, alone or in combination with termination based on computed value changes being less than a predetermined threshold value. In certain embodiments, the Motion Vector Solver 110b solves the equations (3) and (4) below using a damped Newton-Raphson method with suitable initial values used in the computation. In other embodiments, the Motion Vector Solver 110b solves the equations (3) and (4) using bilinear modeling of the displacement field 140. The estimator 100 may provide the derived displacement field 140 for use in a variety of applications, such as video processing using an interpolator to construct one or more additional frames for frame rate up-conversion. In another example, the estimator 100 may provide the displacement field 140 for use with compression processing in a video encoder for selectively dropping certain frames 12 received with the input image sequence 10.

(21) FIGS. 6A and 6B illustrate an exemplary method 300 for processing image data from an input image sequence 10 in accordance with further aspects of the present disclosure. While the method 300 is illustrated and described below in the form of a series of acts or events, it will be appreciated that the various methods of the disclosure are not limited by the illustrated ordering of such acts or events. In this regard, except as specifically provided hereinafter, some acts or events may occur in different order and/or concurrently with other acts or events apart from those illustrated and described herein in accordance with the disclosure. It is further noted that not all illustrated steps may be required to implement a process or method in accordance with the present disclosure, and one or more such illustrated acts may be combined or omitted. The illustrated method 300 may be implemented in hardware, processor-executed software, or combinations thereof, such as in the exemplary motion estimator apparatus 100 described above, and may be embodied in the form of computer executable instructions stored in a computer readable medium, such as in a memory 120 operatively associated with the processor 110 in one example.

(22) At 302 in FIG. 6A, an input image sequence 10 is received including two frames, where the frame 12.sub.t1 includes a plurality of data values at time t.sub.1 I(i, j, t.sub.1), and the second frame 12.sub.t2 includes data values I(i, j, t.sub.2) at time t.sub.2 (e.g., image intensity data I, which itself can be multidimensional (e.g., CMYK, etc.), and can represent optical, thermal, or other image data types, where the input frames 12 can be stored in the estimator memory 120 in certain embodiments. The method 300 also includes presetting initial values of the motion field at an initial reference time (e.g., t.sub.1) at 304. In one possible implementation illustrated in FIG. 6A, seed values for the velocity vector components u.sub.ij(t.sub.1) and v.sub.ij(t.sub.1) are initialized at 304 to be 0.01 for each pixel location i, j, although this is one non-limiting example. A block size integer value (n.sub.0) is preset at 306, which represents an integer number of pixels within an analysis block or tile of the images 12, and in certain embodiments sets the values for n.sub.x and n.sub.y to indicate the number of interpolation points on the orthogonal x and y directions of the image frames 12. In the example of FIG. 6A, any integer value can be used for presetting the block size n.sub.0, such as 5 in one implementation. At 306, moreover, an initial iteration block size n is preset to an integer value that is greater than n.sub.0 (e.g., n is initially set to be n.sub.0+10 in one example). As discussed below, certain embodiments of the displacement field estimator apparatus 100 implements a progressive relaxation algorithm in which the block size n is selectively reduced during iteration to progressively relax an amount of over-constraint of the equation set 130. At 308, a loop count parameter Z is preset to an integer value, such as 10 in one example. In addition, an initial Levenburg-Marquardt factor is preset at 310 for the iteration processing. In one example, this factor is preset to 0.001, although other values can be used.

(23) At 320 in FIG. 6A, a fully constrained nonlinear equation set is provided (e.g., equation set 130 in FIGS. 3-5 above) that includes a forward displaced frame difference equation (e.g., FDFD.sub.ij), a backward displaced frame difference equation (e.g., BDFD.sub.ij), and two or more displacement vector invariant equations (e.g., DVI.sub.ij), which can be forward DVI equations (FIG. 4) or backward DVI equations (FIG. 5). The equation set is then solved for the displacement field (e.g., displacement field 140) at one of the first and second times (t.sub.1, t.sub.2) using known or estimated displacement field values at the other of the first and second times (t.sub.1, t.sub.2). Any suitable numerical technique can be employed at 330 to solve the equation set 130, such as Damped Newton-Raphson or other iterative technique, using bilinear modeling of the displacement field 134, etc., wherein any suitable initial values and termination conditions may be employed.

(24) In the illustrated embodiment, the displacement vectors on node points are solved at 130 (FIG. 6A) by the iteration equations 132 and the displacement field estimate r.sub.ij(t.sub.1) is updated and r.sub.ij(t.sub.2) is solved in the local domain .sub.k1 (FIGS. 8 and 9 below), after which the process 300 proceeds to 360 in FIG. 6B where a cost function is evaluated. A determination is made at 340 (FIG. 6B) as to whether the cost function is minimized or another termination condition is satisfied. If so (YES at 340), a determination is made at 350 as to whether the block size n is equal to the final block size value n.sub.0 preset at 304 in FIG. 6A. If so (YES at 350), the motion field estimation is complete at 360.

(25) If the cost function is not minimized and no other termination conditions are satisfied (NO at 340 in FIG. 6B), the Levenburg-Marquardt factor is multiplied by 10 for the next iteration at 342, and the process 300 returns to again solve the displacement vectors at 320 and the iteration equations at 330 in FIG. 6A. This iterative algorithm continues as described above until the cost function is minimized or another termination condition is satisfied (YES at 340) and the block size n has not yet reached the preset value of n.sub.0 (NO at 350). In this case, a determination is made at 352 as to whether the preset value Z iterations have occurred since the last adjustment to the block size parameter n. If so (YES at 352), the block size n is decremented at 354, and the displacement field estimate r.sub.ij(t.sub.1) is updated at 356, after which the process 300 returns to again solve the iteration equations at 330 in FIG. 6A. Otherwise (NO at 352 and FIG. 6B), the current block size value n is maintained, the displacement field estimate r.sub.ij(t.sub.1) is updated at 356, and the process 300 returns to 330 in FIG. 6A for further iterations. In this manner, the process 300 provides an iterative approach for estimating the r.sub.ij(t.sub.1)={x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)} at the reference time t.sub.1 using the bidirectional DFD and DVI equations and the iteration equations 132 in the system 130.

(26) As noted above, the motion estimator apparatus 100 employs a fully constrained nonlinear system of equations 130 which includes forward and backward displaced frame difference equations (EQ1 and EQ2 in FIG. 3 above) as well as two or more displacement vector invariant equations (EQ3 and EQ4 in FIG. 3). In addition, the displacement vector invariant equations may be either forward DVI equations or backward DVI equations in certain embodiments.

(27) As further discussed below, the inventor has appreciated that such an equation set 130 facilitates motion determination from an image sequence in a variety of computer vision and remote sensing applications such as velocity or displacement estimation from motion, object tracking or recognition, medical imaging, advanced video editing, and ocean surface current estimation. The use of this fully constrained equation set 130 without having to make further approximations or to impose any additional constraints or assumptions provides a novel solution to the inverse problem of motion estimation in two successive image frames 12. In this regard, determination of the instant velocity is ill-posed because the information about the path and rate is lost after a temporal sampling. If the images have enough texture morphologies, both initial and final configurations of a moving particle are recorded by the image sequence. The inverse problem for determination of the displacement field 140, however, is well-posed because both initial and final positions can be determined and observed physically based on the input image sequence 10.

(28) The inventor has appreciated that if the initial and final positions of a particular artifact (e.g., artifact 14 in FIG. 2 above) can be physically observed and determined from an image sequence, a fully constrained system of equations 130 can be constructed for resolving the inverse problem. If the initial and final positions cannot be physically determined or observed, the fully constrained equation system 130 can still be used to solve the inverse problem, although there may exist multiple solutions of the fully constrained (nonlinear) system 130 that are consistent with the physical observation. In addition, all motion fields are defined by the displacement field 130 or the average (conservative) velocity fields denoted by v=v=r/t, where r is the displacement vector and t=t.sub.2t.sub.1 is the time difference between image frames at time t.sub.1 and t.sub.2.

(29) As used herein, conservative velocity is not intended to mean that the physical velocity is a conservative quantity. This concept of equivalent motion fields between the displacement and the conservative velocity is based on the definition. If the initial and final positions of a moving particle or artifact in two successive frames are the only input information, the unique total displacement vector may be associated with many different intermediate dynamic motion states. However, only the motion with the conservative velocity is collinear with the displacement trajectory (a straight line). Focusing strictly on the initial and final results, the conservative velocity as used herein is the resolved velocity field, i.e., the solved velocity is only an average (conservative) velocity. To avoid any misunderstanding on this concept, all derivations in this disclosure are based on the displacement vector field 140 that describes displacement vectors (r.sub.ij(t.sub.1), r.sub.ij(t.sub.2)) with respect to pixel locations (i, j) at the first and second times (t.sub.1, t.sub.2). In this respect, it is understood that there is only limited information available from an input image sequence 10 to study the inverse problem without any way of absolutely knowing the real physical (dynamic) processes or the identity of physical objects, which can be rigid bodies, liquid flow, or deformed objects in an image scene 12, but the above described processed 300 can be advantageously employed to determine motion fields that are consistent with the physical observation.

(30) To illustrate, a function I(r, t) is defined as a scalar property (e.g., intensity) of a continuum of particles that may be expressed as a point function of the coordinates (r=r(t)), where r is a multidimensional position (e.g., x and y positions in a two dimensional example). In other examples, depending on the application, I(r, t) can be defined as the intensity of the optical flow (in computer vision), the tracer concentration (ocean color), or the temperature of heat flow (in geosciences and remote sensing). If the intensity between two images is conserved, then the following equation (1) represents the conservation constraint, regardless of a source term:

(31) $\begin{array}{cc}\frac{I}{t}=\left(\frac{}{t}+v.Math.\right)I=0,& \left(1\right)\end{array}$
where the operator

(32) $\frac{}{t}$
denotes a total derivative with respect to time t, and

(33) $v=\frac{r}{t}=v\left(r\left(t\right),t\right)$
is the velocity vector in Cartesian coordinates. Equation (1) is sometimes referred to as the optical flow or brightness conservation constraint equation (in computer vision), or the heat flow or tracer conservation equation (in geophysics), and is a differential form of conservation constraint. The differential form of the conservation constrained equation (1) which contains linear terms of the components of the velocity holds only for infinitesimal or close infinitesimal motions.
In order to constrain the image scenes at time t=t.sub.1 and t=t.sub.2, equation (1) is an integrated from time t.sub.1 to t.sub.2 as follows:

(34) $\begin{array}{cc}{}_{t_1}^{t_2}\frac{I\left(r\left(t\right),t\right)}{t}t=I\left(t\left(t_2\right),t_2\right)-I\left(r\left(t_1\right),t_1\right)0,& \left(2\right)\end{array}$
where r(t.sub.1) and r(t.sub.2) are the position vectors at time t.sub.1 and t.sub.2. If a displacement vector field is defined by r=r(t.sub.2)r(t.sub.1), then the displaced frame difference (DFD) equation is given by DFD=I(r(t.sub.1)+r, t.sub.2)I(r(t.sub.1),t.sub.1)0. The conservation constraint equation (1) is thus temporally integrated into the path independent equation (2). It is noted that although the DFD equation is an implicit function on the displacement field 140, the two path independent terms in equation (2) correspond to the initial and final states of motion associated with the two successive frames for this conservation system. Employing the DFD equation derived by integrating between time of two successive frames can achieve higher accuracy in comparison with the differential form of the conservation constraint optical flow equation (1) (a first order term of Taylor expansion of the DFD), especially for large scale displacement motion estimation. As the image intensity change at a point due to motion has only one constraint equation (1), while the motion vector at the same point has two components (a projection case), the motion field was previously believed to be not computable without an additional constraint.

(35) Referring also to FIGS. 7A-7C, the displacement vector r(r(t.sub.1),t.sub.1) (forward displacement vector) connects the initial and final positions of the particle 402 as shown in the plot 400 of FIG. 7A and the velocity field is the time rate of change of the property with respect to specific particles of the moving continuum. The same displacement vector on an image scene at time t.sub.1 and t.sub.2 may have two different forms. We define these displacement vectors r(r(t.sub.1),t.sub.1) and r(r(.sub.t2),t.sub.2) as shown in the graph 404 10 of FIGS. 7A and 7B as forward and backward displacement vectors, respectively. FIG. 7C shows forward and backward displacement vectors at a fixed position r at times t.sub.1 and t.sub.2 for two particles (square and round in the figure), where the square particle begins at position 422 (at t.sub.1) and ends at position 424 (at t.sub.2), and the round particle begins at position 424 (at t.sub.1) and ends at position 426 (at t.sub.2) in FIG. 7C. If the displacement vector r(t) is defined by a forward or backward vector on image frames 12 at times t.sub.1 and t.sub.2, the Displacement Vector Invariant (DVI) equation is given by
r(r(t.sub.2),t.sub.2)+r(r(t.sub.2)+r(r(t.sub.2),t.sub.2),t.sub.1)=0,(3)
or
r(r(t.sub.1),t.sub.1)+r(r(t.sub.1)+r(r(t.sub.1),t.sub.2),t.sub.2)=0,
and the Conservative Velocity Constraint (CVC) equation is given by
v(r(t.sub.1)+v(r(t.sub.1),t.sub.1)t,t.sub.2)=v(r(t.sub.1),t.sub.1).(4)
or
v(r(t.sub.2),t.sub.2)=v(r(t.sub.1)v(r(t.sub.2),t.sub.2)t,t.sub.1),
where t=t.sub.2t.sub.1.

(36) The DVI equation (3) establishes an implicit function relationship between the forward and backward displacement fields. Both displacement and conservation velocity fields are equivalent based on the definition. According to the definition of the forward and backward displacement vectors r(r(t.sub.1),t.sub.1) and r(r(.sub.t2),t.sub.2) in FIGS. 7A and 7B, we have the following equations (5):

(37) $\begin{array}{cc}r\left(r\left(t_1\right),t_1\right)=-r\left(r\left(t_2\right),t_2\right)\{\begin{array}{c}r\left(r\left(t_1\right),t_1\right)=r\left(t_2\right)-r\left(t_1\right)=v\left(r\left(t_1\right),t_1\right)\left(t_2-t_1\right)\\ r\left(r\left(t_2\right),t_2\right)=r\left(t_1\right)-r\left(t_2\right)=v\left(r\left(t_2\right),t_2\right)\left(t_1-t_2\right).\end{array}& \left(5\right)\end{array}$

(38) Equations (5) indicate that both the displacement and conservative vectors are equivalent because each vector can be obtained by multiplying and dividing a time different factor (t.sub.2t.sub.1) or (t.sub.1t.sub.2). The DVI equation (3) indicates that the motion fields at time t.sub.1 and t.sub.2 are not equal at the same position, but both fields have a shift from each other for moving objects. The shift vector is the displacement vector of the time differences. Equation (3) or (4) establishes an implicit recursive relationship at time t.sub.1 and t.sub.2 for the displacement vector fields or conservative velocity fields. The inventor has appreciated that since the DVI equation (3) or the CVC equation (4) is a vector equation, the total number of component equations is equal to the number of dimensions of the velocity, and as a result, if one field is given or solved then another corresponding field can be determined by equation (3) or (4) completely.

(39) A fully constrained system of equations 130 can thus be provided in the motion estimator apparatus 100 using the motion compensation concept. The temporal integral form of the conservation constraint equation (2) or the DFD equation indicates that the image at time t.sub.1 can be predicted by motion-compensated prediction with an image at time t.sub.2 and the displacement field 140 at time t.sub.1, thus facilitating interpolation or extrapolation for video coding and/or frame rate up-conversion applications. Also, the image at time t.sub.2 can be predicted also with an image at time t.sub.1 and the displacement field 140 at time t.sub.2. The inventor has thus found that the equation (2) can be described forward and backward at a fixed position based on the displacement fields 140 at different times, t.sub.1 and t.sub.2.

(40) Assuming that the number of pixels in an image frame 12 is equal to N, the total number of data points is equal to 2N for an input image sequence 10 having two frames (e.g., frames 12.sub.t1 and 12.sub.t2 in FIG. 3 above). A fully constraint nonlinear system of equations 130 using forward and backward DFD equations combined with the DVI equations (4) for resolving the displacement field is thus given by the following equation (6):

(41) $\begin{array}{cc}\{\begin{array}{c}\mathrm{FDFD}\left(r,r\left(r,t_1\right)\right)=I\left(r+r\left(r,t_1\right),t_2\right)-I\left(r,t_1\right)=0\\ \mathrm{BDFD}\left(r,r\left(r,t_2\right)\right)=I\left(r,t_2\right)-I\left(r+r\left(r,t_2\right),t_1\right)=0\\ \mathrm{DVI}=r\left(r,t_2\right)+r\left(r+r\left(r,t_2\right),t_1\right)=0,\end{array}& \left(6\right)\end{array}$
where r(r, t.sub.1) and r(r, t.sub.2) are the forward and backward displacement fields 140 at a fixed position r as shown in FIG. 7C. Both the FDFD and BDFD equations at a fixed position r are independent for the displacement fields, r(r, t.sub.1) and r(r, t.sub.2). The total number of two unknown components of r(r, t.sub.1) or r(r, t.sub.2) is equal to 2N. If one of the fields, r(r, t.sub.1) and r(r, t.sub.2), is chosen as an independent motion field, another field can be solved by the DVI equation. The total unknown variables 2N in equations (6) is equal to the total number of constraint equations (6) and as well as the total number of data points in an image sequence. Therefore, both data points and equations match and provide sufficient data information and constraints to solve the two components of the displacement vector 140 without resort to any additional sum assumptions or constraints.

(42) The forward and backward displacement vectors on two image frames 12.sub.t1 and 12.sub.t2 at times t.sub.1 and t.sub.2 as shown in FIGS. 7A and 7B are defined by equations (5), and substituting the position vectors r(t.sub.2) and r(t.sub.1) from equations (5) into equation (2) yields the following two equations:

(43) $\{\begin{array}{c}I\left(r\left(t_1\right)+r\left(r\left(t_1\right),t_1\right),t_2\right)-I\left(r\left(t_1\right),t_1\right)=0\\ I\left(r\left(t_2\right),t_2\right)-I\left(r\left(t_2\right)+r\left(r\left(t_2\right),t_2\right),t_1\right)=0.\end{array}$
Since the above two equations hold for all image scenes, the equations (6) result if all position vectors in above equations are at a fixed position r, and the first and second equations are denoted by FDFD and BDFD. The solution to the inverse problem is therefore based on the motion analysis in physics, the DVI equation, and recognized different displacement fields at time t.sub.1 and t.sub.2 in both FDFD and BDFD equations. If the forward displacement field is chosen as independent variables, the DVI equations (whether forward DVI equations or backward DVI equations are used) link the FDFD and BDFD equations together for solving the forward two component displacement field.

(44) Using the conservative velocity and CVC equation to replace the displacement vector and DVI equation in (6), a set of fully constrained nonlinear system of equations 130 can be determined for solving the conservative velocity field. Defining a function with discrete variables i and j as f.sub.ij(t)=f (i, j, t), the FDFD, BDFD, and the component DVI equations in (6) on a given pixel point i, j are given by the following equations (7) and (8) as the equation set 130:

(45) $\begin{array}{cc}\{\begin{array}{c}{\mathrm{FDFD}}_{\mathrm{ij}}=I\left(+x_{\mathrm{ij}}\left(t_1\right),j+y_{\mathrm{ij}}\left(t_1\right),t_2\right)-I_{\mathrm{ij}}\left(t_1\right)=0\\ {\mathrm{BDFD}}_{\mathrm{ij}}=I_{\mathrm{ij}}\left(t_2\right)-I\left(+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0,\end{array}\mathrm{and}& \left(7\right)\\ \{\begin{array}{c}{x}_{\mathrm{ij}}\left(t_2\right)+x\left(+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0\\ y_{\mathrm{ij}}\left(t_2\right)+y\left(+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0,\end{array}& \left(8\right)\end{array}$
where r.sub.ij(t)=(x.sub.ij(t), y.sub.ij(t)).sup.T. The displacement vector field at time t.sub.2 is a function of the field at time t.sub.1 based on the DVI equation (8), where the same is true of the backward DVI equations discussed below. There are several available techniques for solving nonlinear systems of equations 130. For example, a damped Newton-Raphson method is an efficient technique, but a converged solution may be facilitated by having a good guess for the initial values, especially for a huge dimensional problem.

(46) Referring also to FIG. 10, a unified adaptive framework 600 may be used with simplicity and efficiency for different applications, using techniques to adjust the degree of over-constraint from fully constraint to over-constraint by a block size parameter for solving the nonlinear problem 130. In certain embodiments, the velocity field is expanded by bilinear polynomial functions or two-dimensional B-Spline functions. In order to develop an adaptive framework that allows the estimator to control the resolutions of the displacement field and the degree of over-constraint, and to simplify the computational complexity, bilinear modeling of the displacement field can be used. In general, any two-dimensional function can be approximated by a Lagrange's bilinear function (9):

(47) $\begin{array}{cc}f\left(x,y\right)=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}f\left(p+n_x,q+n_y\right)H_{p+n_x,q+n_y}\left(x,y\right),& \left(9\right)\end{array}$
where the function H.sub.a,b (x, y) is defined by:

(48) 0$H_{a,b}\left(x,y\right)=\frac{1}{n_x{n}_y}\{\begin{array}{cc}\left(n_x-x+p\right)\left(n_y-y+q\right)& \left(a=p.Math.b=q\right)\\ \left(x-p\right)\left(n_y-y+q\right)& \left(a=p+n_x.Math.b=q\right)\\ \left(n_x-x+p\right)\left(y-q\right)& \left(a=p.Math.b=q+n_y\right)\\ \left(x-p\right)\left(y-q\right)& \left(a=p+n_x.Math.b=q+n_y\right),\end{array}$
and where the parameters of block size n.sub.x and n.sub.y are the sampled spaces of the function f on x and y directions as shown in FIG. 5. The quantized indices p and q are functions of x and y and are given by:

(49) $\left\{p,q\right\}=\left\{p\left(x\right),q\left(y\right)\right\}=\left\{n_x.Math.\frac{x}{n_x}.Math.,n_y.Math.\frac{y}{n_y}.Math.\right\},$
where denotes an integer operator. The {p, q} serve as tile indices because the integer operator increments them by unity after an additional n.sub.x or n.sub.y pixels are counted.

(50) Referring also to FIGS. 8 and 9, in this approach, the image domain 500 is partitioned into a number of sub-domains, or tiles, each which contains an n.sub.xn.sub.y array of pixels. In the example of FIGS. 8 and 9, n=n.sub.x=n.sub.y=4 pixels per block. Node points (marked with solid squares) of the bilinear approximation are indexed with p and q. As seen in FIG. 9, the pixel locations are partitioned into blocks or tiles with node points 502 at the corners of the blocks (node points 502 indicated as crosshatched or white squares in FIGS. 8 and 9). During the iterations of the process 300, the current value of the block size parameter n sets the number of interpolation points n.sub.x and n.sub.y on the x and y directions, respectively, and FIG. 9 illustrates the blocks in a local domain proximate a particular node 502 in the image 500.

(51) The displacement field is modeled at 602 in the adaptive framework 600 of FIG. 10. Defining a function f.sub.ij(t)=f (i, j, t) with discrete variables i and j, the displacement fields r.sub.ij(t.sub.1) are represented as two-dimensional bilinear functions controlled by a smaller number of displacement vector estimates r.sub.pq(t.sub.1) which lie on a coarser node point in grid. The two displacement components in an image scene can be expressed by the following bilinear polynomials:

(52) $\begin{array}{cc}{r}_{\mathrm{ij}}=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}{r}_{p+n_x,q+n_y}{H}_{p+n_x,q+n_y}\left(i,j\right),& \left(10\right)\end{array}$
where the displacement vector r.sub.ij=r.sub.ij(t.sub.1)=(x.sub.ij, y.sub.ij). In the special case when block size is unity (n.sub.x=n.sub.y=1), r.sub.ij=r.sub.pq for all indices i and j. All displacement vectors r.sub.ij can be calculated with the bilinear function using the displacement on node points expressed as r.sub.pq. Displacement vectors r.sub.ij off-node are no longer independent variables (for over-constrained case: n=n.sub.x=n.sub.y>1) except on node points (or a fully constrained case: n=1).

(53) In this regard, the block (tile) size parameter n1 can be adjusted to control the number of interpolation points related to the resolution of the displacement field and the degree of the over-constraint. When the parameter n is equal to one, all node points and pixel positions are overlapped together and the system is fully constrained. The system is over-constrained if the parameter n is greater than one.

(54) A nonlinear least-squares model can also be used at 604 in FIG. 10. The displacement field x.sub.pq and y.sub.pq can be resolved when the parameters n1 from the fully or over-constraint system using the nonlinear least-squares model. The presence of possible quantization errors and noise in the intensity measurements suggest the FDFD and BDFD equations in (8) are never identically zero. However, a set of r.sub.ij can be selected for which it is minimized in a least-squares sense. A Mean Square Error (MSE) function that is a sum of the total errors of the FDFD and BDFD in (8) is given by:

(55) $\mathrm{MSE}=\frac{1}{N}\underset{i,j}{.Math.}\left({\mathrm{FDFD}}_{\mathrm{ij}}^2+{\mathrm{BDFD}}_{\mathrm{ij}}^2\right),$
where the range of i and j are the entire (N=N.sub.xN.sub.y) image pixel domain (i[0, N.sub.x1], j[0, N.sub.y1], and the weighting factor can be equal to unity in this case). By minimizing the cost function with respect to the displacement components x.sub.kl and y.sub.kl as variables for given indices k and l on all node points, a fully or over-constrained system of equations (11) for the displacement may be written as follows:

(56) $\begin{array}{cc}\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\mathrm{FDFD}}_{\mathrm{ij}}\left\{\frac{{\mathrm{FDFD}}_{\mathrm{ij}}}{x_{\mathrm{kl}}},\frac{{\mathrm{FDFD}}_{\mathrm{ij}}}{y_{\mathrm{kl}}}\right\}+\underset{i,j{}_{k^{}{l}^{}}}{.Math.}{\mathrm{BDFD}}_{\mathrm{ij}}\left\{\frac{{\mathrm{BDFD}}_{\mathrm{ij}}}{x_{\mathrm{kl}}},\frac{\mathrm{BDFD}}{y_{\mathrm{kl}}}\right\}=0.& \left(11\right)\end{array}$

(57) In addition, the summation domain is reduced from the entire image plane to only a local region .sub.kl or .sub.kl, so that:

(58) $\underset{i,j}{.Math.}->\underset{i,j{}_{}}{.Math.}=\underset{i=-n_x+1}{\overset{+n_x-1}{.Math.}}\underset{j=-n_y+1}{\overset{+n_y-1}{.Math.}},$
where ={k, k} and ={1, 1}, and regions of the summation coverage k and l are defined by {k, l}={k+x.sub.k, I+y.sub.kl}. To obtain equation (14) below, the following is used:

(59) $\frac{x_{\mathrm{ij}}}{x_{\mathrm{kl}}}=\frac{y_{\mathrm{ij}}}{y_{\mathrm{kl}}}=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}{H}_{\mathrm{kl}}\left(i,j\right){}_{k,p+n_x}{}_{I,q+n_y},$
where .sub.ij is the Kronecker-Delta symbol.

(60) Equations (11) obtained by the nonlinear least-squares model are a set of nonlinear system equations with all displacement vectors x.sub.pq and y.sub.pq on node points as variables. To solve the nonlinear system equations (11), an iterative equation can be formulated at 606 in FIG. 10 based on Gauss-Newton and Levenberg-Marquardt methods. The displacement vectors x.sub.kl and y.sub.kl and their surrounding nearest neighbors are the only members now contained in the summations in equation (11). To solve these sparse nonlinear system of equations (11), the FDFD and BDFD can be expanded in a Taylor series using the Gauss-Newton method:

(61) $\left\{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m+1\right)},{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m+1\right)}\right\}\left\{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)},{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}\right\}++\left\{\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}},\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\right\}\left(x_{\mathrm{kl}}^{\left(m+1\right)}-x_{\mathrm{kl}}^{\left(m\right)}\right)++\left\{\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}},\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}}\right\}(y_{\mathrm{kl}}^{\left(m+1\right)}-y_{\mathrm{kl}}^{\left(m\right)},$
where m is a iteration index, and
{FDFD.sub.ij.sup.(m),BDFD.sub.ij.sup.(m)}={FDFR.sub.ij(x.sub.pq.sup.(m),y.sub.pq.sup.(m)),BDFD.sub.ij(x.sub.pq.sup.(m),y.sub.pq.sup.(m))}.

(62) Utilizing these expansions and the above equations (11), the following iterative equations can be used for all indices k and l:

(63) $\begin{array}{cc}\left(\begin{array}{c}{x}_{\mathrm{kl}}^{\left(m+1\right)}\\ y_{\mathrm{kl}}^{\left(m+1\right)}\end{array}\right)=\left(\begin{array}{c}{x}_{\mathrm{kl}}^{\left(m\right)}\\ y_{\mathrm{kl}}^{\left(m\right)}\end{array}\right)-{\left(\begin{array}{cc}\left(+1\right)a_{11}& a_{12}\\ a_{12}& \left(+1\right)a_{22}\end{array}\right)}^{-1}\left(\begin{array}{c}{b}_1\\ b_2\end{array}\right),\mathrm{where}\{\begin{array}{c}{a}_{11}=\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\left(\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\right)}^2+\underset{i,j{}_{\mathrm{kl}}^{\left(m\right)}}{.Math.}{\left(\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\right)}^2\\ a_{12}=\underset{i,j{}_{\mathrm{kl}}}{.Math.}\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}}+\underset{i,j{}_{\mathrm{kl}}^{\left(m\right)}}{.Math.}\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\\ a_{22}=\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\left(\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}}\right)}^2+\underset{i,j{}_{\mathrm{kl}}^{\left(m\right)}}{.Math.}{\left(\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}}\right)}^2,\end{array}\{\begin{array}{c}{b}_1=\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}+\underset{i,j{}_{\mathrm{kl}}^{\left(m\right)}}{.Math.}{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{x_{\mathrm{kl}}^{\left(m\right)}}\\ b_2=\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}\frac{{\mathrm{FDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}}+\underset{i,j{}_{\mathrm{kl}}^{\left(m\right)}}{.Math.}{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}\frac{{\mathrm{BDFD}}_{\mathrm{ij}}^{\left(m\right)}}{y_{\mathrm{kl}}^{\left(m\right)}},\end{array}\underset{i,j{}_{\mathrm{kl}}^{\left(m\right)}}{.Math.}=\underset{i=k^{\left(m\right)}-n_x+1}{\overset{k^{\left(m\right)}+n_x-1}{.Math.}}\underset{j=I^{\left(m\right)}-n_y+1}{\overset{I^{\left(m\right)}+n_y-1}{.Math.}},\mathrm{and}\left\{k^{\left(m\right)},I^{\left(m\right)}\right\}=\left\{.Math.k+x_{\mathrm{kl}}^{\left(m\right)}.Math.,.Math.I+y_{\mathrm{kl}}^{\left(m\right)}.Math.\right\}.& \left(12\right)\end{array}$

(64) The parameter 0 is a Levenberg-Marquardt factor that is adjusted at each iteration to guarantee that the MSE is convergent. A smaller value of the factor can be used, bringing the algorithm closer to the Gauss-Newton method with second order convergence. This Levenberg-Marquardt method can improve converge properties greatly in practice and has become the standard of nonlinear least-squares routines.

(65) All displacement vectors x.sub.pq and y.sub.pq on node points can be obtained by iterating the equations in (12), and the displacement vectors x.sub.ij and y.sub.ij off node points (for n>1) can be calculated by the bilinear functions in (10). Ultimately, all displacement vectors r can be determined, and an optimum solution can be achieved over the large-scale image. In addition, the block size parameter n1 can be adjusted to control the number of interpolation points within a tile, to resolve the displacement field, and to control the degree of the over-constraint.

(66) An optimized motion-compensated predication is also possible. In particular, using an additionally constrained system to estimate motion field, the inverse problem has previously been addressed by minimizing an objective function with a weighting (penalty) parameter. However, there are two major issues with this use of a weighting parameter. The first is the determination of the optimized weighting parameter, and several different values of the weighting parameter have been proposed. However, it is difficult to find a single optimal value of the parameter for realistic applications if the ground truth flow field is unknown. The second issue is that the Peak Signal-to-Noise Ratio (PSNR):

(67) $\mathrm{PSNR}=10{\log }_{10}\left(\frac{{255}^2}{\mathrm{MSE}}\right),$
is not optimized by minimizing the objective function with the weighting parameter. The estimated flow field by this approach cannot always lead to an optimal Motion-Compensated prediction or Interpolation (MCI) image for applications of video compression.

(68) As seen above, however, the iteration equations (12) are derived based on the least-squares principle that leads directly to a solution of the displacement field 140 with a minimized target function MSE or a maximized PSNR (An average PSNR for the Forward and Backward MCP (FMCP and BMCP)). Since both FDFD and BDFD are equivalent physical quantities, the target function PSNR is an optimized function without any additional parameters. Therefore the MCP image using the estimated displacement field 140 based on the fully constrained system 130 is believed to be optimized.

(69) The adjustable block size approach using smaller number of displacement components on nodes to interpolate a full density displacement field 140 (for n>1 case) provides a powerful capability for different applications in both computer vision and remote sensing fields. If the block size shown in FIG. 8 is nn with dimension N.sub.xN.sub.y images, then the number of transmitted or stored motion vectors for both proposed and block-matching is estimators are equal to [(N.sub.x1)/n+1][(N.sub.y1)/n+1] and (N.sub.x/n)(N.sub.y/n), respectively. Using almost the same number of the displacement vectors r in a fixed block size for both approaches, the current framework (the displacement field with C.sup.1 continuity obtained by global optimal strategies) can provide much higher accurate performance than the block-matching algorithm (the displacement field with C.sup.0 continuity obtained by local searching strategies) which is the currently adopted standard for video coding.

(70) With respect to iteration algorithms, certain embodiments of the process 300 and apparatus 100 can start from a set of preset initial values of the displacement field r.sup.(0) at time t.sub.1, then the Motion Vector Solver 110b solves the correspondence field r.sup.(0)(t.sub.2) using all the component equations (8) numerically. An iteration step based on these two displacement fields at time t.sub.1 and t.sub.2 and the iterative equations in (12) can be performed. Furthermore, employing a principle similar to that of the Gauss-Seidel iteration algorithm, updated values of x.sub.pq.sup.(m) and y.sub.pq.sup.(m) can be used on the right-hand side of equation (12) as soon as they become available during the iteration processing. In certain embodiments, all initial displacement field vectors are preset to be equal to 0.01, and the initial Levenberg-Marquardt factor is set to 0.001 and is multiplied by ten in each step if the iteration is not convergent.

(71) The FDFD and BDFD equations on each pixel include two motion-compensated predictions I(ix.sub.ij, jy.sub.ij, t.sub.{2,1}) with variables that may be out of the position on pixels in an image scene. In order to compute the motion-compensated predictions, the general bilinear interpolation function in (4) is utilized for this computation as follows:

(72) 0$I\left(i{x}_{\mathrm{ij}},j{y}_{\mathrm{ij}},t_{\left\{2,1\right\}}\right)=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}{I}_{p+,q+}\left(t_{\left\{2,1\right\}}\right)H_{p+,q+}\left(i{x}_{\mathrm{ij}},j{y}_{\mathrm{ij}}\right),$
where the function H is evaluated when n.sub.x=n.sub.y=1, and {p, q} {p(ix.sub.ij), q(j.sub.ij)}.

(73) The proposed displacement estimation approach leads to a nonlinear system of equations that may have multiple solutions depending on texture morphology in an image sequence. For example, to estimate a displacement field within a featureless region, the displacement field may not be unique because the initial and final positions of a particular particle cannot be physically determined. Even in a texture-rich environment, the realistic motion fields may have multiple possibilities, which satisfy the same equations and are consistent with the same physical observation. The multiple solutions in the inverse problem by solving a nonlinear system are congruent with this physical property.

(74) The inventor has further appreciated that in order to approach a globally minimized solution using the iteration equations (12), an algorithm of PROC that adapts a variable resolution of the displacement structure during the iterations can be employed in this algorithm. In certain embodiments, an initial block size parameter n.sub.0 is selected to be greater than a preset value of the block size n at initial iteration, and it reaches a preset value of n at the end iteration. In certain embodiments, the displacement field can be regularized by changing the block size parameter n from a larger value (higher degree of over-constraint) to a smaller one (lower degree of over-constraint) by one every Nth iteration until it approaches a preset value of n. The inventor has appreciated that the PROC algorithm is helpful for seeking a flow field in which each vector is consistent with its neighbors.

(75) As seen in FIGS. 4 and 5 above, either forward or backward DVI equations can be used in the equation set 130. Algorithms are then derived at 608 in FIG. 9 for solving the DVI equations. Two components of forward and backward DVI equations between time t.sub.1 and t.sub.2 are given by the following equations (A1) and (A2):

(76) $\begin{array}{cc}\{\begin{array}{c}{x}_{\mathrm{ij}}\left(t_2\right)+x\left(i+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0\\ y_{\mathrm{ij}}\left(t_2\right)+y\left(i+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0,\end{array}\mathrm{and}& \left(\mathrm{A1}\right)\\ \{\begin{array}{c}x\left(x,y,t_1\right)+x\left(x+x\left(x,y,t_1\right),y+y\left(x,y,t_1\right),t_2\right)=0\\ y\left(x,y,t_1\right)+y\left(x+x\left(x,y,t_1\right),y+y\left(x,y,t_1\right),t_2\right)=0.\end{array}& \left(A2\right)\end{array}$
If the conservative velocity field {x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)} at time t.sub.1 is given or solved, then the correspondence field {x.sub.ij(t.sub.2), y.sub.ij(t.sub.2)} at time t.sub.2 can be determined (or vice versa) by all component equations in (A1) or (A2). Thus, the processor 110 of the estimator 100 is programmed to solve the equation set 130 using the forward and backward displaced frame difference equations in combination with either the forward DVD equations (A1) (FIG. 4) or the backward DVD equations (A2) (FIG. 5) using the first and second image frames 12.sub.t1 and 12.sub.t2 to estimate the displacement field 140 describing displacement vectors r.sub.ij(t.sub.1) or r.sub.ij(t.sub.2) for pixel locations (i, j) at one of the first and second times (t.sub.1, t.sub.2) using known or estimated displacement field values at the other of the first and second times (t.sub.1, t.sub.2).

(77) In one example, numerically solving the equation set 130 at time t.sub.2 using equations (A1) or (A2) by the Motion Vector Solver 110b if the velocity {x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)} or {x.sub.ij(t.sub.2), y.sub.ij(t.sub.2)} are given involves expanding the field by a bilinear polynomial function, where a bilinear expression of a two-dimensional displacement field is given by equation (10).

(78) The forward and backward DVI equations (A1) and (A2) are implicit recursive functions of the fields ({x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)}) at time t and t.sub.1. Three methods for solving the matrix field {x.sub.ij(t.sub.2), y.sub.ij(t.sub.2)} are described below for the case in which the matrix field {x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)} are given, and it will be appreciated that the converse problem can be solved by similar techniques for solving the equation set 130 at time t.sub.1 where the matrix field is given at time t.sub.2.

(79) An interpolation method with a searching algorithm can be used to solve the equation set 130 using either of the forward DVI equations (FIG. 4) or the backward DVI equations (FIG. 5). In this approach, a lookup table of offsite velocity field at time t with variables that are out of the position on pixels (xi and yj) in an image is generated using the given matrix field {x.sub.ij(t.sub.1), y.sub.ij(t.sub.1)}. The lookup table in one example is an array of records, where each record contains four parameters: position coordinates x=i+x.sub.ij(t.sub.1) and y=j+y.sub.ij(t.sub.1), and the matrix values of x.sub.ij(t.sub.1) and y.sub.ij(t.sub.1) for all i[0,N.sub.x1]j[0,N.sub.y1]. All the values of the matrix field x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2) on the pixel points are then evaluated by a bilinear interpolation method and a searching algorithm. The searching algorithm finds three or four available values of the field x(x, y, t.sub.2) and y(x, y, t.sub.2) in the lookup table for all nearest points around the given point (i, j).

(80) Another technique involves solving the forward DVI equations (A1) by Newton-Raphson method, such as a damped Newton-Raphson method. Assuming that the forward DVI equations (A1) are nonlinear functions with variables x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2), the system of equations 130 can be solved by a Newton-Raphson method. Two-component nonlinear system of equations with variable x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2) are given by:

(81) $\{\begin{array}{c}{\mathrm{.Math.}}_1\left(x_{\mathrm{ij}}\left(t_2\right),y_{\mathrm{ij}}\left(t_2\right)\right)=x_{\mathrm{ij}}\left(t_2\right)+x\left(i+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0\\ {\mathrm{.Math.}}_2\left(x_{\mathrm{ij}}\left(t_2\right),y_{\mathrm{ij}}\left(t_2\right)\right)=y_{\mathrm{ij}}\left(t_2\right)+y\left(i+x_{\mathrm{ij}}\left(t_2\right),j+y_{\mathrm{ij}}\left(t_2\right),t_1\right)=0.\end{array}$

(82) Using the bilinear function in the above equation (A3) to expand the given field at time t.sub.1, all off site values (with non-integer value variables) of the given field at time t.sub.1 are evaluated by the function (10). Since both indexes p and q are functions of the variables x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2), the variables x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2) cannot be solved directly from the above equations. However, these equations are quasi-quadratic, and thus can be solved by a damped Newton-Raphson method. Iteration equations for solving the matrix field at time t for all i and j are given by:

(83) $\begin{array}{cc}\left(\begin{array}{c}{x}_{\mathrm{ij}}^{\left(m+1\right)}\left(t_2\right)\\ y_{\mathrm{ij}}^{\left(m+1\right)}\left(t_2\right)\end{array}\right)=\left(\begin{array}{c}{x}_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)\\ y_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)\end{array}\right)-{\left(\begin{array}{cc}\frac{{\mathrm{.Math.}}_1^{\left(m\right)}}{x_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}& \frac{{\mathrm{.Math.}}_1^{\left(m\right)}}{y_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}\\ \frac{{\mathrm{.Math.}}_2^{\left(m\right)}}{x_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}& \frac{{\mathrm{.Math.}}_2^{\left(m\right)}}{y_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}\end{array}\right)}^{-1}\left(\begin{array}{c}{\mathrm{.Math.}}_1^{\left(m\right)}\\ {\mathrm{.Math.}}_2^{\left(m\right)}\end{array}\right),& \left(\mathrm{A3}\right)\end{array}$
where m is an iteration index. All derivatives with respect variables x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2) in the above equations can be evaluated by the bilinear function (10). Two index variables p and q in function (10) are integer function of the variables x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2), but the derivative of the integer function is equal to zero, and thus:

(84) $\{\begin{array}{c}\frac{x_{\mathrm{pq}}\left(t_1\right)}{x_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}=\frac{x_{\mathrm{pq}}\left(t_1\right)}{y_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}0\\ \frac{y_{\mathrm{pq}}\left(t_1\right)}{x_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}=\frac{y_{\mathrm{pq}}\left(t_1\right)}{y_{\mathrm{ij}}^{\left(m\right)}\left(t_2\right)}0.\end{array}$

(85) A third approach can be used to solve equation set 130 employing the backward DVD equations (A2) (FIG. 5) by a New-Raphson method. In order to solve the velocity field x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2) in equation (A2), the equation can be rewritten as shown in the following equation (A4):

(86) $\begin{array}{cc}\{\begin{array}{c}x\left(x_1,y_1,t_1\right)+x\left(x_1+x\left(x_1,y_1,t_1\right),y_1+y\left(x_1,y_1,t_1\right),t_2\right)=0\\ y\left(x_1,y_1,t_1\right)+y\left(x_1+x\left(x_1,y_1,t_1\right),y_1+y\left(x_1,y_1,t_1\right),t_2\right)=0.\end{array}& \left(\mathrm{A4}\right)\end{array}$
where x.sub.1 and y.sub.1 are two components of a position vector. Two new index variables are introduced as follows:

(87) $\begin{array}{cc}\{\begin{array}{c}i=x_1+x\left(x_1,y_1,t_1\right)\\ j=y_1+y\left(x_1,y_1,t_1\right),\end{array}\mathrm{or}\{\begin{array}{c}{e}_1=x_1+x\left(x_1,y_1,t_1\right)-i=0\\ e_2=y_1+y\left(x_1,y_1,t_1\right)-j=0\end{array}.& \left(\mathrm{A5}\right)\end{array}$
The above equations become:

(88) $\begin{array}{cc}\{\begin{array}{c}{x}_{\mathrm{ij}}\left(t_2\right)=-x\left(x_1,y_1,t_1\right)\\ y_{\mathrm{ij}}\left(t_2\right)=-y\left(x_1,y_1,t_1\right).\end{array}& \left(\mathrm{A6}\right)\end{array}$

(89) All the displacement vector fields x.sub.ij(t.sub.2) and y.sub.ij(t.sub.2) in equation (A6) on all pixel points (i, j) can be determined by off-site field at time t.sub.1 after all position coordinates x.sub.1 and y.sub.1 are solved from equations (A5). According to the bilinear expansion in equation (A3), the displacement field at time t.sub.1 can be expressed by:

(90) $\left\{x\left(x,y,t_1\right),y\left(x,y,t_1\right)\right\}=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}\left\{x_{p+n_x,q+n_y}\left(t_1\right),y_{p+n_x,q+n_y}\left(t_1\right)\right\}{H}_{p+n_x,q+n_y}\left(x,y\right).$

(91) Iterative equations for solving position coordinate x.sub.1 and y.sub.1 are given by

(92) $\left(\begin{array}{c}{x}_1^{\left(m+1\right)}\\ y_1^{\left(m+1\right)}\end{array}\right)=\left(\begin{array}{c}{x}_1^{\left(m\right)}\\ y_1^{\left(m\right)}\end{array}\right)-{\left(\begin{array}{cc}\frac{e_1^{\left(m\right)}}{x_1^{\left(m\right)}}& \frac{e_1^{\left(m\right)}}{y_1^{\left(m\right)}}\\ \frac{e_2^{\left(m\right)}}{x_1^{\left(m\right)}}& \frac{e_2^{\left(m\right)}}{y_1^{\left(m\right)}}\end{array}\right)}^{-1}\left(\begin{array}{c}{e}_1^{\left(m\right)}\\ e_2^{\left(m\right)}\end{array}\right).$

(93) Using the property of the integer function that the derivative is zero yields:

(94) 0$\{\begin{array}{c}\frac{x_{\mathrm{pq}}}{x_1}=\frac{x_{\mathrm{pq}}}{y_1}0\\ \frac{y_{\mathrm{pq}}}{x_1}=\frac{y_{\mathrm{pq}}}{y_1}0.\end{array}$
All position coordinate x.sub.1 and y.sub.1 for given indexes i and j, and the displacement field at time t.sub.1 can be solved by only a few iteration steps, because these equations are quasi-quadratic for this motion model.

(95) In certain implementations, the Motion Vector Solver 110b is programmed to solve the forward displacement vector invariant equations (FIG. 4) in the equation set 130 by a damped Newton-Raphson method using the above iteration equations (A3) for solving the displacement field at time t.sub.2.

(96) Referring now to FIG. 11, an adaptive framework or process 700 is illustrated for solving a system of equations for motion estimation in accordance with further aspects of the present disclosure. At 702 in FIG. 11, a fundamental system of equations is provided for motion estimation. The process 700 can be used for motion estimation equation systems that are linear or nonlinear, and can be applied where the fundamental motion estimation equation system is under-constrained, over-constrained and/or fully constrained. Moreover, process 700 can be used for estimating any type or form of motion field, including without limitation displacement fields and/or velocity fields. For instance, the fundamental system can be provided at 702 as the DFD equations discussed above, or the fundamental equation system can be provided in the form of optical flow equations or heat flow equations and other nonlimiting embodiments. Modeling is then performed at 704, in which a motion field model is created based on the fundamental equation system using a bilinear interpolation function (e.g., 134 above). In this regard, the inventor has appreciated that for each pixel in an image, there may be one equation, or there may be two or more equations, and that motion vectors at individual pixel points are dependent upon motion vectors for neighboring node points. At 706 in FIG. 11, the block size is defined, where n is set to an integer value greater than or equal to 1 so that the number of motion vectors on node points is greater than the number of pixel points. In this manner, motion vectors on node points can be solved without solving all the motion fields for each pixel point. Thus the motion field for node points can be solved using the model (bilinear interpolation function). At 708 a cost function is created based on a linear or nonlinear least-squares principle. In addition, iteration equations are derived at 708 for solving the motion field on node points to convert a fully or over-constrained system to a fully constrained system based on a linear or nonlinear least-squares principle. In this regard, the fundamental equation system at 702 may be over-constrained, fully constrained and/or under-constrained, but the modeling at 704 causes the modeled motion field system to be fully constrained or over-constrained based on the block size definition at 706. Moreover, the inventor has appreciated that using the nonlinear least squares technique at 708, a fully constrained system of iteration equations can be created or derived for solving the motion field on node points, for example, using a Gauss-Newton method or Levenberg-Marquardt technique, where the iteration equations are based on the fundamental equation system provided at 702. At 712, the system of iteration equations is solved for the motion field on node points using an iterative technique with progressive relaxation of over-constraint (PROC), for instance, using a processor.

(97) In an alternative exemplary embodiment of the invention, a method for processing an image sequence can include an under-constrained system with a symmetric displaced frame difference equation (SDFD).

(98) As noted herein, there can be numerous local minima in image data applications when the nonlinear model is used to deal with a problem of global nonlinear minimization with huge number of unknown parameters. A convergent solution can depend on having a good guess for the initial values, especially for huge dimensional problem. However, in an exemplary embodiment of the invention, to avoid obtaining a specious local minimum solution when solving a nonlinear problem using the initial guess far away from the real flow field, a SDFD equation with only half displacement parameters used for resolving the total displacement field can be utilized.

(99) The estimated motion vector from two successive images is a displacement vector (or average velocity in the time interval) between the initial and final positions of a particle in the scene as shown in FIG. 12. FIG. 12 is a plot of a motion particle trajectory, instant velocities at time t.sub.1 and t.sub.2, and a displacement vector between initial and final positions at time t.sub.1 and t.sub.2. The trajectory represented by the curved, upper portion 1205 in FIG. 12 shows a very complex motion for a real particle for a longer temporal range. The motion fields discussed herein are displacement fields or their equivalent average velocity fields defined by v(r(t), t)=r(r(t), t)/t, where r and t=t.sub.2t.sub.1 are the displacement vector and a time difference between image frames at time t.sub.1 and t.sub.2. We denote the average velocity by v=v as the estimated velocity in this paper for convention.

(100) The motion field estimated from successive NOAA Advanced Very High Resolution Radiometer (AVHRR) images or ocean-color imagery spaced several hours apart, and the field is relevant to the important question of surface (or near-surface) drift velocity. For example, one application of this would be the drift of pollution (e.g., an oil plume).

(101) The image intensity I(r, t) measured by brightness (optical flow from visible band images) or temperature (heat flow from thermal images) can be a scalar function of position coordinates r(t) and time t. A brightness or tracer conservation constraint equation can be as follows:

(102) $\begin{array}{cc}\frac{I\left(r\left(t\right),t\right)}{t}=\left(\frac{}{t}+v.Math.\right)I\left(r\left(t\right),t\right)=0& \left(\mathrm{B1}\right)\end{array}$
where operator

(103) $\frac{}{t}$
denotes a total derivative with respect to time t. The differential form of equation (B1) can be linear in the velocity components and can hold only for infinitesimal or nearly infinitesimal motions. Using two successive frames, an estimate of partial derivatives with respect to time in (B1) can be made; however, equation (B1) may not exactly describe the motion between the two successive frames, especially for large-scale displacement motion.

(104) In order to constrain the image scenes at time t=t.sub.1 and t=t.sub.2, equation (B1) can be integrated from time t.sub.1 to t.sub.2

(105) $\begin{array}{cc}{}_{t_1}^{t_2}\frac{I\left(r\left(t\right),t\right)}{t}t=I\left(r\left(t_2\right),t_2\right)-I\left(r\left(t_1\right),t_1\right)0& \left(\mathrm{B2}\right)\end{array}$
where r(t.sub.1) and r(t.sub.2) can be the position vectors at time t.sub.1 and t.sub.2.

(106) In order to estimate a motion field with as small a displacement magnitude as possible, a half displacement vector r.sub.h=r/2=v t.sub.h can be defined along the total displacement vector 1210 as shown in FIG. 12, where t.sub.h=t/2, and v=v(r(t.sub.1+t.sub.h), t.sub.1+t.sub.h) can be the average velocity between t.sub.1 and t.sub.2 at position r(t.sub.1+t.sub.h). The initial and final position vectors can be expressed by

(107) $\{\begin{array}{c}r\left(t_1\right)=r\left(t_1+t_h\right)-v{t}_h\\ r\left(t_2\right)=r\left(t_1+t_h\right)+v{t}_h.\end{array}$
Substituting the position vectors in the above equations into equation (B2), the symmetric displaced frame difference (SDFD) equation can be given by
DFCD=I(r(t.sub.1+t.sub.h)+vt.sub.h,t.sub.2)I(r(t.sub.1+t.sub.h)vt.sub.h,t.sub.1)=0(B3)

(108) The conservation constraint equation (B1) can be integrated into path-independent equation (B3) and all terms in equation (B3) associate with intensity at time t.sub.1 and t.sub.2. It can be clear that employing the SDFD equation (B3) for motion estimations can achieve higher accuracy in comparison with the heat or optical flow equation (B1) especially for large scale displacement estimation. Since the estimated displacement vectors within the terms in equation (B3) are only half of the total displacement vector, convergence to a global optimal solution can be more robust than when the total displacement vector is estimated by an iteration algorithm.

(109) The intensity terms at time t.sub.1 and t.sub.2 in equation (B3) can be calculated from the intensity fields on the image sequence. The inverse problem for solving the velocity vector is under-constrained because two unknown velocity components must be derived from a single SDFD equation (B3) at each of these pixel points.

(110) In an exemplary embodiment of the invention, one of the efficient approaches for resolving the velocity field from the under-constrained system can be to express the velocity field as bilinear polynomial functions or two-dimensional B-Splines functions. Bilinear modeling of the displacement field can be utilized to develop a unified adaptive framework that can allow the estimator to control the resolutions of the displacement field and the degree of over-constraint, and to simplify the computational complexity. The velocity fields u(x, y, t) and v(x, y, t) can be represented as two-dimensional bilinear functions controlled by a smaller number of velocity estimates u(p, q, t) and v(p, q, t) which can lie on a coarser node point in grid as shown in FIG. 13. In FIG. 13, the node points (the small squares, see e.g., 1305, and block size n.sub.x=n.sub.y=4) of the bilinear polynomial functions with indices p and q are mapped on an image scene. All function values of off-node when ip(i) and jq(j) in the image plane can be calculated based on the sampled function values on node points when i=p(i) and j=q(j).

(111) In general, any two-dimensional function can be approximated by a Lagrange's bilinear function

(112) $\begin{array}{cc}f\left(x,y\right)=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}f\left(p+n_x,q+n_y\right)H_{p+n_x,q+n_y}\left(x,y\right)& \left(\mathrm{B4}\right)\end{array}$
where function H.sub.a,b(x, y) can be the bilinear forms. The image domain can be partitioned into sub-domains, or tiles, each of which can contain an n.sub.xn.sub.y array of pixels as represented in FIG. 13. The quantized indices p and q are functions of x and y and can be given by

(113) $\left\{p,q\right\}=\left\{p\left(x\right),q\left(y\right)\right\}=\left\{n_x.Math.\frac{x}{n_x}.Math.,n_y.Math.\frac{y}{n_y}.Math.\right\},$
where denotes an integer operator. The {p, q} serve as tile indices, since the integer operator increments them by unity after an additional n.sub.x or n.sub.y pixels are counted. The bilinear function H.sub.a,b (x, y) can be defined by

(114) $\begin{array}{cc}{H}_{a,b}\left(x,y\right)=\frac{1}{n_x{n}_y}\{\begin{array}{c}\left(n_x-x+p\right)\left(n_y-y+q\right)\left(a=p.Math.b=q\right)\\ \left(x-p\right)\left(n_y-y+q\right)\left(a=p+n_x.Math.b=q\right)\\ \left(n_x-x+p\right)\left(y-q\right)\left(a=p.Math.b=q+n_y\right)\\ \left(x-p\right)\left(y-q\right)\left(a=p+n_x.Math.b=q+n_y\right).\end{array}& \left(\mathrm{B5}\right)\end{array}$

(115) The two component velocity fields on pixels in an image can be approximated by the following discrete forms of the bilinear approximation functions with first order continuity that holds for all N=N.sub.xN.sub.y image globally

(116) $\begin{array}{cc}{v}_{\mathrm{ij}}=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}{v}_{p+n_x,q+n_y}{H}_{p+n_x,q+n_y}\left(i,j\right),& \left(\mathrm{B6}\right)\end{array}$
where a function with discrete variables i and j can be denoted by v.sub.ij=v.sub.ij(t)=v(i, j, t) or .sub.ij(t)=(i, j, t). In a special case, v.sub.ij can be equal to v.sub.pq for all indices i and j when block size parameter n is unity (n=n.sub.x=n.sub.y=1).

(117) All off-node velocities can be calculated by equation (B6) using the on-node velocities expressed as u.sub.pq and v.sub.pq. The SDFD equation in (B3) can become
DFCD.sub.ij=I(i+u.sub.ijt.sub.h,j+v.sub.ijt.sub.h,t.sub.2)+I(iu.sub.ijt.sub.h,jv.sub.ijt.sub.h,t.sub.1)=0,(B7)
where u.sub.ij and v.sub.ij are two components of velocities on each pixel with horizontal and vertical indices {i, j} for the whole image. The SDFD in (B7) is a function of variables of two velocity components that depend on node point velocities in (B6), i.e.
DFCD.sub.ij=DFCD.sub.ij(u.sub.ij,v.sub.ij)=DFCD.sub.ij(u.sub.pq,v.sub.pq).

(118) The SDFD equations (B3) and (B7) are equivalent only when the block size parameters n.sub.x=1 and n.sub.y=1, in this case, p=i and q=j. All independent SDFD.sub.ij equations have only a smaller number of the independent velocities on nodes when the velocity indices i=p(i) and j=q(j). The total number of SDFD equations is N=N.sub.xN.sub.y for an N.sub.xN.sub.y image sequence. The number of node points shown in FIG. 2 is given by

(119) $N_{\mathrm{node}}=\left(.Math.\frac{N_x-1}{n_x}.Math.+1\right)\left(.Math.\frac{N_y-1}{n_y}.Math.+1\right).$

(120) The total number of independent velocity fields with two components u.sub.pq and v.sub.pq is 2N.sub.node. It can be clear that this system is over-constrained because the number of SDFD equations in (B7) for all pixels is greater than the number of independent velocity components u.sub.pq and v.sub.pq (i.e., N>2N.sub.node) if block size n.sub.x>1 and n.sub.y>1.

(121) The presence of possible quantization errors and noise in the remote sensing measurements suggest that the SDFD equation in (B7) is never identically zero; however, we can choose a set of v.sub.ij for which it can be minimized in a least-squares sense. Accordingly, we define a Mean Square Error (MSE) function based on equation (B7) as

(122) 0$\mathrm{MSE}=\frac{1}{N}\underset{i,j}{.Math.}{\mathrm{DFCD}}_{\mathrm{ij}}^2,$
where i and j go over all pixels in N.sub.xN.sub.y image (i[0,N.sub.x1]j[0,N.sub.y1]). Minimizing the cost function with parameters u.sub.kl and v.sub.kl as variables for given indices k and l on all node points in an image represented in FIG. 13, a set of over-constrained system of equations for all cases n.sub.x>1 and n.sub.y>1 is given by

(123) $\begin{array}{cc}\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\mathrm{DFCD}}_{\mathrm{ij}}\left\{\frac{{\mathrm{DFCD}}_{\mathrm{ij}}}{u_{\mathrm{kl}}},\frac{{\mathrm{DFCD}}_{\mathrm{ij}}}{v_{\mathrm{kl}}}\right\}=0\mathrm{where}& \left(\mathrm{B8}\right)\\ \{\begin{array}{c}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}}{u_{\mathrm{kl}}}=t_h\frac{u_{\mathrm{ij}}}{u_{\mathrm{kl}}}\left[\begin{array}{c}\frac{I\left(x,j+v_{\mathrm{ij}}{t}_h,t_2\right)}{x}{|}_{x=i+u_{\mathrm{ij}}{t}_h}+\\ +\frac{I\left(x,j-v_{\mathrm{ij}}{t}_h,t_1\right)}{x}{|}_{x=i-u_{\mathrm{ij}}{t}_h}\end{array}\right]\\ \frac{{\mathrm{DFCD}}_{\mathrm{ij}}}{v_{\mathrm{kl}}}=t_h\frac{v_{\mathrm{ij}}}{v_{\mathrm{kl}}}\left[\begin{array}{c}\frac{I\left(i+u_{\mathrm{ij}}{t}_h,y,t_2\right)}{y}{|}_{y=j+v_{\mathrm{ij}}{t}_h}+\\ +\frac{I\left(i-u_{\mathrm{ij}}{t}_h,y,t_2\right)}{y}{|}_{y=j-v_{\mathrm{ij}}{t}_h}\end{array}\right]\end{array}& \left(\mathrm{B9}\right)\end{array}$
To obtain equation (B8), the following is used

(124) $\frac{u_{\mathrm{ij}}}{u_{\mathrm{kl}}}=\frac{v_{\mathrm{ij}}}{v_{\mathrm{kl}}}=\underset{=0}{\overset{1}{.Math.}}\underset{=0}{\overset{1}{.Math.}}{H}_{\mathrm{kl}}\left(i,j\right){}_{k,p+n_x}{}_{l,q+n_y},$
where symbol .sub.ij is the Kronecker-Delta symbol, and the summation domain is reduced from the entire image plane to only a local region .sub.kl, so that

(125) $\underset{i,j}{.Math.}->\underset{i,j{}_{\mathrm{kl}}}{.Math.}=\underset{i=k-n_x+1}{\overset{k+n_x-1}{.Math.}}\underset{j=l-n_y+1}{\overset{l+n_y-1}{.Math.}}.$
The two independent equations in (B8) on a node point can be degraded back to a single SDFD equation (B7) when block size n.sub.x and n.sub.y are equal to unity. Of the original pixel set, the velocities u.sub.kl and v.sub.kl and their nearest neighbors are the only members contained in the summations in equation (B8).

(126) To solve the nonlinear system equations (B8), the SDFD can be expanded in Taylor series using Gauss-Newton method

(127) ${\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m+1\right)}{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}+\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{u_{\mathrm{kl}}^{\left(m\right)}}\left(u_{\mathrm{kl}}^{\left(m+1\right)}-u_{\mathrm{kl}}^{\left(m\right)}\right)+\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{v_{\mathrm{kl}}^{\left(m\right)}}\left(v_{\mathrm{kl}}^{\left(m+1\right)}-v_{\mathrm{kl}}^{\left(m\right)}\right)$
where m is an iteration index, and
DFCD.sub.ij.sup.(m)=DFCD.sub.ij(u.sub.pq.sup.(m),v.sub.pq.sup.(m))
Substituting SDFD expansion into equation (B8), the following iterative equations of two component velocity u.sub.kl and v.sub.kl can be found for all indexes k and l on node points

(128) $\begin{array}{cc}{v}_{\mathrm{kl}}^{\left(m+1\right)}=v_{\mathrm{kl}}^{\left(m\right)}-{\left(A_{\mathrm{kl}}^{\left(m\right)}\right)}^{-1}{B}_{\mathrm{kl}}^{\left(m\right)}\mathrm{where}{A}_{\mathrm{kl}}^{\left(m\right)}=\left(\begin{array}{cc}\left(+1\right)\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\left(\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{u_{\mathrm{kl}}^{\left(m\right)}}\right)}^2& \underset{i,j{}_{\mathrm{kl}}}{.Math.}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{u_{\mathrm{kl}}^{\left(m\right)}}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{v_{\mathrm{kl}}^{\left(m\right)}}\\ \underset{i,j{}_{\mathrm{kl}}}{.Math.}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{u_{\mathrm{kl}}^{\left(m\right)}}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{v_{\mathrm{kl}}^{\left(m\right)}}& \left(+1\right)\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\left(\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{u_{\mathrm{kl}}^{\left(m\right)}}\right)}^2\end{array}\right)\mathrm{and}{B}_{\mathrm{kl}}^{\left(m\right)}=\left(\begin{array}{c}\underset{i,j{}_{\mathrm{kl}}}{.Math.}{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{u_{\mathrm{kl}}^{\left(m\right)}}\\ \underset{i,j{}_{\mathrm{kl}}}{.Math.}{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}\frac{{\mathrm{DFCD}}_{\mathrm{ij}}^{\left(m\right)}}{v_{\mathrm{kl}}^{\left(m\right)}}\end{array}\right)& \left(10\right)\end{array}$
where 0 is a Levenberg-Marquardt factor that can be adjusted at each iteration to guarantee that the objective function MSE is convergent. A smaller value of the factor can be used, bringing the algorithm closer to the Gauss-Newton method with second order converging. This Levenberg-Marquardt method can improve converge properties efficiently in practice and has become the standard of nonlinear least-squares routines.