Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Fast way to implement 2D convolution in C

Tags:

c

image

filter

2d

I am trying to implement a vision algorithm, which includes a prefiltering stage with a 9x9 Laplacian-of-Gaussian filter. Can you point to a document which explains fast filter implementations briefly? I think I should make use of FFT for most efficient filtering.

like image 368
Atilla Filiz Avatar asked Mar 24 '09 09:03

Atilla Filiz


People also ask

What is the fastest way to take the convolution of an image?

FFT is the fastest technique known for convolving signals, and FFTW is the fastest free library available for computing the FFT.

How do you do convolution in 2D?

The 2D convolution is a fairly simple operation at heart: you start with a kernel, which is simply a small matrix of weights. This kernel “slides” over the 2D input data, performing an elementwise multiplication with the part of the input it is currently on, and then summing up the results into a single output pixel.

How many for loops are needed to implement a 2D convolution?

C++ Algorithm for Convolution 2D We need 4 nested loops for 2D convolution instead of 2 loops in 1D convolution. The above snippet code is simple and easiest way to understand how convolution works in 2D.

How do you find the convolution of two arrays?

Given two arrays A[] and B[] consisting of N and M integers respectively, the task is to construct a convolution array C[] of size (N + M – 1). The convolution of 2 arrays is defined as C[i + j] = ∑(a[i] * b[j]) for every i and j.


3 Answers

Are you sure you want to use FFT? That will be a whole-array transform, which will be expensive. If you've already decided on a 9x9 convolution filter, you don't need any FFT.

Generally, the cheapest way to do convolution in C is to set up a loop that moves a pointer over the array, summing the convolved values at each point and writing the data to a new array. This loop can then be parallelised using your favourite method (compiler vectorisation, MPI libraries, OpenMP, etc).

Regarding the boundaries:

  • If you assume the values to be 0 outside the boundaries, then add a 4 element border of 0 to your 2d array of points. This will avoid the need for `if` statements to handle the boundaries, which are expensive.
  • If your data wraps at the boundaries (ie it is periodic), then use a modulo or add a 4 element border which copies the opposite side of the grid (abcdefg -> fgabcdefgab for 2 points). **Note: this is what you are implicitly assuming with any kind of Fourier transform, including FFT**. If that is not the case, you would need to account for it before any FFT is done.

The 4 points are because the maximum boundary overlap of a 9x9 kernel is 4 points outside the main grid. Thus, n points of border needed for a 2n+1 x 2n+1 kernel.

If you need this convolution to be really fast, and/or your grid is large, consider partitioning it into smaller pieces that can be held in the processor's cache, and thus calculated far more quickly. This also goes for any GPU-offloading you might want to do (they are ideal for this type of floating-point calculation).

like image 158
Phil H Avatar answered Oct 21 '22 06:10

Phil H


Here is a theory link http://hebb.mit.edu/courses/9.29/2002/readings/c13-1.pdf

And here is a link to fftw, which is a pretty good FFT library that I've used in the past (check licenses to make sure it is suitable) http://www.fftw.org/

All you do is FFT your image and kernel (the 9x9 matrix). Multiply together, then back transform.

However, with a 9x9 matrix you may still be better doing it in real coordinates (just with a double loop over the image pixels and the matrix). Try both ways!

like image 36
Greg Reynolds Avatar answered Oct 21 '22 08:10

Greg Reynolds


Actually you don't need to use a FFT size large enough to hold the entire image. You can do a lot of smaller overlapping 2d ffts. You can search for "fast convolution" "overlap save" "overlap add".

However, for a 9x9 kernel. You may not see much advantage speedwise.

like image 2
Mark Borgerding Avatar answered Oct 21 '22 07:10

Mark Borgerding