Speedup

Readers may recognize Eq. 5 (or Step 1 of the modified algorithm in Section 2.2) as the two-dimensional, digital counterpart of the convolution operation. A well-known property of the convolution of two functions, $f(t)$ and $g(t)$, is that it is the inverse Fourier transform of the product of the Fourier transforms of the two signals (For an explanation of this theorem, consult any introductory engineering mathematics text, such as (Kreyszig, 1983)), i.e.

$$(f * g)(t) = \mathcal{F}^{-1} ( \mathcal{F}(f(t)) . \mathcal{F}(g(t)) )$$ (7)

where $\mathcal{F}( f(t) )$ is the Fourier transform of $f(t)$. Denoting the Digital Fourier Transform of an image X by X, we can rewrite Eq. 5 as:

$${\bf F_k}_{r,s} = {\bf E_k}_{r,s} . {\bf D}_{r,s}$$ (8)

where the new indices represent the idea that we are operating on discrete spatial frequencies. Comparing this with Eq. 6, we see that we can obtain the filtered image, $F$ as follows:

$$F = \bigvee_{k,\forall i,j} \{ \mathcal{F}^{-1} ( \frac{1}{N_k}{\bf E_k}_{r,s} . {\bf D}_{r,s} ) \}$$ (9)

Note that the indices are back to $i,j$, since the inverse transform puts us back into spatial coordinates. The two Fourier Transforms (${\bf E_k}$ and ${\bf D}$) must be the same size if we are to do an element-by-element multiplication. So, one will have to pad the elliptical filters with zeros until they are the same size as the image.

Eq. 9 is the entire filtering process. Described in words:

1. Do once: Compute the elliptical filters (1's in the cross-hatched areas in Figure 1 and 0 everywhere else) and compute the $N_k$, the number of 1's in the orientation. Divide the filters by $N_k$'s. Then, pad the filters with zeros until they are of size $N\times N$. Compute the DFTs, ${\bf E_k}$'s of the elliptical filters.
2. Compute the DFT, ${\bf D}$, of the $N\times N$ image, $D$, to be filtered.
3. Do $q$ times: For each orientation of the filter, obtain the element-by-element product of the DFTs: ${\bf E_k . D}$. Then, take the inverse DFT of the product.
4. Assign to each point in the filtered image the maxima of the results from each orientation's image at that point.

Fast algorithms exist to compute the DFT for multiples of small prime numbers. ⁶Even if the original image is not of a convenient size, it can be padded with zeros to make it of such a size. These algorithms typically take of the order of NlogN operations to compute the DFT of a sequence of length N.

Since the images are of size $N\times N$, the computational overhead of step 2 in our technique is $(N^2)*log_2(N^2)$. In step 3, we need to do $N^2$ multiplications and take an inverse DFT, making the overhead $N^2 + N^2*log_2(N^2)$. Since step 3 has to be done $q$ times, the total computational overhead of the modified technique is only $(3q+2)*N^2*log_2(N)$. As in Section 2.3, we don't count the maximum operation in our computational overhead calculation. Note that our computational overhead does not depend on $2p+1$, the size of the filter. This is because we pad our filters to match the size of the images.

Again using a $512\times512$ image and a $15\times64$ elliptical filter in 10-degree increments, we have $N=512$, $p=32$ and $q=18$. The computational overhead of the modified technique is of the order of 132 million operations. Recall that the original technique took on the order of 19 billion operations. We get an improvement of more than two orders of magnitude. In practice, we won't achieve such a high improvement because the original technique can be implemented with integer arithmetic but the DFT technique requires floating point operations, which are slower on most computers.

Comparing the orders of magnitude calculations, we expect an improvement of $1.33*p^2 / log_2(N)$. ⁷ For a $512\times512$ image with a filter size of $15\times64$, the improvement should, in theory, be about 150 times. In practice (as can be seen from Table 1), we can do the filtering about 50 times faster by using the technique described in this paper. It is also worth noting that for small filter sizes ($p$ less than 5), the technique actually degrades performance.

The idea behind the speedup in this paper, that convolution can be performed faster on computers using Fourier Transforms, is more than 30 years old; Cooley and Tukey (1965) introduced the fast algorithm that forms the basis for most FFT algorithms today (including the one we used, by Frigo and Johnson (1998)). By 1967, there was a special issue of the IEEE Transactions on Audio and Electroacoustics devoted to Fast Fourier Transform methods.