Take the 2-minute tour ×
Programmers Stack Exchange is a question and answer site for professional programmers interested in conceptual questions about software development. It's 100% free, no registration required.

I am working with 2d correlation for image processing techniques (pattern recognition etc...). I was wondering if there is a theoretical approach on how to tell when to use multiplication in frequency space over correlation in time space. For sizes of 2x frequency space is obviously faster but how about small, prime sizes like e.g. 11?

share|improve this question

1 Answer 1

up vote 10 down vote accepted

I'll assume this is being done on a conventional CPU, one core, executing one simple thread, no fancy hardware. If there is more than that going on, it can probably be accounted for with adjustments to the reasoning for a simpler system. Not much more can be said without either a specific system to discuss, or a whole textbook or research paper to cover a range of possibilities.

I wouldn't worry about power-of-two sizes. It doesn't matter. FFT algorithms with the butterfly units and all that exist for factors of 3, or any small number, not just 2. There are clever algorithms for prime-sized data series, too. I don't like quoting wikipedia due to its impermanent nature, but anyway: "there are FFTs with O(N log N) complexity for all N, even for prime N" (http://en.wikipedia.org/wiki/Fast_Fourier_transform) Implementations of FFTs for arbitrary N can be found in the GPL'd library FFTW (http://fftw.org/)

The only trustworthy way in terms of serious engineering is to build and measure, but we certainly can get an idea from theory, to see relationships between variables. We need estimates of how many arithmetic operations are involved for each method.

Multiplying is still slower than addition on most CPUs, even if the difference has shrunken tremendously over the years, so let's just count multiplications. Accounting also for addition takes a bit more thinking and keeping track of stuff.

A straightforward convolution, actually multiplying and adding using the convolution kernel, repeating for each output pixel, needs (W^2)(K^2) multiplications, where W is the number of pixels along one side of the image (assuming square for simplicity), and K is the size of the convolution kernel, as pixels along one side. It takes K^2 multiplications to compute one output pixel using the kernel and same-size portion of the input image. Repeat for all output pixels, which number the same as in the input image.

N_mults_direct = (W^2)*(K^2)

To do the job in Fourier space, we must Fourier transform the image. This is done by applying an FFT to each column separately, and then to each row. The FFT for N data points takes about 2N*log(N) multiplications; we want N to be W, the length of one column or row. All logarithms here are base two. There are W rows and W columns, so after all the FFTs are done, we have done 2W*(2W*log(W)) multiplications. Double that, because after we multiply by the Fourier transform of the kernel, we have to inverse-Fourier the data to get back to sensible image. That's 8*W^2*log(W). Of course, multiplying by the Fourier transform of the kernel has to be done, another W^2 multiplications. (Done once, not once per output pixel, per row or anything.) These are complex multiplications, so that's 4*W^2 real multiplications.

So, unless I goofed up (and I probably did) we have

N_mults_fouriermethod = 4*W^2*(2*log(W)+1) 

When do we want to do things the direct way? When K is sufficiently small to make (W^2)*(K^2) smaller than 4*W^2*(2*log(W)+1). A common factor of W^2 is easily factored out. We can probably drop the "+1" since we're dealing with idealized estimates. The +1 is likely lost in errors relative to actual implementations, from not counting additions, loop overheads and so on. That leaves:

K^2 < 8*log(W)

This is the approximate condition for choosing a direct approach over a frequency space approach.

Note that correlation of two same-size images is just like convolving with a kernel of size K=W. Fourier space is always the way to do it.

This can be refined and argued over to account for overhead, pipelining of opcodes, float vs. fixed-point, and thrown out the window with GPGPU and specialized hardware.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.