Search code examples
c++wolfram-mathematicagaussianconvolutionfftw

Blur a matrix using Fast Fourier Transforms


I want to blur values in matrix, so that in neighboring elements we'll not have sharp transitions.

From the Wikipedia page Gaussian Blur I've found some info on Gaussian blurring. I've tried it with the most simple algorithm, and, hence the run-time was too long. Frankly, I'm not sure if my implementation is correct, as on boundary tiles sharp transition still exists.

I've noticed that this blurring could be done with discrete Fourier transforms which is much faster, but I couldn't figure it out.

So, the idea is that we can get blurred matrix with the formulas below:

blurredMatrix = IFFT[FFT[initialMatrix]FFT[weightingFunction]]

Where FFT/IFFT are Fast Fourier Transform/Inverse Fast Fourier Transform.

Currently I'm trying to do some testing on Wolfram Mathematica to make sure that this kind of approximation with Fourier transforms is correct.

I'm using GaussianMatrix as weightingFunction.

I need 2d blurring, so I've created Gaussian matrix as below:

Suppose our initial matrix has nxn sizes, where n = 2k+1

G = Chop[GaussianMatrix[k] GaussianMatrix[k], 10^6]

And, then, I've tried to create blurredMatrix as below:

blurredMatrix = Chop[FourierDCT[(FourierDCT[G]) (FourierDCT[initialMatrix]), 3], 10^-6]

But I'm getting zeros as a result.

Seems I'm doing it all wrong.

Also, I've tried another approach:

f[xi_, yj_] := 1/(2 \[Pi] \[Sigma]^2) Exp[-(((xi^2) + (yj^2) )/(2 \[Sigma]^2))];<br/>
[Sigma] = 3;<br/>
G = Chop[N[Table[f[i, j], {i, 1, 100}, {j, 1, 100}]]]; <br/>
Tavg = Chop[ 1000 InverseFourier[(Fourier[G]) (Fourier[T]) ], 10^-6]; <br/>

With this approach the picture looks fine (the image is blurred), but there is big difference between between values of blurredMatrix and initialMatrix.

Seems there are some normalization or other problems.

I need to write the code in C/C++, there is a library FFTW library in C which supports discrete Fourier transforms.

Please let me know if this is a wrong way of blurring and there are other possibilities to do what I want.


Solution

  • Using FFT to do convolutions is only efficient when you have very large convolution kernels. In most blurring applications the kernel is much much smaller than the image, e.g. 3x3, so FFT would be significantly slower.

    There are many implementations for doing small-kernel convolutions.
    Most modern hardware supports such intrinsic operations (MMX, SSE, GPUs...).
    FFT is probably not the way to go in your case.

    In C++, OpenCV supports both cross-platform and hardware accelerated image convolutions. Convolutions really are one of (if not) the most foundational operation of any image and signal processing package.