Methods and apparatus for motion estimation in motion picture processing
The U.S. government has rights in this invention pursuant to NSF Grant No. DC1-8407285. 
    
    
     BACKGROUND OF THE INVENTION 
     The technical field of this invention is signal processing and, in particular, the estimation of motion in motion picture frames. 
     A motion picture is composed of a sequence of still frames which are displayed in rapid succession. The frame rate necessary to achieve proper motion rendition in typical visual scenes is sufficiently high that there is a great deal of temporal redundancy among adjacent frames. Most of the variation from one frame to the next is due to object motion. This motion may occur within the scene or relative to the camera which generates the sequence of still frames. 
     There are a wide variety of applications where one desires to manipulate a motion picture by exploiting the temporal redundancy. In order to do this, it is necessary to account for the presence of motion. A number of systems are known for determining the movement of objects within the sequence of still frames. The process of determining the movement of objects within image sequences is known as motion estimation. 
     In general, prior systems for motion estimation have encountered three primary problems. First, such systems often have difficulty in accurately estimating motion in noisy images. Many systems are explicitly formulated under the assumptions of high signal-to-noise level. As a consequence, if the systems are applied to noisy pictures, the motion estimation errors are typically large. Most motion-compensated, image-processing systems require very accurate motion estimates in order to maintain adequate picture quality. 
     Second, objects within the image sequence which move with large velocities are difficult to track. In real-life motion pictures, the velocity field is a complicated function of spatio-temporal position; different parts of a motion picture move with different velocities. Therefore, most motion estimation systems are based on local operations (i.e., motion is estimated over a small region of the image). One of the problems with this approach is that typically only small velocity fields can be estimated reliably. 
     Third, the need to reduce computational complexity is an ever-present problem. Many applications of motion compensation require real-time operation. For real-time operation to be feasible, it is necessary for the system to be computationally efficient. Even in those applications where real-time operation is not required, computational complexity is an important characteristic which affects the cost of implementing a specific motion estimation system. 
     There exists a need for better motion estimation systems. Systems that address the problems of background noise, large velocities and computational complexity would satisfy a long-felt need in the art. It is an object of the present invention to provide better motion estimation systems for image processing purposes, including speed alteration, still image enhancement, picture coding and various other applications. 
     In the description of the present invention, matrix and vector notations are used. Matrices are represented with upper case symbols (A, B, etc.) and vectors are represented with either upper or lower case symbols with a bar over the symbol (a, b, S etc.). For example, a set of linear equations is written as 
     
         Ax=b. 
    
     The inverse of a matrix A is written as A -1 , and the transpose is written as A T . Entries to a matrix are referred to with indexed notation. Therefore, A ij  or A(i,j) refers to the i th  row and j th  column of matrix A. 
     All vectors are column vectors. Entries of a vector are referred to with subscripted notation. Therefore, b i  refers to the i th  element of vector b. 
     The signals are either single images or sequences of images which comprise a motion picture. The luminance of an image is a function of two variables, x and y. For the sake of notational convenience, the tuple (x,y) will be written as x in many occasions. Therefore, the image s(x,y) is equivalent to the image s(x). Continuous sequences of images are written as s(x,y,t)=s(x,t). Therefore, s(x,t o ) refers to the frame at the time instant t o . 
     SUMMARY OF THE INVENTION 
     A motion picture defined by a plurality of picture frames can be processed in a variety of ways by first determining the movement of objects within the image sequence. Motion estimation techniques for determining the velocity field of predetermined regions of the image are disclosed herein. These estimation techniques are applicable to noise reduction, frame rate conversion, and picture coding, 
     Let I(x,y,t) denote the intensity of a motion picture frame at time t and spatial position (x,y). In the present invention, it is assumed that a motion picture frame I(x,y,t o ) can be mapped into I(x,y,t) with the mapping functions a x  (x,y,t) and a y  (x,y,t) as follows: ##EQU1## The mapping functions a x  and a y  determine corresponding points between frame I(x,y,t o ) and frame I(x,y,t). In effect, the motion estimation problem is to determine these mapping functions. An important consequence of the relationship indicated in Equation (1) above is that the velocity field is related to I(x,y,t) as follows: ##EQU2## 
     In many cases of practical interest, the velocity field is sufficient to relate corresponding points between two adjacent frames. 
     The motion estimation methods and apparatus of the present invention are directed to the determination of the components of velocity field, v x  and v y  at arbitrary points between two given frames. 
     By assuming that the velocity components v x  and v y  are constant over a region Ψ, equation (2) above, is reduced to: ##EQU3## The function set forth on the left hand side of equation 2A is defined herein as the directional gradient residual (&#34;DGR&#34;). In the present invention, a least squares estimate of v x  and v y  is obtained by minimizing the average value of (DGR) 2  over the region Ψ: ##EQU4## 
     Equation (3) forms the basis for the motion estimation method disclosed herein. There are five novel features of this method which are discussed in the following sections: 
     (1) Signal models are used to solve Equation (3) given only image intensity samples (as opposed to continuous functions). 
     (2) Numerical stability is achieved by incorporating an eigenvalue/eigenvector decomposition in the linear equations implied by Equation (3). 
     (3) Motion estimation accuracy is improved by using weighted averaging of individual velocity estimates. 
     (4) Region incompatibility is detected by examination of the least squares residual. 
     (5) Large velocities are estimated by applying the method using multigrid techniques. 
     Because the least squares estimation approach involves quadratic functions of the parameter values, the optimal velocity components v x  and v y  of the velocity field v can be determined from a set of linear equations: 
     
         Wv=γ                                                 (4) 
    
     where W is a covariance matrix and γ is a correlation vector. W and γ are given by the following equations: ##EQU5## 
     As will be described in more detail below, W and γ can be derived solely from the image intensity or luminance by least squares signal estimation techniques. The solution technique applied to equation (4) is referred to hereinafter as &#34;eigenvalue/eigenvector decomposition.&#34; 
     In another aspect of the invention, signal estimation techniques are disclosed for obtaining a set of signal coefficients S i  from samples of I(x,y,t). A three dimensional parametric signal model is employed which assumes that the signal in region Ψ can be represented by the following equation: ##EQU6## where S i  are a set of signal coefficients and Φ i  (x,y,t) are predetermined basis functions. A least squares estimate of S i  is obtained by minimizing the error function: ##EQU7## A linear relationship between the model parameters and the signal exists, which greatly simplifies the calculations. With this approach, the available samples are used to estimate the signal coefficients. These signal coefficients are used to compute the entries for the covariance matrix W and the correlation vector γ, which are then used to estimate the velocity field. 
     In yet another aspect of the invention, techniques are described below for addressing the problem of uniform luminance (background) values throughout sampled region and the special problem of motion estimation in the presence of edges. In both of these instances the eigenvectors and eigenvalues of the covariance matrix W can be used to resolve the ambiguity. In particular, for the case of edges the eigenvectors of W point in the directions of minimum and maximum contrast. By applying a singular value decomposition (SVD) to the matrix W, an estimate of the velocity field in the direction orthogonal to the edge can be generated. 
     The invention can also be employed in a multigrid fashion to progressively refine the velocity field estimates. This technique is particularly useful in determining large velocities, with both high resolution and high accuracy. In the multigrid approach, the velocity field is first determined over the entire image on a coarse grid. The coarse grid is obtained from the original frames by down-sampling the images. Down-sampling the images has the property of contracting the velocity field; large velocities in the original frames become small velocities in the down-sampled frames. 
     The coarse velocity field is interpolated to generate initial estimates of the velocity field at a finer grid. The velocity field is then calculated by the methods of the present invention again using a finer subset of samples (i.e. a smaller down-sampling factor). The process can be repeated at successively finer grids until the sampling factor is unity. 
     Also disclosed herein are techniques for improving velocity estimates by weighted averaging; techniques for selecting the region Ψ on each of the comparison frames from an initial velocity estimate; and techniques for detecting regions of inaccurate motion estimates. 
     The invention will be described next in connection with certain preferred embodiments. In particular, the overall apparatus and methods will be described and applied to the specific problems of image synthesizing (i.e. between frame images), noise reduction and picture coding. However, it should be clear that various changes, additions and subtractions can be made without departing from the spirit or scope of the invention. For example, the method can be applied to color images by letting the image intensity values be equal to the luminance component. Moreover, although a set of three dimensional polynomials are used herein to provide a continuous signal representation of the image samples, other representations can also be employed. The continuous representation need only be a linear combination of predefined basis functions which are differentiable and linearly independent. 
     Background materials and further details on particular implementations can be found in Appendix A, a doctorate thesis by D. Martinez, a co-inventor herein. Appendix A is submitted herewith and incorporated by reference. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     FIG. 1 is a schematic illustration of a motion estimation apparatus according to the invention. 
     FIG. 2 is a more detailed schematic illustration of the velocity estimating component of the motion estimation apparatus of FIG. 1. 
     FIG. 3 is a flow diagram of a multigrid velocity estimation technique according to the present invention. 
     FIG. 4 is a schematic illustration of an image synthesizing apparatus according to the invention. FIG. 4A is a schematic illustration of the operation performed by the system of FIG. 4. 
     FIG. 5 is a schematic illustration of an image noise reduction apparatus according to the invention. FIG. 5A is a schematic illustration of the operation performed by the system of FIG. 5. 
     FIG. 6 is a schematic illustration of an image coding apparatus according to the invention. FIG. 6A is a schematic illustration of the operation performed by the system of FIG. 6. 
    
    
     DETAILED DESCRIPTION 
     FIG. 1 is a block diagram of a motion estimation apparatus 10 according to the present invention including an analog-to-digital (A/D) converter 12, a signal sampler 14, a signal estimator 16, a velocity estimator 18, and a velocity interpolator 20. The A/D converter 12 serves to digitize the image in any given frame into individual pixels (e.g. a 512×512 grid of image pixels) and assigns to each pixel a luminance value (e.g. ranging from 0 &#34;black&#34; to 255 &#34;white&#34;). The signal sampler 14 down-samples the pixels by a predefined factor. For example, by down-sampling by a factor of 2, one out of four pixels is sampled; similarly, down-sampling by a factor of four results in one out of sixteen pixels being sampled. 
     The signal estimator 16 takes the sampled digital values and creates continuous functions, which are used in turn by the velocity estimator 18 to obtain the velocity estimates over particular regions (e.g., 5×5 sample regions) of the image from one frame to the next. The velocity estimates are then used by interpolator 20 to obtain velocity estimates for the unsampled pixels. The process can be repeated with successively smaller down-sampling factors to obtain refined motion estimates. 
     FIG. 2 is a more detailed diagram of an implementation of the motion estimating apparatus of the FIG. 1. As illustrated in FIG. 2, first and second frame memories, 30 and 32 respectively, provide digitized stores of image pixels, which are sampled by first and second samplers, 34 and 36 respectively, to obtain two sets of down-sampled digitized image data, I d  (x,y,t o ) and I d  (x,y,t 1 ). 
     At this point it is necessary to select the samples from frame I(x,y,t o ) and frame I(x,y,t 1 ) which will define the region Ψ to be analyzed. For each point in the image, an initial velocity vector stored in the initial velocity memory 38 is employed to project a region Ψ onto the first and second frames. The velocity interpolator 40 insures that the projectory is appropriately scaled in accord with the down-sampling factor. The end result of this process is a set of samples from the first and second frames which are initially presumed to represent the movement of region Ψ of the image from time t o  to time t 1 . This set is organized in vector format by sample selector 42 to establish a sample vector, I, and subjected to matrix/vector multiplication by multiplier 44 to obtain signal coefficients, S i . Quadratic evaluator 46 derives a covariance matrix W and a correlation vector γ from the signal coefficient vector S, and from these parameters, velocity computer 48 computes a new velocity estimate, v, which is stored in store 22 and used to update the initial velocity memory 38. 
     The sample selector 42 defines the region Ψ over which motion is to be estimated. At each point (x,y,t) an intial velocity estimate v*(x,y,t) is obtained from velocity memory 38 and interpolator 40. (For system initialization, the velocity typically is assumed to be zero). The selector 42 selects frames I(x,y,t o ) and I(x,y,t 1 ) from frame memory elements 30 and 32. The frames are selected such that t o  is less than t 1 , and ##EQU8## From these two displacement vectors, two sets of LxL image samples (typically 5×5 matrices) are obtained from the frames, with each window centered about the displacement vector so that: ##EQU9## where L=2M+1. 
     Multiplier 44 performs a matrix/vector multiplication using the samples defined by selector 42 to obtain the signal coefficients S i . As noted above the system assumes that the signal can be represented by the equation (7) set forth above, and that a least squares estimate of S i  can be obtained by minimizing the following error function set forth in equation (8) above. 
     In multiplier 44 the minimizing operation is performed by matrix multiplication such that: ##EQU10## where Q st  =(A st   T  A st ) -1  A st  and ##EQU11## 
     The basis functions Φ i  (x,y,t) are selected to perform signal smoothing and should be easily differentiable. The rank of A st  is equal to P, the number of basis functions. One exemplary choice of basis functions is the set of three dimensional polynomials: ##EQU12## 
     Once multiplier 44 has generated the signal coefficients S i , evaluator 46 can employ this continuous function to obtain W and γ. Thus, W and γ can be computed as follows: ##EQU13## 
     Once W and γ are determined, the velocity components v x  and v y  of the velocity field are defined by the linear equation Wv=γ. 
     The velocity estimates can be further refined by weighted averaging using the eigenvalues of the covariance matrix W. This averaging technique is based upon the observation that high contrast images yield better velocity estimates and the eigenvalues of W are proportional to the contrast. Let λ 1  and λ 2  be the eigenvalues of W at some point, such that λ 1  ≧λ 2 . An averaged velocity field can be derived for each point according to the following equation: ##EQU14## where  and λ 1  (i) and λ 2  (i) are the eigenvalues of W at the ith position and the summation from J=1 to J=N represents the eigenvalues for the jth position in vicinity of the desired point. 
     The eigenvalues and eigenvectors of the covariance matrix W can also be used to address particular problems. If the eigenvalues λ 1  and λ 2  and eigenvectors α 1  and α 2  of the symmetric 2×2 matrix W are computed, then a uniform background is observed whenever neither λ 1  nor λ 2  exceed a threshold value. In such cases a zero velocity field can be assumed. Moreover, an edge is detected when λ 1  greatly exceeds λ 2 . 
     In the case of edge detection, a more effective equation for estimating the velocity vector can be employed as follows: ##EQU15## In all other cases the normal velocity vector computation, i.e. v=W -1 .sub.γ, can be employed. 
     In FIG. 3 a flow diagram of a multigrid velocity estimation technique is shown for refinement of velocity estimates. This technique is particularly useful in determining large velocities, with both high resolution and high accuracy. In the multigrid approach, the velocity field is first determined over the entire image on a coarse grid. The coarse grid is obtained from the original frames by down-sampling the images. Down-sampling the images has the property of contracting the velocity field; large velocities in the original frames become small velocities in the down-sampled frames. For example, as shown in FIG. 3, the images can first be down-sampled by a factor of four (i.e. every sixteenth pixel). The velocity estimate computed in this fashion replaces the initial velocity estimate. 
     The coarse velocity field is interpolated to generate an estimate of the velocity field at a finer grid. The velocity field is then calculated by the methods of the present invention again using a finer subset of samples (i.e. as shown in FIG. 3, down-sampling by a factor of two). The process can be repeated at a successively finer grids until the sampling factor is unity and a final velocity field is obtained. 
     In FIG. 4 an apparatus for generating synthetic frames is shown. This embodiment is particularly useful for time scale modification of motion pictures. By creating synthetic frames and adding or substituting such frames for all or a portion of the original frames, one can lengthen or shorten the time taken to play or broadcast a motion picture. The synthetic frame generating system 50 of FIG. 4 includes a first frame memory 52 and a second frame memory 54, a motion estimator 56 and frame interpolator 58. Motion estimator 56, similar to the system 10 shown in FIGS. 1 and 2, extracts samples from the frame memories 52 and 54 and generates a velocity field. Frame interpolator 58 then produces a synthetic frame output from the velocity field and the frame memories. 
     The synthetic frame generating system 50 computes an intermedate velocity field v(x o  y o ,t) between original frames I(x,y,t o ) and I(x,y,t 1 ) in accordance with the methods disclosed above. Displacement vectors, d o  and d 1 , are again defined as described above by the equations: ##EQU16## From the displacement vectors d o  and d 1 , frame interpolator 50 computes the points I(x o  -d o ,t o ) and I(x 1  -d 1 ,t 1 ). To obtain intermediate frame values, the following formula can be employed: ##EQU17## FIG. 4A is a schematic illustration of the operation performed by system 50 of FIG. 4. 
     In FIG. 5 a noise reduction system 60 is shown, particularly useful in preparing still &#34;snapshots&#34; from individual motion picture frames. The noise reduction system 60 of FIG. 5 includes first, second and third frame memories 62 ,64, 66, respectively, motion estimators 68 and 70, frame interpolators 72 and 74, and filter 76. 
     Motion estimator 68 estimates a velocity vector v o  from samples extracted in the first and second frame memories 62, 64. Similarly, motion estimator 70 estimates a velocity vector v 1  from sample stored in the second and third frame memories 64, 66. 
     In frame interpolator 72, velocity vector v o  is used to compute a synthetic version of the middle (second) frame. In a similar fashion, frame interpolator 74 generates another synthetic version of the second frame based in this instance on velocity vector v 1 . As before, displacement vectors are computed from the velocity vectors, viz: ##EQU18## 
     The two synthetic versions and the original frame data from the second frame memory are then processed by filter 76 (e.g. by averaging or median selection) to obtain the final, noise-reduced, output. This technique can be used to obtain enhanced snapshots of any particular frame, or can be employed to provide continuous, real-time processing prior to display, transmission, etc. FIG. 5A is a schematic illustration of the operation preformed by the system 60 of FIG. 5. 
     In FIG. 6, a motion picture coding system 80 employing the present invention is shown. Transmitter 82 includes a two frame memory 84, frame predictor 86, quantizer 88, and summers 90 and 92. Receiver 94 includes another two frame memory 96, frame predictor 98 and summer 100. In the transmitter 82, frame predictor 86 estimates a velocity v from the two frame memory 84 and uses it to synthesize a predicted frame. The predicted frame is subtracted from the next actual frame in summer 90 and the differences only are processed by quantizer 88 for transmission. (The actual frame can be reconstructed by adding the predicted frame and the quantizer output using summer 92 in order to update frame memory 84.) 
     In a similar fashion, frame predictor 98 in receiver 94 estimates velocity v from two frame memory 96 and likewise uses it to synthesize a predicted frame. The predicted frame is added to the transmitted difference signals in summer 100 to obtain a picture output and update the receiver&#39;s frame memory 96. 
     Both frame predictor 86 and frame predictor 98 predict a frame at time t from frames at t o  and t 1 , where t is greater than t 1 , and t 1  is greater than t o . As before, a velocity vector is computed from the image data and a displacement vector d is calculated as follows 
     
         d=v(x.sub.0, y.sub.0, t)·t.sub.1 -t)              (9C) 
    
     The predicted value in frame t is 
     
         I(x,t)=I(x.sub.0 -d,t.sub.1)                               (18A) 
    
     FIG. 6A is schematic illustration of the frame prediction operations performed by predictors 86 and 98 in FIG. 6. 
     In each of the above-described implementations, motion estimation errors can be detected by examination of the least squares residual function: ##EQU19## 
     If the residual exceeds a predetermined threshold, the estimate should not be used. For frame interpolation applications, a higher-than-threshold residual value, can trigger the system to project zero velocity. In noise reduction applications, a high residual value can trigger suppression of the frame filter. For motion picture coding, the high residual value can suppress frame prediction.