DCT-based Motion Estimation Method

This application makes reference to the following U.S. application Ser. No. 08/937,100, entitled "Optimal Unified Architectures for the Real-Time Computation of Time-Recusive Discrete Sinusoidal Transforms," filed Aug. 31, 1992, and issued as U.S. Pat. No. 5,339,265 on Aug. 16, 1994. This patent is hereby incorporated by reference.

1. Field of the Invention

The present invention relates generally to a motion estimation system and method utilizing sinusoidal orthogonal principals to estimate displacements of moving objects in the transform domain based upon the concept of pseudo phases, and more particularly to an improved discrete cosine transform (DCT) based motion estimation system and method.

2. Description of the Prior Art

In recent years, great interests have been found in motion estimation due to its various promising applications in high-definition television (HDTV), multimedia, video telephony, target detection and tracking, computer vision, etc. Extensive research has been conducted over many years in developing new methods and designing cost-effective and massively parallel hardware architectures suitable for current VLSI technology.

The most commonly used motion estimation method in video coding is the Full Search Block Matching Method (BMA-ME) which searches for the best block among all the blocks in a search area of larger size in terms of either the mean-square error or the mean of the absolute frame difference. The computational complexity of this approach is very high, i.e., O(N.sup.4) for a N.times.N block. Even so, BMA-ME has been successfully implemented on VLSI chips. To reduce the number of computations, a number of suboptimal fast block matching methods have been proposed. However, these methods require three or more sequential steps to find suboptimal estimates. Recently, a correlation-based approach (CLT-ME) using Complex Lapped Transform (CLT) to avoid the block effect was proposed, but it still requires searching over a larger search area and, thus, results in a very high computational burden. Moreover, motion estimation using the CLT-ME is accurate on moving sharp edges but not on blurred edges.

In addition to block-based approaches, pel-based estimation methods such as the Pel Recursive Method (PRA-ME) and the Optical Flow Approach (OPA-ME), are very vulnerable to noise by virtue of their involving only local approaches and may suffer from the instability problem. For multiframe motion detecting, 3D-FFT has been successfully used to estimate motion in several consecutive frames based on the phenomenon that the spatial and temporal frequencies of a moving object lie on a plane of spatiotemporal space. However, this requires processing of several frames rather than two, and the Fast Fourier Transform (FFT) operates on complex numbers and is not used in most video standards.

It is therefore an object of the present invention to provide a DCT motion based estimation system which reduces the computational complexity from O(N.sup.4) to O(N.sup.2) for a N.times.N block.

It is a further object to provide a highly parallel pipelined architecture implementation of a DCT video coding-based system.

It is yet another object to provide improve system throughput by eliminating the need for an inverse-discrete cosine transform (IDCT).

It is yet another object to remove the discrete cosine transform from the feedback loop in order to provide a digital video system for high-end real-time applications.

According to one broad aspect of the present invention, there is provided a method for estimating the translational motion of an object, the method comrising the steps of:

computing DCT-I and DST-I coefficients of x.sub.1 (n) and the DCT-II and DST-II coefficients of x.sub.2 (n);

computing the pseudo phase g.sub.m.sup.S (k) for k=, . . . , N by solving this equation: ##EQU1## feeding the computed pseud phase, {g.sub.m.sup.s (k); k=1, . . . , N}, into an IDST-II decoder to produce an output {d(n); n=0, . . . , N-1}, searching for the peak value; and

estimating the displacement m by substituting the peak position into: ##EQU2## where i.sub.p =arg max.sub.n .linevert split.d(n).linevert split. is the index at which the peak value is located.

According to another broad aspect of the invention, there is provided a method for estimating the translational motion of an object, the method comrising the steps of:

providing previous frame x.sub.t-1 and the current frame x.sub.t information into 2D-DCT-II and 2D-DCT-1 coders, respectively;

computing four coefficient, DCCTII, DCSTII, DCSTII and DSSTII, each of which is defined as a two-dimensional separable function formed by 1D-DCT/DST-II kernels: ##EQU3## or symbolically,

computing four coefficients, the two-dimensional DCT coefficients of the first kind (2D-DCT-I) are calculated based on 1D-DCT/DST-I kernels: ##EQU4## or symbolically,

determining the respective pseudo phase functions from the following equations: ##EQU5## passing these two pseudo phase functions through 2D-IDCT-II coders (IDCSTII and IDSCTII) to generate two functions, DCS(.cndot.,.cndot.) and DSC(.cndot.,.cndot.): ##EQU6## searching for the peak values; and estimating the displacement by substituting the peak positions into: ##EQU7##

Other objects and features of the present invention will be apparent from the following detailed description of the preferred embodiment.

The invention will be further described in conjunction with the accompanying drawings, in which:

FIG. 1a is a prior art MPEG-compliant coder;

FIG. 1b is a motion-compensated DCT hybrid coder with motion estimation constructed in accordance with a preferred embodiment of the invention;

FIGS. 2a and 2b illustrate a right and left shift motion of an object, without noise, respectively;

FIG. 3a is a logic flow diagram of the method associated with a 1-dimensional DCT motion estimation system constructed in accordance with a preferred embodiment of the invention;

FIGS. 3a to 3f illustrate a right shift motion of a random signal, with noise;

FIGS. 3g to 3l illustrate a left shift motion of a random signal with noise;

FIG. 4 is an illustration of an object which moves translationally by m.sub.1 in X direction and n.sub.1 in Y direction as viewed on a camera plane and within the scope of a camera in a noiseless environment;

FIG. 5a is a logic flow diagram of the method associated with a 2 dimensional DCT motion estimation system constructed in accordance with a preferred embodiment of the invention;

FIG. 5b is a two-dimensional motion-compensated DCT hybrid coder with motion estimation constructed in accordance with a preferred embodiment of the invention;

FIGS. 6a to 6d are graphs of DCS functions which represent the associated pseudo phase functions;

FIGS. 6e to 6h are graphs of DSC functions representing the associated pseudo phase functions;

FIGS. 7a and 7b are images of a rectangularly-shaped moving object with arbitrary texture;

FIGS. 7c and 7d are illustrations of FIGS. 7a and 7b, respectively, which are corrupted by additive white Gaussian noise at SNR=10 dB;

FIGS. 7e and 7f are graphs of the resulted pseudo phase functions f and g;

FIGS. 7g and 7h are graphs of the DCS and DSC functions;

FIGS. 8a and 8b are graphs of the resulted pseudo phase functions, f and g, for motion in the direction (5, -3) in a uniformly bright background (c=3).

FIGS. 8c and 8d are graphs of the DCS and DCS functions for motion in the direction (5, -3) in a uniformly bright background (c=3);

FIGS. 8e and 8f are the graphs of DCS and DSC for the motion (8, 7) with background pixel level, c=255.

FIG. 9a is a block diagram of an Extended DXT-ME Estimator (E-DXT-ME);

FIG. 9b is a block diagram of a Simplified Extended DXT-ME Estimator (SE-DXT-ME);

FIGS. 10a, 10b, and 10c are graphs of the reference block and search area, the search range and DCS/DSC for the overlapping approach, respectively;

FIGS. 11a to 11p are comparisons of DXT-ME of size 64 by 64 with a Full Search Block matching Method (BMA-ME) of block size (bs=16 pels) but different search areas on a noisy small car with the noise being either 10 dB or 0 dB, respectively;

FIGS. 12a, 12b, and 12c are graphs of the original, the edge extracted and the frame differentiated images of Frame 57 in the sequence "Flower Garden";

FIGS. 13a and 13b are comparisons of SE-DXT-ME with BMA-ME on the sequence "Flower Garden" for different block sizes and preprocessing functions;

FIGS. 14a to 14d are comparisons of SE-DXT-ME and Overlapping SE-DXT-ME with BMA-ME, respectibely, on the small "Flower Garden" for different block sizes and preprocessing functions;

FIGS. 15a, 15b, and 15c are comparisons of the original, the edge-extracted and the frame-differentiated imaegs in the sequence "Infrared Car";

FIGS. 16a to 16d are comparisons of SE-DXT-ME and Overlapping SE-DXT-ME with BMA-ME, respectively, on the sequence "Infrared Car" for different block sizes and preprocessing functions;

FIGS. 17a to 17d are comparisons of SE-DXT-ME and Overlapping SE-DXT-ME with BMA-ME, respectively, of the sequence "Miss America" for different block sizes and preprocessing functions; and

FIGS. 18a and 18b are comparisons of SE-DXT-ME and Overlapping SE-DXT-ME with BMA-ME, respectively, on the sequence "Table Tennis" for different block sizes and preprocessing functions.

The invention will now be described in conjunction with the accompanying drawings wherein like numerals represent like parts.

In most international video coding standards such as CCITT H.261, MPEG, as well as the proposed HDTV standard, Discrete Cosine Transform (DCT) and block-based motion estimation are the essential elements to achieve spatial and temporal compression, respectively. Most implementations of a standard-compliant coder adopt the structure as shown in FIG. 1a. The DCT is located inside the loop of temporal prediction, which also includes an Inverse DCT (IDCT) and a spatial-domain motion estimator (SD-ME) which is usually the BMA-ME. The IDCT is needed solely for transforming the DCT coefficients back to the spatial domain in which the SD-ME estimates motion vectors and performs motion compensated prediction. This is an undesirable coder architecture for the following reasons. In addition to the additional complexity added to the overall architecture, the DCT and IDCT must be put inside the feedback loop which has long been recognized as the major bottleneck of the entire digital video system for high-end real-time applications. The throughput of the coder is limited by the processing speed of the feedback loop, which is roughly the total time for the data stream to go through each component in the loop. Therefore the DCT (or IDCT) must be designed to operate at least twice as fast as the incoming data stream. A compromise is to remove the loop and perform open-loop motion estimation based upon original images instead of reconstructed images in sacrifice of the performance of the coder.

An alternative solution without degradation of the performance is to develop a motion estimation method which may work in the DCT transform domain. In this way, the DCT may be moved out of the loop as depicted in FIG. 1b and thus the operating speed of this DCT may be reduced to the data rate of the incoming stream. Moreover, the IDCT is removed from the feedback loop which now has only two simple components Q and Q.sup.-1 (the quantizers) in addition to the transform-domain motion estimator (TD-ME). This not only reduces the complexity of the coder but also resolves the bottleneck problem without any trade-off of the performance. Furthermore, different components may be jointly optimized if they operate in the same transform domain.

A novel method for motion estimation, called the DCT-Based Motion Estimation (DXT-ME) to estimate motion in the discrete cosine transform or discrete sine transform (DCT/DST or DXT, for short) domain is disclosed herein. DXT-ME is based on the principle of orthogonality of sinusoidal functions. This new method has certain merits over conventional methods. It has very low computational complexity (on the order of N2 compared to N4 for BMA-ME) and is robust even in a noisy environment. This method takes DCT coefficients of images as input to estimate motions and therefore may be incorporated efficiently with the DCT-based coders used for most current video compression standards. As explained before, this combination of both the DCT and motion estimation into a single component reduces the coder complexity and, at the same time, increases the system throughput. Finally, due to the fact that the computation of pseudo phases is inherently highly local operation, a highly parallel pipelined architecture for this method is possible.

As is well known, Fourier transform (FT) of a signal, x(t) is related to T of its shifted (or delayed if t represents time) version, x(t-4) by this equation:

where F{.cndot.} denotes Fourier transform. The phase of Fourier transform of the shifted signal contains the information about the amount of the shift r, which may easily be extracted. However, Discrete Cosine Transform (DCT) or its counterpart, Discrete Sine Transform, (DST) do not have any phase components as usually found in discrete Fourier transform (DFT), but DCT (or DST) coefficients of a shifted signal also carry this shift information. To facilitate explanation of the idea behind DXT-ME, one must first consider the case of one-dimensional discrete signals. Suppose that the signal {x.sub.1 (n); n.di-elect cons.{0, . . . , N-1}} is right shifted by an amount m (in our convention, a right shift mean that m>0) to generate another signal {x.sub.2 (n); n.di-elect cons.{0, . . . , N-1}}. The values of x.sub.1 (n) are all zeros outside the support region S(x.sub.1). Therefore, ##EQU8## implies that both signals have resemblance to each other except that the signal is shifted. It may be shown that, for k=1, . . . , N-1, ##EQU9## Here X.sub.2.sup.S and X.sub.2.sup.C are DST (DST-II) and DCT (DCT-II) of the second kind of x.sub.2 (n), respectively, whereas Z.sub.1.sup.S and Z.sub.1.sup.C are DST (DST-I) and DCT (DCT-I) of the first kind of x.sub.1 (n), respectively as defined as follows: ##EQU10## where ##EQU11## The displacement, m is embedded solely in the terms ##EQU12## which are called pseudo phases analogous to phases in Fourier transform of shifted signals. To find out m, we first solve equations (3) and (4) for the pseudo phases and then use the sinusoidal orthogonal principles as follows: ##EQU13## Here .delta.(n) is the discrete impulse function.

Indeed, if we replace ##EQU14## by the computed sine and cosine pseudo phase components, g.sub.m.sup.s (k) and g.sub.m.sup.c (k) in (9) and (10), both equations simply become IDST-II and IDCT operations on g.sub.m.sup.s (k) and g.sub.m.sup.c (k): ##EQU15## the notation g is used to distinguish the computed pseudo phase from the one in a noiseless situation (i.e., ##EQU16## A closer look at the right-hand side of equation (9) tells us that .delta.(m-n) and .delta.(m+n+1) have opposite signs. This property will help us detect the directions of a motion. If we perform an IDST-II operation on the pseudo phases found, then the observable window of the index space in the inverse DST domain will be limited to {0, . . . , N-1}. As illustrated in FIG. 2a, for a right-shift motion, one spike (generated by the positive .delta. function) is pointing upwards at the location n=m in the gray region (i.e., the observable index space), while the other .delta. point downwards at the n=-(m+1) outside the gray region. In contrary, for a left-shift motion, FIG. 2b, the negative spike at n=-(m+1)>0 falls in the gray region, but the positive .delta. function at n=m stays out of the observable index space. It may easily be seen that a positive peak value in the gray region implies a right shift and a negative one means a left shift. This enables one to determine from the sign of the peak value the direction of the movement of a signal.

The concept of pseudo phases plus the application of the sinusoidal orthogonal principles leads to a new approach to estimate the translational motion as depicted in FIG. 3:

1. Compute the DCT-I and DST-I coefficients of x.sub.1 (n) and the DCT-II and DST-II coefficients of x.sub.2 (n).

2. Compute the pseudo phase g.sub.m.sup.s (k) for k=, . . . , N by solving this equation: ##EQU17## 3. Feed the computed pseud phase, {g.sub.m.sup.s (k); k=1, . . . , N}, into an IDST-II decoder to produce an output {d(n); n=0, . . . , N-1}, and search for the peak value. The estimated displacement m may be found by: ##EQU18## where i.sub.p =arg max.sub.n .linevert split.d(n).linevert split. is the index at which the peak value is located.

In Step 1, the DCT and DST may be generated simultaneously with only 3N multipliers and the computation of DCT may be easily obtained from DCT-II with minimal overhead as will be shown later. In Step 2, if noise is absent and there is only purely translational motion, g.sub.m (k) will be equal to ##EQU19## The output d(n) will then be an impulse function in the observation window. This procedure is illustrated by two examples in FIGS. 3a to 3f and 3g to 3l with a randomly generated signal input at SNR=20 dB. From these two examples, it is obvious that the motion estimation is accurate even though strong noise is present.

The concept of how to extract shift values from the pseudo phases of one-dimensional signals, as explained above, may be extended to provide the basis of a new DCT-based motion estimation scheme for two-dimensional images, which is called DXT-ME Motion Estimation Scheme. Before going into details of this new method, let us confine the problem of motion estimation to the scenario in which an object moves translationally by m.sub.1 in X direction and n.sub.1 in Y direction as viewed on a camera plane and within the scope of a camera in a noiseless environment as shown in FIG. 4. Then one may extract the motion vector out of the two consecutive frames of the images of that moving object by making use of the sinusoidal orthogonal principles (9) and (10).

The method of this DCT-based motion estimation scheme is depicted in FIG. 5a. The previous frame x.sub.t-1 and the current frame x.sub.t are fed into 2D-DCT-II and 2D-DCT-1 coders, respectively. A 2D-DCT-II coder computes four coefficient, DCCTII, DCSTII, DCSTII and DSSTII, each of which is defined as a two-dimensional separable function formed by 1D-DCT/DST-II kernels: ##EQU20## or symbolically,

In the same fashion, the two-dimensional DCT coefficients of the first kind (2D-DCT-I) are calculated based on 1D-DCT/DST-I kernels: ##EQU21## or symbolically,

Similar to the one dimensional case, assuming that only translational motion is allowed, one may derive a set of equations to relate DCT coefficients of x.sub.t-1 (m,n) with those of x.sub.t (m,n) in the same way as in (3) and (4).

where N={1, . . . , N-1}, ##EQU22##

Here, Z.sub.t-1 (k,l).epsilon.R.sup.4.times.4 is the system matrix of DXT-ME at (k,l). At the boundaries of each block in the transform domain, the DCT coefficients of x.sub.t-1 (m,n) and x.sub.t (m,n) have a one-dimensional relationship as given below: ##EQU23##

In a two-dimensional space, an object may move in four possible directions: northeast (NE: m.sub.1 >0, n.sub.1 >0); northwest (NW: m.sub.1 <0, n.sub.1 <0); southeast (SE: m.sub.1 >0, n.sub.1 >0); and southwest (SW: m.sub.1 <0, n.sub.1 <0). As explained above, the orthogonal equation for the DST-II kernel in Equation (9) may be applied to the pseudo phase g.sub.m.sup.s (k) to determine the sign of m (i.e., the direction of the shift). In order to detect the signs of both m.sub.1 and n.sub.1 (or equivalently the direction of motion), from the observation in the one-dimensional case, it is necessary to compute the pseudo phases g.sub.m.sbsb.1.sub.n.sbsb.1.sup.CS (.cndot.,.cndot.) and g.sub.m.sbsb.1.sub.n.sbsb.1.sup.CS (.cndot.,.cndot.) so that the signs of m.sub.1 and n.sub.1 may be determined from g.sub.m.sbsb.1.sub.n.sbsb.1.sup.SC (.cndot., .cndot.) and g.sub.m.sbsb.1.sub.n.sbsb.1.sup.CS (.cndot.,.cndot.), respectively. By taking the block boundary equations (27)-(32) into consideration, we define the two pseudo phase functions as follows: ##EQU24##

These two pseudo phase functions pass through 2D-IDCT-II coders (IDCSTII and IDSCTII) to generate two functions, DCS(.cndot.,.cndot.) and DSC(.cndot.,.cndot.) in view of the orthogonal property of DCT-II and DST-II in (9) and (10): ##EQU25##

By the same argument as in the one-dimensional case, the 2D-IDCT-II coders limit the observable index space {(i, j):i, j=0, . . . , N-1} of DCS and DSC to the first quadrant of the entire index space as shown as gray regions in FIGS. 6a and 6d and 6e to 6h, respectively which depict (35) and (36). Similar to the one-dimensional case, if m.sub.1 is positive, the observable peak value of DSC(m, n) will be positive regardless of the sign n.sub.1 since DSC(m, n)=.delta.(m-m.sub.1).cndot.›.delta.(n-n.sub.1)+.delta.(n+n.sub.1 +1)! in the observable index space. Likewise, if m.sub.1 is negative, the observable peak value of DSC(m, n) will be negative because DSC(m, n)=.delta.(m+m.sub.1 +1).cndot.›.delta.(n-n.sub.1)+.delta.(n+n.sub.1 +1)! in the gray region. As a result, the sign of the observable peak value of DSC determines the sign of m.sub.1. The same reasoning may apply to DCS in the determination of the sign of n.sub.1. The estimated displacement, d=(m.sub.1,n.sub.1), may thus be found by locating the peaks of DCS and DSC over {0, . . . , N-1}.sup.2 or over an index range of interest, usually, .PHI.={0, . . . , n/2}.sup.2 for slow motion. How the peak signs determine the movement is summarized in Table I. Once the direction is found, .sup.d a may be estimated accordingly: ##EQU26##

Normally, these two peak indices are consistent, but in noisy circumstances they may not agree. In this case, an arbitration rule must be made to pick the best index (i.sub.D, j.sub.D) in terms of minimum non-peak-to-peak ratio (NPR): ##EQU28##

This index (i.sub.D, j.sub.D) will then be used to determine .sup.d by (37) and (38). Here NPR is defined as the ratio of the average of all absolute non-peak values to the absolute peak value. Thus 0.ltoreq.NPR.ltoreq.1, and for a pure impulse function, NPR=0. Such an approach to choose the best index among the two indices is found empirically to improve the noise immunity of this estimator.

In situations where slow motion is preferred, it is better to search the peak value in a zigzag way as widely used in DCT-based hybrid video coding. Starting from the index (0,0), zigzagly smay all the DCS or (DSC) values and mark the point as the new peak index if the value at that point (i, j) is larger than the current peak value by more than a present threshold .theta.:

In this way, large spurious spikes at the higher index points will not affect the performance of the estimator and thus improve its noise immunity further.

FIGS. 7a to 7h demonstrate this DXT-ME method. Images of a rectangularly-shaped moving object with arbitrary texture are generated as in FIG. 7a and corrupted by additive white Gaussian noise at SNR=10 dB as in FIGS. 7c and 7d. The resulted pseudo phase functions f and g, as well as DCS and DSC, are depicted in FIGS. 7e to 7h. Large peaks may be seen clearly in FIGS 7g and 7h on rough surfaces caused by noise in spite of noisy input images. The positions of these peaks give us an accurate motion estimate (5, -3).

The stability of this motion estimator depends upon the property of the determinant of the system matrix .linevert split.Z.sub.t-1 (k, l).linevert split. in (23). A zero or near-zero value of .linevert split.Z.sub.t-1 (k, l).linevert split. may jeopardize the performance. However, it may be shown analytically that this determinant will rarely be zero. Some algebraic manipulations on the determinant of Z.sub.t-1 (k, l) give us this close form of .linevert split.Z.sub.t-1 (k, l).linevert split.: ##EQU29## in which Z.sub.t-1 (k, l).linevert split.=0 implies that: Z.sub.t-1.sup.CS (k,l)=.+-.Z.sub.t-1.sup.SC (k,l), Z.sub.t-1.sup.CC (k,l)=.+-.Z.sub.t-1.sup.SS (k,l), Z.sub.t-1.sup.CC (k,l)Z.sub.t-1.sup.CS (k,l)=-Z.sub.t-1.sup.SC (k,l)Z.sub.t-1.sup.SS (k,l), Z.sub.t-1.sup.CC (k,l)Z.sub.t-1.sup.SC (k,l)=-Z.sub.t-1.sup.CS (k,l)Z.sub.t-1.sup.SS (k,l).

However, the last two equalities allow only two possible conditions:

Alternatively, in their explicitly compact forms, ##EQU30##

Here the minus signs in (47) and (48) correspond to the first condition in (45) and the plus signs correspond to (46). Satisfying either condition requires x.sub.t-1 (m, n)=0. Therefore, it is very unlikely that this determinant is zero and as such, the DXT-ME is a stable estimator. Even so, if Z.sub.t-1 (k, l)=0 really happens or Z.sub.t-1 (k, l) is less than a threshold, then we may let f(k, l)=g(k, l)=1, which is equivalent to the situation when x.sub.t-1 (m, n)=0. In this way, the catastrophical effect of computational precision of a certain implementation on the stability of DXT-ME will be kept to a minimum or even eliminated.

What if an object is moving in a uniformly bright background instead of a completely dark environment? It may be shown analytically and empirically that uniformly bright background introduces only very small spikes which does not affect the accuracy of the estimate. Suppose that {x.sub.t-1 (m, n)} and {x.sub.t (m, n)} are pixel values of two consecutive frames of an object displaced by (m.sub.1, n.sub.1) on a uniformly bright background. Then let y.sub.t (m, n) and y.sub.t-1 (m, n) be the pixel value of x.sub.t (m, n) and x.sub.t-1 (m, n) subtracted by the background pixel value c(c>0) respectively:

In this way, {x.sub.t-1 (m, n)} and {x.sub.t (m, n} may be considered as the images of an object moving in a dark environment. Denote Z.sub.t-1 (k, l) as the system matrix of the input image x.sub.t-1 and U.sub.t-1 (k, l) as that of y.sub.t-1 for k, l.epsilon.N. Also let x.sub.t (k,l)be the vector of the 2D-DCT-II coefficients of x.sub.t and y.sub.t (k,l) be the vector for y.sub.t. Applying the DXT-ME method to both situations, we have, for k, l.epsilon.N,

Here .O slashed..sub.m.sbsb.1.sub.n.sbsb.1 (k,l) is the vector of the computed pseudo phases for the case of dark background and thus ##EQU31##

Starting from the definition of each element in Z.sub.t-1 (k, l) and x.sub.t (k,l), we obtain

where D(k, l) is the system matrix with {d(m, n)=1, .A-inverted.m, n={0, . . . , N-1}} as input and c(k,l) is the vector of the 2D-DCT-II coefficients of d(m, n). Substituting (53) and (54) into (52) we get

Since c(k,l)=D(k,l).cndot..O slashed..sub.00 (k,l), (55) becomes

provided that .linevert split.Z.sub.t-1 (k, l).linevert split..noteq.0. Similar results may also be found at block boundaries. Referring to (24), we know that D(k, l) is composed of D.sup.CC (k, l), D.sup.CS (k, l), D.sup.SC (k, l), and D.sup.SS (k, l), each of which is a separable function made up by ##EQU32## From the above equations, we may see that D.sup.C (k)=D.sup.S (k)=0 if k is even, and for odd k>0, D.sup.c (k)=2/N while ##EQU33## Hence, D.sup.CC (k, l)=D.sup.CS (k, l)=D.sup.SC (k, l)=D.sup.SS (k, l)=0 if either k or l is even. As a result, .THETA..sub.m.sbsb.1.sub.n.sbsb.1 (k,l)=.phi..sub.m.sbsb.1.sub.n.sbsb.1 (k,l) if either k or l is even. For odd indices k and l, it is possible to find a constant s and a matrix N(k, l).epsilon.R.sup.4.times.4 such that U.sub.t-1 (k, l)=s›D(k, l)-N(k, l)! and .linevert split.N(k, l)D.sup.-1 (k, l).linevert split.<1 for .linevert split.D(k, l).linevert split..noteq.0. Therefore, for ##EQU34##

If we lump all the high-order terms of ##EQU35## in one term H(k, l), then ##EQU36## Usually, 0.ltoreq.c, s.ltoreq.255 for the maximum gray level equal to 255. Typically, s=1. For moderately large c, H(k, l) is very small. Define the subsampled version of the pseudo phase function .phi..sub.ab (k,l) as ##EQU37## Recall that a 2D-IDCT-II operation on .phi..sub.m.sbsb.1.sub.n.sbsb.1 (k,l) or .phi..sub.00 (k,l) produces .delta..sub.m.sbsb.1.sub.n.sbsb.1 or .delta..sub.00, respectively, where ##EQU38## Therefore, ##EQU39## where n is the noise term contributed from 2D-DCT-II{H(k,l)›.lambda..sub.00 -.lambda..sub.m.sbsb.1.sub.n.sbsb.1 !}. Because .lambda..sub.ab is equivalent to downsampling .phi..sub.ab in a 2D index space and it is known that downsampling produces in the transform domain mirror images of magnitude only one-fourth of the original and of sign depending on the transform function, we obtain ##EQU40## where diag(.cndot.) is the diagonal matrix of a vector and .cedilla.(i=1,2,3) is a vector consisting of .+-.1. A similar expression may also be established for 2D-DCT-II{.lambda..sub.00 }. In conclusion, ##EQU41## The above equation predicts the presence of a very small noise term n and several small spikes, E.sub.00 and E.sub.m.sbsb.1.sub.n.sbsb.1, of magnitude moderated by ##EQU42## which are much smaller than the displacement peak, as displayed in FIGS. 8c to 8f where n for the case of c=3 in (b) is observable but very small and may be regarded as noise whereas n is practically absent as in (c) when c=255.

The block diagram in FIG. 5a shows that a separate 2D-DCT-I is needed in addition to the standard DCT (2D-DCT-II). This is undesirable from the complexity viewpoint. However, this problem may be circumvented by considering the point-to-point relationship between 2D-DCT-I and 2D-DCT-II coefficients in the frequency domain for k, l.epsilon.N: ##EQU43## where X.sub.t-1.sup.CC, X.sub.t-1.sup.CS, X.sub.t-1.sup.SC, and X.sub.t-1.sup.ss are the 2D-DCT-II coefficients of the previous frame. Similar relation also exists for the coefficients at block boundaries. This observation results in the simple structure in FIG. 5b, where Block T is a coefficient transformation unit realizing (65).

If the DCT has computational complexity O.sub.dct, the overall complexity of DXT-ME is O(N.sup.2)+O.sub.dct with the complexity of each component summarized in Table II. The computational complexity of the pseudo phase computation component is only O(N.sup.2) for an N.times.N block and so is the unit to determine the displacement. For the computation of the pseudo phase functions f(.cndot.,.cndot.) in (33) and g(.cndot.,.cndot.) in (34), DSCT, DCST and DSST coefficients (regarded as DST coefficients) must be calculated in addition to DCCT coefficients (i.e., the usual 2D DCT). However, all of these coefficients may be generated with little overhead in the course of computing 2DDCT coefficients. As a matter of fact, a parallel and fully-pipelined 2d DCT lattice structure has been developed to generate 2D DCT coefficients at a cost of O(N) operations. This DCT coder computes DCT and DST coefficients dually due to its internal lattice architecture. These internally generated DST coefficients may be output to the DXT-ME module for pseudo phase computation. This same lattice structure may also be modified as a 2D IDCT which also has O(N) complexity. To summarize, the computation complexity of this DXT-ME is only O(N.sup.2), much lower than the O(N.sup.4) complexity of BMA-ME.

A closer look reveals that the operations of pseudo phase computation and coefficient transformation are performed independently at each point (k, l) in the transform domain and therefore are inherently highly parallel operations. Since most of the operations in the DXT-ME method involve mainly pseudo phase computation and coefficient transformation in addition to DCT and Inverse DCT operations which have been studied extensively, the DXT-ME method may easily be implemented on highly parallel array processors or dedicated circuits. This is very different from BMA-ME which requires shifting of pixels and summation of differences of pixel values and hence discourages parallel implementation.

For complicated video sequences in which objects may move across the border of blocks in non-uniform background, preprocessing may be employed to enhance the features of motion objects and avoid violation of the assumption made for DXT-ME before feeding the images into the DXT-ME motion estimator. Intuitively speaking, the DXT-ME method tries to match the features of any object on two consecutive frames to that any translation motion may be estimated regardless of the shape and texture of the object as long as these two frames contain significant energy level of the object features. Due to this feature matching property of the DXT-ME method, effective preprocessing will improve the performance of motion estimation if preprocessing may enhance the object features in the original sequence. In order to keep the computational complexity of the overall motion estimator low, the chosen preprocessing function must be simple but effective in the sense that unwanted features will not affect the accuracy of estimation. Our study found that both edge extraction and frame differentiation are simple and effective schemes for extraction of motion information.

Edges of an object may represent the object itself in motion estimation as its features and contain the information of motion without violating the assumption for DXT-ME. The other advantage of edge extraction is that any change in the illumination condition does not alter the edge information and, in turn, makes no false motion estimates by the DXT-ME method. Since we only intend to extract the main features of moving objects while keeping the overall complexity low, we employ a very simple edge detection by convolving horizontal and vertical Sobel operators of size 3.times.3 ##EQU44## with the image to obtain horizontal and vertical gradients respectively and then combine both gradients by taking the square root of the sum of the squares of both gradients. Edge detection provides us the features of moving objects but also the features of the background (stationary objects) which is undesirable. However, if the features of the background have smaller energy than those of moving objects within every block containing moving objects, then the background features will not affect the performance of DXT-ME. The computational complexity of this preprocessing step is only O(N.sup.2) and thus the overall computational complexity is still O(N.sup.2).

Frame differentiation generates an image of the difference of two consecutive frames. This frame differentiated image contains no background objects but the difference of moving objects between two frames. The DXT-ME estimator operates directly on this frame differentiated sequence to predict motion in the original sequence. The estimate will be good if the moving objects are moving constantly in one direction in three consecutive frames. For thirty frames per second, the standard NTSC frame rate, objects may usually be viewed as moving at a constant speed in three consecutive frames. Obviously this step also has only O(N.sup.2) computational complexity.

Alternatively, instead of using only one preprocessing function, we may employ several simple operators in the preprocessing step to extract features of images as shown in FIG. 9a, in which four DXT-ME estimators generate four candidate estimates of which one may be chosen as the final estimated displacement based upon either the mean squared error per pixel (MSE) or the mean of absolute differences per pixel (MAD) criteria.

Preferably, a simple decision rule similar to the one used in the MPEG-standard 21, as depicted in FIG. 9b is used to choose among the DXT-ME estimate and no motion. This simplified extended DXT-ME motion estimator works very well as will be shown in the next section.

Above, it was mentioned that peaks of DSC and DCS are searched over a fixed index range of interest .PHI.={0, . . . , N/2}.sup.2. However, if we follow the partitioning approach used in BMA-ME, then we may dynamically adjust .PHI.. At first, partition the whole current frame into bs.times.bs nonoverlapping reference blocks shown as the shaded area in FIG. 10a. Each reference block is associated with a larger search area (of size sa) in the previous frame (the dotted region in the same figure) in the same way as for BMA-ME. From the position of a reference block and its associated search area, a search range D={(u, v):-u.sub.1 .ltoreq.u.ltoreq.u.sub.2, -v.sub.1 .ltoreq.v.ltoreq.v.sub.2 } may then be determined as in FIG. 10b. Differing from BMA-ME, DXT-ME requires that the reference block size and the search area size must be equal. Thus, instead of using the reference block, we use the block of the same size and position in the current frame as the search area of the previous frame. The peak values of DSC and DCS are searched in a zigzag way as described in Section III over this index range .PHI.={0, . . . , max(u.sub.2, u.sub.1 -1)}.times.{0, . . . , max(v.sub.2, v.sub.1 -1)}. In addition to the requirement that the new peak value must be larger than the current peak value by a preset threshold, it is necessary to examine if the motion estimate determined by the new peak index lies in the search region D. Since search areas overlap on one another, the SE-DXT-ME architecture utilizing this approach is called Overlapping SE-DXT-ME.

Simulations have been performed on a number of video sequences with different characteristics. The performance of the DXT-ME scheme is evaluated in terms of MSE (mean squared error per pel) defined as ##EQU45## and compared with that of the full search block matching method (BMA-ME), which minimizes the MSE function over the whole search area: ##EQU46##

The MSE values of two consecutive frames without motion compensation (DIF) are also computed for comparison. The MSE value of DXT-ME is expected to be upper bounded by the MSE value without motion compensation.

To test the performance of DXT-ME on noisy images, an image of a small car (SCAR.1) is manually shifted to produce the second frame (SCAR.2) with a known displacement and additive Gaussian noise is added to attain a desired signal-to-ratio (SNR) level. Since the object (small car) moves within the boundary of the frame in a completely darkened background, no preprocessing is required. As may be seen in FIGS. 11a to 11d and 11i to 11l, DXT-ME is performed on the whole image of block size 64.times.64 and estimates the motion correctly at SNR level even down to 0 dB, whereas BMA-ME as shown in FIGS. 11e to 11h and 11m to 11p, produces some wrong motion estimates for boundary blocks and blocks of low signal energy. The values of MAD also indicate better overall performance of DXT-ME over BMA-ME for this sequence. Furthermore, DXT-ME may perform on the whole frame while BMA-ME needs division of the frame into sub-blocks due to the requirement of larger search areas than reference blocks. This is one of the reasons that BMA-ME does not work as well as DXT-ME because smaller block size makes BMA-ME more susceptible to noise and operation of DXT-ME on the whole frame instead of on smaller blocks lends itself to better noise immunity. Even though the Culmen filtering approach 30 may also estimate velocity accurately for a sequence of noisy images, it requires iterative complicated computations while DXT-ME may estimate motion based upon two consecutive frames in one step, requiring low-complexity computations.

The first sequence is the "Flower Garden" (FG) sequence where the camera is moving before a big tree and a flower garden in front of a house. Each frame has 720.times.486 pixels. Simple preprocessing is applied to this sequence: edge extraction on frame differentiation. Simulation of SE-DXT-ME of block size 8, 16, 32 as well as 64, was performed on 100 frames of this sequence. The results were compared with those of BMA-ME of block sizes 8, 16, 32 with the search area being twice as large as the block as shown in FIGS. 13a and 13b.

As may be seen in FIG. 12b, the edge extracted frames contain significant features of moving objects in the original frames so that DXT-ME may estimate the movement of objects based upon the information provided by the edge extracted frames. Because the camera is moving at a constant speed in a one direction, the moving objects occupy almost the whole scene. Therefore, the background features do not interfere with the operation of DXT-ME. The frame differentiated images of the "Flower Garden" sequence, one of which is shown in FIG. 12c, have the residual energy strong enough for DXT-ME to estimate the motion directly on this frame differentiated sequence due to the constant movement of the camera.

Observable in FIGS. 13a and 13b, a large reference block hampers the performance of BMA-ME indicated by high MSE values whereas increasing block size may boost up the performance of DXT-ME with smaller MSE values. The reason is that a block of larger size for DXT-ME contains more features of objects, which enables DXT-ME to find a better estimate due to its feature matching property, and also a block of larger size means a larger search area because the size of a search area is the same as the block size for DXT-ME. As a matter of fact, BMA-ME is supposed to perform better than DXT-ME if both methods use the same block size because BMA-ME requires a larger search area and thus BMA-ME has more information available before processing than DXT-ME. Therefore, it is not faire to compare BMA-ME with DXT for the same block size. Instead, it is more reasonable to compare BMA-ME with DXT-ME of the block size equal to the size of the search area of BMA-ME. As shown in FIGS. 13a and 13b, the MSE values for DXT-ME of block size 32 (DXT32) preprocessed by either edge extraction or frame differentiation are comparable with or even better than those for BMA-ME of block size 32 and search size 64.

If the sequence is shrunk to half size (352.times.224 pixels) and forms the small "Flower Garden" sequence, then each block will capture more features than a block of the same size in the original "Flower Garden" sequence. The simulation results are plotted in FIGS. 14a and 14c for SE-DXT-ME, and in FIGS. 14b and 14d for Overlapping SE-DXT-ME discussed in Section IV. As expected, the MSE values for SE-DXT-ME of block size 16 now become as small as those for BMA-ME of block size 32 as shown. Some points of the curve for SE-DXT-ME of block size 32 (DXT32) with either preprocessing step are below the points of BMA-ME of search size 32 (BKM16:32). If the overlapping approach is adopted in determining the search region, the results suggest that both Overlapping SE-DXT-ME and BMA-ME have comparable performance for either preprocessing method. The fact that smaller frame size of the contents improves the performance of DXT-ME in terms of smaller MSE values recommends the hierarchical motion estimation technique combined with DXT-ME (SE-DXT-ME or Overlapping SE-DXT-ME).

Another simulation is done on the "Infrared Car" sequence which has the frame size 192.times.224 and one obvious moving object, the car moving along the curved road towards the camera fixed on the ground. In FIG. 15b, the feature of both the car and the background are captured in the edge extracted frames. Even though the background features are not desirable, the simulation for SE-DXT-ME of various block sizes shows that the estimates of SE-DXT-ME produce low MSE values compared to the result of BMA-ME, especially in certain frames such as the 6th to 13th frames shown in FIG. 16a. This may be explained by taking a closer look at the edge extracted frame in FIG. 15b where the background features, such as the trees and poles, are far away from the car, especially in the 6th to 13th frames. For the first few frames, the features of the roadside behind the car mix with the features of the car and affect the performance of DXT-ME but then the car moves away from the roadside towards the camera so that the car features are isolated from the background features and so DXT-ME may estimate motion more accurately. As to the frame differentiated images as shown in FIG. 15c, the residual energy of the moving car is completely separated from the rest of the scene in most of the preprocessed frames and, therefore, lower MSE values are obtained with this preprocessing function than with edge extraction. As may be seen in FIGS. 16a and 16c, the MSE value of SE-DXT-ME of block size 32 at Frame 9 is smaller than that of BMA-ME of search area 32. In other words, SE-DXT-ME is better than BMA-ME in this particular case. However, the results of Overlapping SE-DXT-ME in FIGS. 16b and 16d show little improvement over SE-DXT-ME with fixed .PHI.={0, . . . , N/2}.sup.2, in contrary to the large gain in the performance of Overlapping SE-DXT-ME over SE-DXT-ME on the small "Flower Garden" sequence.

Simulation is also performed on the "Miss America" sequence, in which a lady is talking to the camera. Each frame has 352.times.288 pixels. This sequence has only translational motion of the head and shoulders, but mainly the mouth and eyes open and close. This makes the task of motion estimation difficult for this sequence but the DXT-ME method may still perform reasonably well compared to the BMA-ME method, as may be found in FIGS. 17a and 17c. Especially for the Overlapping SE-DXT-ME method with block size 32, the MSE values are very close to those of BMA-ME of different block size as shown in FIGS. 17b and 17d.

The last sequence is the small "Table Tennis" which has the frame size 352.times.224 pixels. In this sequence, the first twenty-three frames contain solely a bouncing ball with rough texture of the wall behind. Thus, the edge extracted frames capture a lot of background features mixed with the ball features. This influences negatively the estimation of DXT-ME. However, the frame differentiated images have no such background features and, as a result, Overlapping SE-DXT-ME of block size 32 (DXT32:64) with frame differences as input performs much better than Overlapping SE-DXT-ME after edge extraction (see FIG. 18b) even though the ball is not moving in a constant speed and its residual energy after frame subtraction is weak. After the 23rd frame, the camera is zooming out quickly, making no method predict well. Then the zooming action slows down. In this situation, the MSE values of SE-DXT-ME go down suddenly to as low as those of BMA-ME as shown in FIGS. 18a and 18b.

To further improve the coder gain, motion estimation with subpixel accuracy is essential because movements in a video sequence are not necessarily multiples of the sampling grid distance. Most currently existing subpixel motion estimation schemes estimate displacements from interpolated images. However, interpolation not only increases the complexity and data flow of a coder but also may adversely affect the accuracy of motion estimates from interpolated images. It is more desirable that subpixel accuracy of motion estimates can be obtained without interpolating the images. Based on the sinusoidal orthogonal principles at the subpixel level, it is possible to estimate subpixel motion without interpolation to a low computational cost in the DCT domain. Studies of pseudo phases and the sinusoidal orthogonal principles at the subpixel level enable us to provide half-pel or even quarter-pel DCT-based motion estimation schemes with minor modification of the algorithm proposed here.

Although the present invention has been fully described in connection with the preferred embodiment thereof with reference to the accompanying drawings, it is to be noted that various changes and modifications are apparent to those skilled in the art. Such changes and modifications are to be understood as included within the scope of the present invention as defined by the appended claims, unless they depart therefrom.