This is how the DFT may be computed efficiently.

has to be evaluated for

It is possible to calculate the DFT more efficiently than this, using the

We shall assume for simplicity that

If we define to be the root of unity given by , and set

This can be split apart into two separate sums of alternate terms from the original sum,

Now, since the square of a root of unity is an root of unity, we have that

and hence

If we call the two sums demarcated above and respectively, then we have

Note that each of and for is in itself a discrete Fourier transform over

How does this help us?

Well

and we can also write

Thus, we can compute an

- The first half of
*F*(*u*) for can be found from Eqn. 28, - The second half for can be found simply be reusing the same terms differently as shown by Eqn. 30.
- This is obviously a divide and conquer method.

To show how many operations this requires, let

the first term on the right hand side coming from the two transforms of half the original size, and the second term coming from the multiplications of by . Induction can be used to prove that

A similar argument can also be applied to the number of additions required, to show that the algorithm as a whole takes time .

Also Note that the same algorithm can be used with a little modification to perform the inverse DFT too. Going back to the definitions of the DFT and its inverse,

and

If we take the complex conjugate of the second equation, we have that

This now looks (apart from a factor of 1/

Thus to compute an inverse DFT,

- take the conjugate of the Fourier space data,
- put conjugate through a forward DFT algorithm,
- take the conjugate of the result, at the same time multiplying each value by
*N*.

The same fast Fourier transform algorithm can be used -- applying the separability property of the 2D transform.

Rewrite the 2D DFT as

The right hand sum is basically just a one-dimensional DFT if

So we can compute a two-dimensional DFT by

- performing a one-dimensional DFT for each value of
*x*,*i.e.*for each column of*f*(*x*,*y*), then - performing a one-dimensional DFT in the opposite direction (for each row) on the resulting values.

This requires a total of 2

Although it is well known that cross correlation can be efficiently implemented
in the transform domain, the normalized form of cross correlation preferred for
feature matching applications does not have a simple frequency domain expression.
Normalized cross correlation has been computed in the spatial domain for this reason.
This short paper shows that unnormalized cross correlation can be efficiently normalized
using precomputing integrals of the image and image
over the search window.

The correlation between two signals (cross correlation) is a standard approach to feature detection [6, 7] as well as a component of more sophisticated techniques (e.g. [3]).

Textbook presentations of correlation describe the convolution theorem and the attendant possibility of efficiently computing correlation in the frequency domain using the fast Fourier transform. Unfortunately the normalized form of correlation (correlation coefficient) preferred in template matching does not have a correspondingly simple and efficient frequency domain expression. For this reason normalized cross-correlation has been computed in the spatial domain (e.g., [7], p. 585). Due to the computational cost of spatial domain convolution, several inexact but fast spatial domain matching methods have also been developed [2]. This paper describes a recently introduced algorithm [10] for obtaining normalized cross correlation from transform domain convolution. The new algorithm in some cases provides an order of magnitude speedup over spatial domain computation of normalized cross correlation (Section 5).

Since we are presenting a version of a familiar and widely used algorithm no attempt will be made to survey the literature on selection of features, whitening, fast convolution techniques, extensions, alternate techniques, or applications. The literature on these topics can be approached through introductory texts and handbooks [16, 7, 13] and recent papers such as [1, 19]. Nevertheless, due to the variety of feature tracking schemes that have been advocated it may be necessary to establish that normalized cross-correlation remains a viable choice for some if not all applications. This is done in section 3.

In order to make the paper self contained, section 2 describes normalized cross-correlation and section 4 briefly reviews transform domain and other fast convolution approaches and the phase correlation technique. These sections can be skipped by most readers. Section 5 describes how normalized cross-correlation can be obtained from a transform domain computation of correlation. Section 6 presents performance results.

(where

the term is constant. If the term is approximately constant then the remaining cross-correlation term

is a measure of the similarity between the image and the feature.

There are several disadvantages to using (1) for template matching:

- If the image energy varies with position, matching using (1) can fail. For example, the correlation between the feature and an exactly matching region in the image may be less than the correlation between the feature and a bright spot.
- The range of
*c*(*u*,*v*) is dependent on the size of the feature. - Eq. (1) is not invariant to changes in image amplitude such as those caused by changing lighting conditions across the image sequence.

The

where is the mean of the feature and is the mean of

It is clear that normalized cross-correlation (NCC) is not the ideal approach to feature tracking since it is not invariant with respect to imaging scale, rotation, and perspective distortions. These limitations have been addressed in various schemes including some that incorporate NCC as a component. This paper does not advocate the choice of NCC over alternate approaches. Rather, the following discussion will point out some of the issues involved in various approaches to feature tracking, and will conclude that NCC is a reasonable choice for some applications.

The SSDA algorithm is simple and provides a significant speedup over spatial domain cross-correlation. It has the disadvantage that it does not guarantee finding the maximum of the correlation surface. SSDA performs well when the correlation surface has shallow slopes and broad maxima. While this condition is probably satisfied in many applications, it is evident that images containing arrays of objects (pebbles, bricks, other textures) can generate multiple narrow extrema in the correlation surface and thus mislead an SSDA approach. A secondary disadvantage of SSDA is that it has parameters that need to determined (the number of terms used to form an estimate of the correlation coefficient, and the early termination threshold on this estimate).

Successful gradient descent search requires that the interframe translation be less than the radius of the basin surrounding the minimum of the matching error surface. This condition may be satisfied in many applications. Images sequences from hand-held cameras can violate this requirement, however: small rotations of the camera can cause large object translations. Small or (as with SSDA) textured templates result in matching error surfaces with narrow extrema and thus constrain the range of interframe translation that can be successfully tracked. Another drawback of gradient descent techniques is that the search is inherently serial, whereas NCC permits parallel implementation.

Each of the approaches discussed above has been advocated by various authors, but there are fewer comparisons between approaches. Reference [19] derives an optimal feature tracking scheme within the gradient search framework, but the limitations of this framework are not addressed. An empirical study of five template matching algorithms in the presence of various image distortions [4] found that NCC provides the best performance in all image categories, although one of the cheaper algorithms performs nearly as well for some types of distortion. A general hierarchical framework for motion tracking is discussed in [1]. A correlation based matching approach is selected though gradient approaches are also considered.

Despite the age of the NCC algorithm and the existence of more recent techniques that address its various shortcomings, it is probably fair to say that a suitable replacement has not been universally recognized. NCC makes few requirements on the image sequence and has no parameters to be searched by the user. NCC can be used `as is' to provide simple feature tracking, or it can be used as a component of a more sophisticated (possibly multi-resolution) matching scheme that may address scale and rotation invariance, feature updating, and other issues. The choice of the correlation coefficient over alternative matching criteria such as the sum of absolute differences has also been justified as maximum-likelihood estimation [18]. We acknowledge NCC as a default choice in many applications where feature tracking is not in itself a subject of study, as well as an occasional building block in vision and pattern recognition research (e.g. [3]). A fast algorithm is therefore of interest.

Consider the numerator in (2) and assume that we have images and in which the mean value has already been removed:

For a search window of size and a feature of size (3) requires approximately additions and multiplications.

Eq. (3) is a convolution of the image with the reversed feature

where is the Fourier transform. The complex conjugate accomplishes reversal of the feature via the Fourier transform property .

Implementations of the FFT algorithm generally require that

When

There are several well known ``fast'' convolution algorithms that do not use transform domain computation [13]. These approaches fall into two categories: algorithms that trade multiplications for additional additions, and approaches that find a lower point on the characteristic of (one-dimensional) convolution by embedding sections of a one-dimensional convolution into separate dimensions of a smaller multidimensional convolution. While faster than direct convolution these algorithms are nevertheless slower than transform domain convolution at moderate sizes [13] and in any case they do not address computation of the denominator of (2).

Because (4) can be efficiently computed in the transform domain, several transform domain methods of approximating the image energy normalization in (2) have been developed. Variation in the image energy under the template can be reduced by high-pass filtering the image before cross-correlation. This filtering can be conveniently added to the frequency domain processing, but selection of the cutoff frequency is problematic--a low cutoff may leave significant image energy variations, whereas a high cutoff may remove information useful to the match.

A more robust approach is

Examining again the numerator of (2), we note that the mean of the feature can be precomputed, leaving

Since

Examining the denominator of (2), the length of the feature vector can be precomputed in approximately operations (small compared to the cost of the cross-correlation), and in fact the feature can be pre-normalized to length one.

The problematic quantities are those in the expression . The image mean and local (RMS) energy must be computed at each

These quantities can be efficiently computed from tables containing the integral (running sum) of the image and image square over the search area, i.e.,

with when either

and similarly for the image sum under the feature.

The problematic quantity can now be computed with very few operations since it expands into an expression involving only the image sum and sum squared under the feature. The construction of the tables requires approximately operations, which is less than the cost of computing the numerator by (4) and considerably less than the required to compute at each

This technique of computing a definite sum from a precomputed running sum has been independently used in a number of fields; a computer graphics application is developed in [5]. If the search for the maximum of the correlation surface is done in a systematic row-scan order it is possible to combine the table construction and reference through state variables and so avoid explicitly storing the table. When implemented on a general purpose computer the size of the table is not a major consideration, however, and flexibility in searching the correlation surface can be advantageous. Note that the

The performance of this algorithm will be discussed in the context of special effects image processing. The integration of synthetic and processed images into special effects sequences often requires accurate tracking of sequence movement and features. The use of automated feature tracking in special effects was pioneered in movies such as

The algorithm described in this paper was developed for the movie

The relative performance of our algorithm is a function of both the search window size and the ratio of the feature size to search window size. Relative performance increases along the window size axis (Fig. 1); a higher resolution plot would show an additional ripple reflecting the relation between the search window size and the bounding power of two. The property that the relative performance is greater on larger problems is desirable. Table 1 illustrates the performance obtained in a special effects feature tracking application. Table 2 compares the performance of our algorithm with that of a high-end commercial image compositing package.

Note that while a small (e.g. ) feature size would suffice in an ideal digital image, in practice much larger feature sizes and search windows are sometimes required or preferred:

- The image sequences used in film and video are sometimes obtained from moving cameras and may have considerable translation between frames due to camera shake. Due to the high resolution required to represent digital film, even a small movement across frames may correspond to a distance of many pixels.
- The selected features are of course constrained to the available features in the image; distinct ``features'' are not always available at preferred scales and locations.
- Many potential features in a typical digitized image are either out of focus or blurred due to motion of the camera or object (Fig. 2). Feature match is also hindered by imaging noise such as film grain. Large features are more accurate in the presence of blur and noise.

The fast algorithm in some cases reduces high-resolution feature tracking from an overnight to an over-lunch procedure. With lower (``proxy'') resolution and faster machines, semi-automated feature tracking is tolerable in an interactive system. Certain applications in other fields may also benefit from the algorithm described here.

**1**-
P. Anandan, ``A Computational Framework and an Algorithm for the Measurement
of Visual Motion'',
*Int. J. Computer Vision*, 2(3), p. 283-310, 1989.

**2**-
D. I. Barnea, H. F. Silverman, ``A class of algorithms for fast
digital image registration'',
*IEEE Trans. Computers*, 21, pp. 179-186, 1972.

**3**-
R. Brunelli and T. Poggio, ``Face Recognition: Features versus Templates",
*IEEE Trans. Pattern Analysis and Machine Intelligence*, vol. 15, no. 10, pp. 1042-1052, 1993.

**4**-
P. J. Burt, C. Yen, X. Xu, ``Local Correlation Measures for
Motion Analysis: a Comparitive Study'',
*IEEE Conf. Pattern Recognition Image Processing*1982, pp. 269-274.

**5**-
F. Crow, ``Summed-Area Tables for Texture Mapping'',
*Computer Graphics*, vol 18, No. 3, pp. 207-212, 1984.

**6**-
R. O. Duda and P. E. Hart,
*Pattern Classification and Scene Analysis*, New York: Wiley, 1973.

**7**-
R. C. Gonzalez and R. E. Woods,
*Digital Image Processing*(third edition), Reading, Massachusetts: Addison-Wesley, 1992.

**8**-
A. Goshtasby, S. H. Gage, and J. F. Bartholic, ``A Two-Stage
Cross-Correlation Approach to Template Matching'',
*IEEE Trans. Pattern Analysis and Machine Intelligence*, vol. 6, no. 3, pp. 374-378, 1984.

**9**-
C. Kuglin and D. Hines, ``The Phase Correlation Image Alignment Method,''
Proc. Int. Conf. Cybernetics and Society, 1975, pp. 163-165.

**10**-
J. P. Lewis, ``Fast Template Matching'',
*Vision Interface,*p. 120-123, 1995.

**11**-
A. R. Lindsey, ``The Non-Existence of a Wavelet Function Admitting a Wavelet
Transform Convolution Theorem of the Fourier Type'', Rome Laboratory Technical Report
C3BB, 1995.

**12**-
B. D. Lucas and T. Kanade, ``An Iterative Image Registration Technique
with an Application to Stereo Vision'',
*IJCAI*1981.

**13**-
S. K. Mitra and J. F. Kaiser,
*Handbook for Digital Signal Processing*, New York: Wiley, 1993.

**14**-
A. V. Oppenheim and R. W. Schafer,
*Digital Signal Processing*, Englewood Cliffs, New Jersey: Prentice-Hall, 1975.

**15**-
D. Polk, ``Product Probe'' - Panasonic PV-IQ604,
*Videomaker*, October 1994, pp. 55-57.

**16**-
W. Pratt,
*Digital Image Processing*, John Wiley, New York, 1978.

**17**-
Qi Tian and M. N. Huhns, ``Algorithms for Subpixel Registration'',
*CVGIP*35, p. 220-233, 1986.

**18**-
T. W. Ryan, ``The Prediction of Cross-Correlation Accuracy in Digital
Stereo-Pair Images'', PhD thesis, University of Arizona, 1981.

**19**-
J. Shi and C. Tomasi, ``Good Features to Track'', Proc.
*IEEE Conf. on Computer Vision and Pattern Recognition*, 1994.

**20**-
*Flame*effects compositing software, Discreet Logic, Montreal, Quebec.

**21**-
*Advance*effects compositing software, Avid Technology, Inc., Tewksbury, Massachusetts.

**22**-
*After Effects*effects compositing software, Adobe (COSA), Mountain View, California.