Forward FFT image and inverse FFT image to get the same result

Question

I am trying to perform an FFT image using the library http://www.fftw.org/ so that I can perform convolution in the frequency domain. But I can’t figure out how to make it work. To understand how to do this, I try to forward the FFT image as an array of pixelcolors, and then back to the FFT to get the same array of pixels. That's what I'm doing:

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]; } } 

Can someone please show me an example of how to forward an FFT image and then back an FFT image using FFTW to get the same result? I looked at a lot of examples showing how to use FFTW for FFT, but I can’t understand how this relates to my situation when I have an array of pixels representing the image.

Answer 1

It is important to note that when you perform an FFT transfer followed by a reverse FFT, this usually results in the scaling factor N being applied to the final result, i.e. the resulting pixel values ​​of the image must be divided by N in order to match the original pixel values. (N is the size of the FFT.) So your output loop should 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 will probably want to make a real complex FFT, followed by a complex and real IFFT (somewhat more efficient in terms of both memory and performance). At the moment, although it looks like you are doing a complex complex in both directions, and this is normal, but you are not populating your input arrays correctly. If you are going to stick with a complex complex, 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. pixel values ​​fall into the real parts of the complex input values, and the imaginary parts must be zeroed out.

Another thing to note: when you eventually get this job, you will find that the performance is terrible - it takes a lot of time to create a plan for the time it takes for the actual FFT. The idea is that you create a plan only once, but use it to execute many FFTs. Thus, you will want to separate the creation of the plan from the actual FFT code and place it in the initialization procedure or constructor or something else.

Answer 2

But if you use realToComplex or ComplexToRealFunction, note that the image will be stored in the [height x (width / 2 +1)] size matrix, and if you want to do some intermediate calculations in the frequency domain, they will become a little more complicated .. .