Forward FFT an image and backward FFT an image to get the same result

Question

I am trying to FFT an image using the library from http://www.fftw.org/ so that I can do a convolution in the frequency domain. But I can't figure out how to make it work. To understand how to do this I am trying to forward FFT an image as an array of pixelcolors and then backward FFT it to get the same array of pixelcolors. Here's what I do:

fftw_plan planR, planG, planB;
fftw_complex *inR, *inG, *inB, *outR, *outG, *outB, *resultR, *resultG, *resultB;

//Allocate arrays.
inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        inR[y * width + x][0] = pixelColors[currentIndex];
        inG[y * width + x][0] = pixelColors[currentIndex + 1];
        inB[y * width + x][0] = pixelColors[currentIndex + 2];
    }
}

//Forward plans.
planR = fftw_plan_dft_2d(width, width, inR, outR, FFTW_FORWARD, FFTW_MEASURE);
planG = fftw_plan_dft_2d(width, width, inG, outG, FFTW_FORWARD, FFTW_MEASURE);
planB = fftw_plan_dft_2d(width, width, inB, outB, FFTW_FORWARD, FFTW_MEASURE);

//Forward FFT.
fftw_execute(planR);
fftw_execute(planG);
fftw_execute(planB);

//Backward plans.
planR = fftw_plan_dft_2d(width, width, outR, resultR, FFTW_BACKWARD, FFTW_MEASURE);
planG = fftw_plan_dft_2d(width, width, outG, resultG, FFTW_BACKWARD, FFTW_MEASURE);
planB = fftw_plan_dft_2d(width, width, outB, resultB, FFTW_BACKWARD, FFTW_MEASURE);

//Backward fft
fftw_execute(planR);
fftw_execute(planG);
fftw_execute(planB);

//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        pixelColors[currentIndex] = resultR[y * width + x][0];
        pixelColors[currentIndex + 1] = resultG[y * width + x][0];
        pixelColors[currentIndex + 2] = resultB[y * width + x][0];
    }
}

Could someone please show me an example of how to forward FFT an image and then backward FFT the image using FFTW to get the same result? I have been looking at a lot of examples showing how to use FFTW to FFT, but I can't figure out how it applies to my situation where I have an array of pixelcolors representing an Image.

Answer 1

One important thing to note when you do forward FFT followed by inverse FFT is that this normally results in a scaling factor of N being applied to the final result, i.e. the resulting image pixel values will need to be divided by N in order to match the original pixel values. (N being the size of the FFT.) So your output loop should probably look something like this:

//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        pixelColors[currentIndex] = resultR[y * width + x][0] / (width * height);
        pixelColors[currentIndex + 1] = resultG[y * width + x][0] / (width * height);
        pixelColors[currentIndex + 2] = resultB[y * width + x][0] / (width * height);
    }
}

Also note that you probably want to do a real-to-complex FFT followed by a complex-to-real IFFT (somewhat more efficient in terms of both memory and performance). For now though it looks like you're doing complex-to-complex in both directions, which is fine, but you're not filling your input arrays correctly. If you're going to stick with complex-to-complex then you probably want to change your input loop to something like this:

//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        inR[y * width + x][0] = (double)pixelColors[currentIndex];
        inR[y * width + x][1] = 0.0;
        inG[y * width + x][0] = (double)pixelColors[currentIndex + 1];
        inG[y * width + x][1] = 0.0;
        inB[y * width + x][0] = (double)pixelColors[currentIndex + 2];
        inB[y * width + x][1] = 0.0;
    }
}

i.e. the pixel values go into the real parts of the complex input values and the imaginary parts need to be zeroed.

One more thing to note: when you eventually get this working you'll find that performance is terrible - it takes a long time to create a plan relative to the time taken for the actual FFT. The idea is that you create the plan just once, but use it to perform many FFTs. So you'll want to separate out the plan creation from the actual FFT code and put it in an initialisation routine or constructor or whatever.

Answer 2

But if you use the realToComplex or the ComplexToRealFunction pay attention to the fact that the image will be stored in a matrix of dimensions [height x (width/2 +1)] and if you want to do some intermediate calculations in the frequency domain, they will become a bit harder...

Answer 3

The reason it didn't work is that fftw_plan_dft_2d() does some bench-marking to find the best algorithm and changes input data in the process, so you have to fill the input data after fftw_plan_dft_2d(), not before it.