Estimation of Subpixel Motion Using Bispectrum

Abstract

Motion estimation techniques are widely used in today's video processing systems. The frequently used techniques are frequency-domain motion estimation methods, most notably phase correlation (PC). If the image frames are corrupted by Gaussian noises, then cross-correlation and related techniques do not work well. In this paper, however, we have studied this topic from a viewpoint different from the above. Our scheme is based on the bispectrum method for sub-pixel motion estimation of noisy image sequences. Experimental results show that our proposed method performs significantly better than PC technique.

1. Introduction

Image frames are generated by scanning a scene several times a second where each frame, generally, consists of two regions. One region, referred to as stationary background, is virtually the same as the previous frame. The other region, referred to as moving data, has moved with respect to the previous frame. Estimating the motion between image frames has long been a problem of interest in areas such as video compression, robot vision, and biomedical engineering [1].

Many motion estimation schemes have been developed. They can be classified into spatial-domain and frequency-domain approaches. The spatial-domain algorithms consist of matching algorithms and gradient-based algorithms. The frequency domain algorithms consist of phase correlation algorithms, wavelet transform-based algorithms, and DCT-based algorithms [2]. The vast majority of these algorithms consider noise-free data, although in [3] the displacement vector is estimated from noisy data using the generalized maximum likelihood criterion. If the image frames are corrupted by Gaussian noises, then cross-correlation and related techniques do not work well. In this circumstance, higher-order spectra (HOS) in general and the bispectrum in particular have recently been widely used as an important tool for signal processing. The classical methods based on the power spectrum are now being effectively superseded by the bispectral ones due to some definite disadvantages of the former. These include the inability to identify systems fed by non-Gaussian noise (NGN) inputs and nonminimum phase (NMP) systems, and by identification of system nonlinearity. In these cases, the autocorrelation-based methods offer no answer. Out of all these, the identifiability of NMP systems has received the maximum attention from researchers.

HOS-based methods have already been proposed to estimate motion between image frames [4–9]. In [6], the displacement vector is obtained by maximizing a third-order statistics criterion. In [8], the global motion parameters obtained by a new region recursive algorithm. In [4, 5], several algorithms are developed based on a parametric cumulant method, a cumulant-matching method, and a mean-kurtosis error criterion. In this correspondence, a novel algorithm for the detection of motion vectors in video sequences is proposed. The algorithm uses bispectrum method for subpixel motion estimation of noisy image sequences to obtain a measure of content similarity for temporally adjacent frames and responds very well to scene motion vectors.

2. Bispectrum-Based Image Motion Estimation

The problem of motion estimation can be stated as follows: “given an image sequence, compute a representation of the motion field that best aligns pixels in one frame of the sequence with those in the next” [9]. This is formulated as 𝑔𝑔_𝑘−1(𝑥,𝑦)=𝑓_𝑘−1(𝑥,𝑦)+𝑛_𝑘−1(𝑥,𝑦),_𝑘(𝑥,𝑦)==𝑓𝑓_𝑘(𝑥_𝑘−1,𝑦)+𝑛_𝑘(𝑥,𝑦)=𝑔_𝑘−1𝑥−𝑑^𝑘𝑥,𝑦−𝑑_𝑘𝑦+𝑛_𝑘(𝑥,𝑦)𝑥−𝑑^𝑘𝑥,𝑦−𝑑_𝑘𝑦+𝑤_𝑘(𝑥,𝑦),(1)where 𝑤_𝑘(𝑥,𝑦)=𝑛_𝑘(𝑥,𝑦)−𝑛_𝑘−1(𝑥−𝑑_𝑘𝑥,𝑦−𝑑_𝑘𝑦); (𝑥,𝑦) denotes spatial image position of a point; 𝑔_𝑘(𝑥,𝑦) and 𝑔_𝑘−1(𝑥,𝑦) are observed image intensities at instant 𝑘 and 𝑘−1, respectively; 𝑓_𝑘(𝑥,𝑦) and 𝑓_𝑘−1(𝑥,𝑦) are the noise-free frames; 𝑛_𝑘(𝑥,𝑦) and 𝑛_𝑘−1(𝑥,𝑦) are assumed to be spatially and temporally stationary, zero-mean image Gaussian noise sequences with unknown covariance; (𝑑_𝑘𝑥,𝑑_𝑘𝑦) is the displacement vector of the object during the time interval [𝑘−1,𝑘].

The third-order autocumulants and cross-cumulants of a zero-mean from 2D random field 𝑔_𝑘(𝑥,𝑦) are defined, respectively, as follows:=𝐸𝑔_𝑘𝐶^{𝑔𝑘𝑔𝑘𝑔𝑘}₃𝑟¹,𝑟₂;𝑠₁,𝑠₂(𝑥,𝑦)𝑔_𝑘𝑥+𝑟¹,𝑦+𝑟₂𝐶⋅𝑔^𝑘𝑥+𝑠¹,𝑦+𝑠₂_{𝑔𝑘𝑔𝑘−1𝑔𝑘3},=𝐸𝑔_𝑘𝑟¹,𝑟₂;𝑠₁,𝑠₂(𝑥,𝑦)𝑔_𝑘−1𝑥+𝑟¹,𝑦+𝑟₂⋅𝑔^𝑘𝑥+𝑠¹,𝑦+𝑠₂,(2)where 𝐸{⋅} represents the expectation operator.

The discrete bispectrum of the frame 𝑔_𝑘(𝑥,𝑦) is defined as follows: =ℱ𝐶𝐵^{𝑔𝑘𝑔𝑘𝑔𝑘}_{3𝑔𝑘𝑔𝑘𝑔𝑘3}𝑢¹,𝑢₂;𝑣₁,𝑣₂=𝐺_𝑔𝑘𝑟¹,𝑟₂;𝑠₁,𝑠₂𝑢¹,𝑢₂𝐺^𝑔_𝑘𝑣¹,𝑣₂𝐺^∗𝑔_𝑘𝑢¹+𝑣₁,𝑢₂+𝑣₂,(3)where ℱ denotes the Fourier transform operation; 𝐺_𝑔𝑘(𝑢) corresponds to the DFT of the frame 𝑔_𝑘(𝑥,𝑦); ∗ indicates the complex conjugate; (𝑢₁,𝑢₂) and (𝑣₁,𝑣₂) are the frequency co-ordinates for the 2D Fourier transform.

Due to the shift-invariance of the bispectrum=𝐵_{𝑔𝑘−1𝑔𝑘−1𝑔𝑘−13}𝐵^{𝑔𝑘𝑔𝑘𝑔𝑘}₃𝑢¹,𝑢₂;𝑣₁,𝑣₂𝑢¹,𝑢₂;𝑣₁,𝑣₂+𝐵^{𝑤𝑘𝑤𝑘𝑤𝑘3}𝑢¹,𝑢₂;𝑣₁,𝑣₂,(4)where 𝐵^{𝑔𝑘−1𝑔𝑘−1𝑔𝑘−1}₃(𝑢₁,𝑢₂;𝑣₁,𝑣₂) and 𝐵^{𝑤𝑘𝑤𝑘𝑤𝑘}₃(𝑢₁,𝑢₂;𝑣₁,𝑣₂) denote bispectrum of the frame 𝑔_𝑘−1(𝑥,𝑦) and noise, respectively. If the probability density function of the noise is symmetrical, that is, 𝑝(𝑛)=𝑝(−𝑛), or at least not skewed, that is, ∫𝑛³𝑝(𝑛)𝑑𝑛=0, then the term 𝐵^{𝑤𝑘𝑤𝑘𝑤𝑘}₃(𝑢₁,𝑢₂;𝑣₁,𝑣₂) is negligible which renders the triple-correlation very effective in detecting a signal embedded in noise [10]. Then𝐵^{𝑔𝑘𝑔𝑘𝑔𝑘}₃𝑢¹,𝑢₂;𝑣₁,𝑣₂=𝐵^{𝑔𝑘−1𝑔𝑘−1𝑔𝑘−13}𝑢¹,𝑢₂;𝑣₁,𝑣₂.(5)The cross-bispectrum is defined as follows:=ℱ𝐶𝐵^{𝑔𝑘𝑔𝑘−1𝑔𝑘}_{3𝑔𝑘𝑔𝑘−1𝑔𝑘3}𝑢¹,𝑢₂;𝑣₁,𝑣₂=𝐺_𝑔𝑘𝑟¹,𝑟₂;𝑠₁,𝑠₂𝑢¹,𝑢₂𝐺^𝑔_𝑘−1𝑣¹,𝑣₂𝐺^∗𝑔_𝑘𝑢¹+𝑣₁,𝑢₂+𝑣₂.(6)Then,=𝐵_{𝑔𝑘−1𝑔𝑘−1𝑔𝑘−13}𝐵^{𝑔𝑘𝑔𝑘−1𝑔𝑘}₃𝑢¹,𝑢₂;𝑣₁,𝑣₂𝑢¹,𝑢₂;𝑣₁,𝑣₂⋅𝑒^{−𝑗2𝜋(𝑣1𝑑𝑘𝑥+𝑣2𝑑𝑘𝑦)}.(7)Using the relation in (5), we can transform (7) as=𝐵_{𝑔𝑘𝑔𝑘𝑔𝑘3}𝐵^{𝑔𝑘𝑔𝑘−1𝑔𝑘}₃𝑢¹,𝑢₂;𝑣₁,𝑣₂𝑢¹,𝑢₂;𝑣₁,𝑣₂⋅𝑒^{−𝑗2𝜋(𝑣1𝑑𝑘𝑥+𝑣2𝑑𝑘𝑦)}.(8)

As we can see from (3), the bispectrum has two vector arguments containing totally four scalar frequency variables. Assuming that 𝐺_𝑔𝑘(𝑢₁,𝑢₂) is an N-by-N discrete Fourier transform of 𝑔_𝑘(𝑥,𝑦), the bispectrum becomes a four-dimensional N-by-N-by-N-by-N matrix. It is therefore not practical to evaluate the whole bispectrum. A better solution is to take 2D slices of the 4D spectrum. There are basically various ways of defining these slices, but we will only consider the case where𝐵𝐵^{𝑔𝑘𝑔𝑘𝑔𝑘}_3,𝑙𝑢¹,𝑢₂=𝐵^{𝑔𝑘𝑔𝑘𝑔𝑘3}𝑢¹,𝑢₂;𝑙𝑢₁,𝑙𝑢₂∀𝑙∈𝑅,_{𝑔𝑘𝑔𝑘−1𝑔𝑘3,𝑙}𝑢¹,𝑢₂=𝐵^{𝑔𝑘𝑔𝑘−1𝑔𝑘3}𝑢¹,𝑢₂;𝑙𝑢₁,𝑙𝑢₂∀𝑙∈𝑅.(9)

Although we have now taken only a small portion of the whole spectrum [11], it can be shown that the motion vector is still possible and no essential information has been lost.

Thus, the third-order hologram, ℎ_𝑙(𝑟₁,𝑟₂), is defined byℎ_𝑙𝑟¹,𝑟₂=ℱ⁻¹[𝐵^{𝑔𝑘𝑔𝑘−1𝑔𝑘3,𝑙}𝑢¹,𝑢₂|||𝐵_{𝑔𝑘𝑔𝑘𝑔𝑘3,𝑙}𝑢¹,𝑢₂|||]=𝛿𝑟¹−𝑑_𝑘𝑥,𝑟₂−𝑑_𝑘𝑦.(10)

As a result, by finding the location of the pulse in (10) we are able to tell the displacement, which is the motion vector. Since third-order statistics are used, the method is insensitive (in theory) to both spatially and temporally corrupted by noise which is symmetrically distributed (e.g., Gaussian). In practice, the motion vector is not an impulse; hence, we estimate (𝑑_𝑘𝑥,𝑑_𝑘𝑦) as the index (𝑟₁,𝑟₂), which maximizes |ℎ_𝑙(𝑟₁,𝑟₂)|.

The co-ordinates (𝑟_1𝑚,𝑟_2𝑚) of the maximum of the real-valued array ℎ_𝑙(𝑟₁,𝑟₂) can be used as an estimate of the horizontal and vertical components of motion between 𝑔_𝑘(𝑥,𝑦) and 𝑔_𝑘−1(𝑥,𝑦) as follows:𝑟^1𝑚,𝑟_2𝑚=argmaxℎ^𝑙𝑟¹,𝑟₂.(11)

3. Subpixel Accuracy

Subpixel performance is a critical element of the proposed algorithm. With reference to our previously published work [12, 13], we are introducing a number of important new features, which improve the accuracy of the motion estimates.

Subpixel accuracy of motion measurements is obtained dy variable-separable fitting performed in the neighborhood of the maximum using one-dimensional quadratic function. Using the notation in (11) above, prototype functions are fitted to the tripletsℎℎ^𝑙𝑟^1𝑚−1,𝑟_2𝑚,ℎ^𝑙𝑟^1𝑚,𝑟_2𝑚,ℎ^𝑙𝑟^1𝑚+1,𝑟_2𝑚,(12)_𝑙𝑟^1𝑚,𝑟_2𝑚−1,ℎ_𝑙𝑟^1𝑚,𝑟_2𝑚,ℎ^𝑙𝑟^1𝑚,𝑟_2𝑚+1,(13)(𝑑_𝑘𝑥,𝑑_𝑘𝑦)that is, the maximum peak of the phase correlation surface and its two neighboring values on either side, vertically and horizontally.

The location of the maximum of the fitted function provides the required subpixel motion estimate 𝑑_𝑘𝑥. Fitting a parabolic function horizontally to the data triplet (12) yields a closed-from solution for the horizontal component of the motion estimate 𝑑_𝑘𝑥=ℎ^𝑙𝑟^1𝑚+1,𝑟_2𝑚−ℎ^𝑙𝑟^1𝑚−1,𝑟_2𝑚2𝐻,(14) as follows:𝐻=ℎ_𝑙𝑟^1𝑚+1,𝑟_2𝑚−2ℎ^𝑙𝑟^1𝑚,𝑟_2𝑚+ℎ^𝑙𝑟^1𝑚−1,𝑟_2𝑚.where 𝑑_𝑘𝑦

The fractional part of the vertical component can be obtained in a similar way using (13) instead of (12).

Finally, the horizontal and vertical components of the subpixel accurate motion estimate are obtained by computing the location of the maxima of each of the above fitted quadratics.

In [14], it is shown that half-pixel accuracy motion vectors leads to a very significant improvement when compared to one-pixel accuracy, where as a higher precision results in negligible changes. Therefore, a half-pixel accuracy was chosen in our simulations.

4. Experimental Results

To prove the feasibility of the proposed method, we compared it to a PC technique implemented in a similar manner as our approach. In this section, we examine a few examples and compare the performance, efficiency, and complexity of the two methods. In our experiments, we used the well-known test sequences: foreman (176 pixels by 144 lines), mother-daughter and Stefan (352 pixels by 288 lines), table tennis (352 pixels by 240 lines). Although the original sequences are in color, only the luminance (brightness) component is used to estimate the motion vectors.

To assess the performances of the different motion estimation techniques, the following comparisons were made. First, the subjective quality of the estimated motion field was evaluated, showing the capability of the algorithm to estimate the true motion in the scene. Second, the PSNR of motion compensated was measured, giving insight about the quality of the prediction. Results obtained using the foreman, mother-daughter and Stefan sequences are shown in Figure 1. All image sequences are degraded with additive zero-mean Gaussian noise to a signal-to-noise ratio (SNR) of 10 dB. Our results demonstrate that the proposed method outperforms PC, achieving higher precision and a significantly smaller corresponding measurement error. This confirms the motion that the proposed technique of an image is a superior feature selector utilizing the portions of the image spectrum most likely to contribute to reliable motion estimation.

Figure 1 

PSNR obtained for noisy sequences (SNR = 10 dB).

The ability of the bispectrum method to accurately estimate the displacement vector field from a degraded sequence is demonstrated in Figure 2. In this Figure, the estimated motion vector fields for the mother-daughter sequence using the two aforementioned motion estimation methods. The motion vectors estimated between the frames 126 and 127 are shown for the mother-daughter sequence. We can see that the estimates from the PC seem very random, but the bispectrum technique gives better results, producing the same motion vectors. Thus, the motion fields estimated by the “our approach” tend to be very smooth due to the smoothness constraint. Because of the noise-resistant property of the bispectrum, it produces more reliable estimates. Therefore, the proposed method motion estimation results globally in motion fields more representative of the true motion in the scene.

(a) PC

(b) Bispectrum

Figure 2 

Motion field for the mother-daughter sequence in the presence of noise.

In terms of motion compensated images, from mother-daughter sequence, we observe better compensated images by the proposed method. We also observe that the motion compensated images for the “our method” are much closer to the original images. Thus, the “our scheme” is able to measure the motion vector more accurately and is more robust in general. Overall, the bispectrum typically offers better visual quality images than the PC method. Figure 3 gives examples of motion compensated images.

(a) Original image

(b) PC

(c) Bispectrum

Figure 3 

Prediction for frame 3 of the mother-daughter.

Comparisons of the PC and bispectrum methods indicate that the bispectrum is a robust technique for motion vector. Results of these comparisons are shown for different noise levels and video sequences. Additive Gaussian noise (AGN) was added to image sequences with an input SNR varying from 15 dB to 28 dB. Figure 4 shows the PSNR of the motion compensated prediction error for noise power of 15 dB for table tennis sequence, confirming that our scheme is consistently more immune to noise. Experiments with other levels of noise power demonstrated that similar performance gains are achievable. Similar experiments were performed for the other sequences. The results further confirm that the “our method” consistently outperforms PC. In terms of complexity, this is measured by the computation time. All the computations are performed on Intel(R) Pentium(R)D CPU 3.4 GHz with Windows XP. The two algorithms have been implemented using a prototype written in MATLAB 6.5 R13. The comparison between the “our method” and the PC confirmed that the two methods have the complexity on the same order. This is shown in Table 1.

Table 1 

The comparison between two methods for the computation time.

Figure 4 

PSNR versus frame number for motion compensated prediction of the table tennis sequence.

5. Conclusion

In this paper, the bispectrum method for subpixel motion estimation of noisy image sequences in frequency-domain was presented. The “our” proposed method provides an advantage over the PC algorithm in the presence of AGN. With “our method,” the displacement vector field is smoother, providing a more accurate measure of object motion. At relatively low noise levels, the bispectrum performance is comparable to its performance in the noise-free environment. At high noise levels SNR around 10 dB, the PC fails, yet even under these extreme conditions, the bispectrum provides improvement in performance over the PC algorithm. In addition to its PSNR performance, the bispectrum also yields smooth motion fields.