Fast 2D convolution for DSP

Question

I want to implement some image-processing algorithms which are intended to run on a beagleboard. These algorithms use convolutions extensively. I'm trying to find a good C implementation for 2D convolution (probably using the Fast Fourier Transform). I also want the algorithm to be able to run on the beagleboard's DSP, because I've heard that the DSP is optimized for these kinds of operations (with its multiply-accumulate instruction).

I have no background in the field so I think it won't be a good idea to implement the convolution myself (I probably won't do it as good as someone who understands all the math behind it). I believe a good C convolution implementation for DSP exists somewhere but I wasn't able find it?

Could someone help?

EDIT: Turns out the kernel is pretty small. Its dimensions are either 2X2 or 3X3. So I guess I'm not looking for an FFT-based implementation. I was searching for convolution on the web to see its definition so I can implement it in a straight forward way (I don't really know what convolution is). All I've found is something with multiplied integrals and I have no idea how to do it with matrices. Could somebody give me a piece of code (or pseudo code) for the 2X2 kernel case?

Answer 1

What are the dimensions of the image and the kernel ? If the kernel is large then you can use FFT-based convolution, otherwise for small kernels just use direct convolution.

The DSP might not be the best way to do this though - just because it has a MAC instruction doesn't mean that it will be more efficient. Does the ARM CPU on the Beagle Board have NEON SIMD ? If so then that might be the way to go (and more fun too).

For a small kernel, you can do direct convolution like this:

// in, out are m x n images (integer data)
// K is the kernel size (KxK) - currently needs to be an odd number, e.g. 3
// coeffs[K][K] is a 2D array of integer coefficients
// scale is a scaling factor to normalise the filter gain

for (i = K / 2; i < m - K / 2; ++i) // iterate through image
{
  for (j = K / 2; j < n - K / 2; ++j)
  {
    int sum = 0; // sum will be the sum of input data * coeff terms

    for (ii = - K / 2; ii <= K / 2; ++ii) // iterate over kernel
    {
      for (jj = - K / 2; jj <= K / 2; ++jj)
      {
        int data = in[i + ii][j +jj];
        int coeff = coeffs[ii + K / 2][jj + K / 2];

        sum += data * coeff;
      }
    }
    out[i][j] = sum / scale; // scale sum of convolution products and store in output
  }
}

You can modify this to support even values of K - it just takes a little care with the upper/lower limits on the two inner loops.

Answer 2

I know it might be off topic but due to the similarity between C and JavaScript I believe it could still be helpful. PS.: Inspired by @Paul R answer.

Two dimensions 2D convolution algorithm in JavaScript using arrays

function newArray(size){
    var result = new Array(size);
    for (var i = 0; i < size; i++) {
        result[i] = new Array(size);
    }

    return result;
}

function convolveArrays(filter, image){
    var result = newArray(image.length - filter.length + 1);

    for (var i = 0; i < image.length; i++) {
        var imageRow = image[i];
        for (var j = 0; j <= imageRow.length; j++) {

            var sum = 0;
            for (var w = 0; w < filter.length; w++) {
                if(image.length - i < filter.length) break;

                var filterRow = filter[w];
                for (var z = 0; z < filter.length; z++) {
                    if(imageRow.length - j < filterRow.length) break;
                    sum += image[w + i][z + j] * filter[w][z];
                }
            }

            if(i < result.length && j < result.length)
                result[i][j] = sum;
        }   
    }

    return result;
}

You can check the full blog post at http://ec2-54-232-84-48.sa-east-1.compute.amazonaws.com/two-dimensional-convolution-algorithm-with-arrays-in-javascript/