1
$\begingroup$

Sorry if this question is too simplistic, but I am new to image processing and after successfully writing two programs, a convolution program and a low pass filtering program, I noticed that for approximately the same output, the low pass filtering method takes 4 times the time!

Lets assume that the image is a square of side N. When convolving we can use an arbitrarily small kernel (out to $3\sigma$ for example) and put a Gaussian in it, say of size M=$3\sigma$. We then have to simply pad the image and the kernel into an N+M-1 side square array and run DFT onit, multiply it and run IDFT.

But for low pass filtering we have to create a padded square array of width 2*N. There is one less DFT required, but the one DFT and the IDFT on the array with 4 times the elements will take a significantly longer amount of time because the padded array is huge (look at Gonzalez&Woods Fig.4.36 as an example)! Since the DFT of a Gaussian kernel is also a Gaussian, then why should we use a Gaussian low pass filtering at all when it takes nearly 4 times the time (look below)?

Source code: I have put the source code in github.

I am adding some timings in my program to show what I mean. The image was 233*233 pixels. My CPU is clocked at 3.07GHz.

For a low pass filter, note that the FFTs and IFFTs are done on an array of size 466*466:

Gaussian low pass filter timing:

Padding:                0.002564 (seconds)
FFT:                    0.018210 (seconds)
Applying filter:        0.006686 (seconds)
IFFT:                   0.017254 (seconds)
Time for final step:    0.000550 (seconds)

Total:                  0.047424 (seconds)

For a convolution of the same image with a kernel of size 9*9, note that the three FFTs and IFFTs are done on arrays of size 241*241:

Gaussian convolution timing:

Making the kernel:      0.000077 (seconds)
Pad both to same size:  0.001389 (seconds)
Padded image FFT:       0.002699 (seconds)
Padded kernel FFT:      0.001866 (seconds)
Multiplying:            0.000384 (seconds)
IFFT:                   0.002241 (seconds)
Crop from padded:       0.000574 (seconds)

Total:          0.011149 (seconds)
$\endgroup$
18
  • $\begingroup$ Have you tried using FFT instead of raw DFT? $\endgroup$ Aug 22, 2013 at 11:20
  • $\begingroup$ Yes, I am using FFTW in C. $\endgroup$
    – makhlaghi
    Aug 22, 2013 at 11:22
  • $\begingroup$ If for both methods you use same code structures then comparing complexity of two methods would be the key answer to your question. $\endgroup$ Aug 22, 2013 at 11:31
  • $\begingroup$ I added some information about the timings. The logic is very simple: For a LPF you need FFT to work on N*N pixels. While for convolution you just need to work on two N+M-1 arrays. It is clear that LPF will be much more time consuming as long is M is small. $\endgroup$
    – makhlaghi
    Aug 22, 2013 at 11:34
  • 1
    $\begingroup$ @makhlaghi Check out this link (and the intro) for convolution (in the spatial domain). ftp.isy.liu.se/pub/colour/misc/convolution.txt Conceptually, convolution is filtering. In your case both ways could be called 'low pass filtering'. In both, you are filtering via multiplication the Fourier domain so its not technically convolution. Convolution is kind of a shift-multiply-add of the spatial kernel and the image. Using your LPF method (the second one) if your kernel is well approximated when restricted to M x M pixels, then padding the image to N+M-1, not 2N, will be fine. $\endgroup$ Aug 23, 2013 at 11:16

1 Answer 1

4
$\begingroup$

In practice gaussian low-pass filters are written using F.I.R. or I.I.R. filters, first the filter runs left then right on each horizontal line, then the filter runs up then down on each vertical line. It's much cheaper than doing FFTs and IFFTs.

Your two operations are not performing the same operation. The effects of the filtering around the edges will be very different. You are doing a circular convolution. For the Discrete Fourier Transform the convolution theorem says that the multiplication of fourier transforms is equivalent to convolving the image infinitely repeated/tiled in the plane. So around the left edge, for example you are blurring together some components from the right side of the picture. (For example try performing your operations on a picture that is black on the left half and white on the right half.)

Your "low pass" filter is padding the image and the kernel to try to reduce some of the edge effects. It's more expensive, but it is closer to what you would get if you did a "real" convolution of a gaussian with your image.

In your application the circular convolution effect of blurring in data from the other side of the image may be okay, but in many application it is unacceptable.

As to the reason that the FFT of an $N\times N$ image is 10 times faster than the FFT of a $2N \times 2N$ image. A 2d FFT is done by taking 1d FFTs of each horizontal line, then 1d FFTs of each vertical line. So in the $N \times N$ case you are doing $2N$ 1-d FFTs, each of which take $O(N \log N)$ for a total of $O(2N^2 \log N)$. In the $2 N \times 2 N$ case you are doing $4N$ 1-d FFTs each of which take $O(2N \log 2N)$ for a total of (a little more than) $O(8N^2 \log N)$. Add in cache effects for the larger image and 10x slower seems quite reasonable.

$\endgroup$
9
  • $\begingroup$ Thanks. I was not familiar with IIR filters! Do you know of its implementation in C? Or anywhere explaining the details of its algorithm so I can write it my self? $\endgroup$
    – makhlaghi
    Aug 22, 2013 at 12:08
  • 1
    $\begingroup$ The code for an I.I.R. filter is trivial. The hard parts are (a) choosing the coefficients and (b) dealing with the transients at the edges. In general Matlab is used to help design the filter. In the case of gaussians the filters are based on the following paper: citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.12.2826. A nice tutorial is cwp.mines.edu/Meetings/Project06/cwp546.pdf‎. Ways of dealing with transients are described in hal.inria.fr/inria-00548616/en. $\endgroup$ Aug 22, 2013 at 12:21
  • $\begingroup$ I understand the edge problem you explained and saw it in my programs. I am reading about IIR filters now. Based on your explanations (..."but in many application it is unacceptable."), I understand that the edge-effect is removed with these implementations, am I correct? $\endgroup$
    – makhlaghi
    Aug 22, 2013 at 12:32
  • $\begingroup$ That last paper (the one from hal.inria.fr) is specifically about dealing with the edge effects. It's actually a bit tricky to do it efficiently. $\endgroup$ Aug 22, 2013 at 12:37
  • 1
    $\begingroup$ "In practice gaussian low-pass filters are written using I.I.R. filters". I have never heard of this before, can you provide a link about this? I thought gaussian filters are always FIRs. $\endgroup$ Aug 22, 2013 at 14:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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