I will analyze the algorithmic complexity of the convolution algorithms described in my previous posts.

To make things simpler, let’s assume that the dimensions of the image are *>=* the dimensions of the convolution kernel, and that both are square, with dimensions *SxS* and *sxs*, respectively.

**1- Naive algorithm — O(n^4)**

*“w·h operations for each of the W·H image pixels”.*

i.e., *S·S·s·s* operations. This is *quadratic with a heavy constant* for tiny kernels, but quickly becomes *quartic* for medium-to-large kernels.

The auxiliary buffer can be identical in size and bit-depth to the original image. So the memory usage factor is 2x.

**2- Separable convolution — O(n^3)**

*“One 1D convolution for each row + One 1D convolution for each column”.*

i.e., *2·S·S·s* operations. This is *quadratic with a bad constant* for small kernels, but becomes *cubic* for large kernels.

Remarkably, *the total running time depends on the dimensions of the image -and- the dimensions of the kernel*.

Again, the auxiliary buffer can be identical in size and bit-depth to the original image. So the memory usage factor is 2x.

**3- Convolution theorem — O(n^2·log(n))**

*“Two FFTs + One point-wise product + One IFFT”.*

Let’s call *S’* to the closest power-of-2 such that *S’>=S*. Then a proper implementation of the FFT/IFFT does (approx.) *2·S’·S’·log(S`)* operations, while the point-wise product does *S’·S’* operations. This makes the algorithm *O(S’·S’·log(S’))* with some heavy (quadratic) overhead due to the memory copying, padding, and the point-wise product.

Remarkably, *the total running time is independent of the size of the kernel*.

This algorithm is quite memory hungry, though, because two *S’xS’* complex-number buffers are required. This means two floating-point numbers per entry (i.e., the real/imaginary coefficients). The algorithm starts by copying the image/kernel to the real part of the corresponding complex-number buffer, leaving the imaginary coefficient and the surface excess filled with zeroes. Then the FFTs/product/IFFT happen *in-place*.

So the auxiliary memory required is *4·S’·S’* floating-point numbers.

In the worst case where *S* is a power-of-2 plus 1, *S’* gets nearly twice as large as *S*. If the original image is 8-bit and we are using single-precision floating-point math for the FFT/IFFT, this means a memory usage factor of 64x. In the case of an HDR (single-precision floating-point) grayscale image, we are speaking of a worst case scenario of 16x. In average, however, the memory usage factor is around 8x. If *S* is a power-of-2, then the memory usage factor goes down to 4x.

This image with heavy glare has been output with some Arion-related experimental tonemapping code I am working on these days.

**Conclusions:**

Assuming that we are only interested in sheer performance:

1- The FFT-based method is (by far) the fastest for large images/kernels. Interestingly, the algorithm is not affected by the size of the kernel, which can be as large as the (padded) image itself without a penalty.

2- The FFT-based method becomes even faster if the same kernel is applied multiple times. The kernel FFT can be calculated just once, and then be re-used.

3- Due to the heavy setup overhead in the FFT-based method, the separable method can be faster for small (separable) kernels where *s* is in the range of *log(S’)*.

Last, but certainly not least, **there is a much faster and more light-weight algorithm for the special case of Box/Gaussian Blur**. I will talk about this in a separate blog entry.