The discussion revolves around finding a suitable free library for performing two-dimensional Fourier Transforms (FFT or DFT) in C or C++. Participants share their experiences with existing libraries, particularly FFTW, and discuss code-related issues encountered during implementation.
Participants express varying opinions on the use of FFTW, with some advocating for its efficiency while others report issues. There is no consensus on the best approach or library, and the discussion remains unresolved regarding the specific problems encountered.
Participants mention potential issues with memory allocation and the order of operations in the code, but these remain unverified. The discussion reflects a range of experiences with FFTW and alternative methods without reaching a definitive solution.
...
#include <complex>//Includes C++ complex numbers (uses templates)
#include <fftw3.h> // Library to perform FFT transform (GPL license)
const int BLOCKSIZE= 10;
...
// Applies an MTF to an image. This filters the image. It is done in the frequency domain
// image_in[BLOCKSIZE][BLOCKSIZE] : Input image
// double MTF[BLOCKSIZE][BLOCKSIZE] : Modulation Transfer function (MTF)
// double image_out[BLOCKSIZE][BLOCKSIZE]: Output image after applying the MTF
//void applyMTF(const double image_in[][BLOCKSIZE], const double MTF[][BLOCKSIZE], double image_out[][BLOCKSIZE]){
void applyMTF(){
fftw_complex *image_v_spat; //Image vector in spatial domain
fftw_complex *image_v_freq; //Image vector in frequency domain
fftw_complex *MTF_v; //Modulation Transfer Function (in frequency domain by definition)
fftw_plan p; //Plan used by the fftw library
int i= 0, j= 0, counter= 0; //Loop counters
//Input parameters temporarily defined inside the function
double image_in[BLOCKSIZE][BLOCKSIZE];
double image_out[BLOCKSIZE][BLOCKSIZE];
double MTF[BLOCKSIZE][BLOCKSIZE];
//Fill image and MTF with dummy values:
for(i= 0; i< BLOCKSIZE; i++){
for(j= 0; j< BLOCKSIZE; j++){
image_in[i][j]= (double)(i*i - i + j);
MTF[i][j]= (double)(j);
}
}
//Allocate memory
image_v_spat = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCKSIZE * BLOCKSIZE);
image_v_freq = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCKSIZE * BLOCKSIZE);
MTF_v = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCKSIZE * BLOCKSIZE);
//Expand the image and MTF arrays (2-dimensional real) into vectors (1-dimensional complex array). Fill the imaginary part with zeros.
counter= 0;
for(i= 0; i< BLOCKSIZE; i++){
for(j= 0; j< BLOCKSIZE; j++){
counter++;
//Image
image_v_spat[counter][0]= image_in[i][j]; //Real part
image_v_spat[counter][1]= 0; //Imaginary part
//Same for MTF
MTF_v[counter][0]= MTF[i][j]; //Real part
MTF_v[counter][1]= 0; //Imaginary part
}
}
//Create the plan for the DFT
p= fftw_plan_dft_2d(BLOCKSIZE, BLOCKSIZE, image_v_spat, image_v_freq, FFTW_FORWARD, FFTW_ESTIMATE);
//Execute the plan (perform the Fourier transform of the image)
fftw_execute(p);
//Multiply the image in frequency domain times the MTF taking into account that they are complex numbers: (a0 + a1 i) * (b0 + b1 i) = (a0*b0 - a1*b1) + (a0*b1 + a1*b0)i
//The real part is stored in index 0 and the complex part in index 1
for(counter=0; counter < BLOCKSIZE*BLOCKSIZE; counter++){
image_v_freq[counter][0]= image_v_freq[counter][0] * MTF[counter][0] - image_v_freq[counter][1] * MTF[counter][1]; //Real part
image_v_freq[counter][1]= image_v_freq[counter][0] * MTF[counter][1] + image_v_freq[counter][1] * MTF[counter][0]; //Imaginary part
}
//Do the inverse Fourier transform
//First create the plan for the DFT
p= fftw_plan_dft_2d(BLOCKSIZE, BLOCKSIZE, image_v_freq, image_v_spat, FFTW_BACKWARD, FFTW_ESTIMATE);
//Then execute the plan (perform the inverse Fourier transform of the image)
fftw_execute(p);
//Do the absolute value and store the resulting image in a real 2D matrix instead of a complex 1D vector. This image is then returned by the function by ref.
counter= 0;
for(i= 0; i< BLOCKSIZE; i++){
for(j= 0; j< BLOCKSIZE; j++){
counter++;
image_out[i][j]= sqrt(image_v_spat[counter][0]*image_v_spat[counter][0] + image_v_spat[counter][1]*image_v_spat[counter][1]);
}
}
//Deallocate memory
fftw_destroy_plan(p);
fftw_free(image_v_spat);
fftw_free(image_v_freq);
fftw_free(MTF_v);
return;
}#include <iostream>
#include <cstdlib>
#include <fftw3.h>
#include <fstream>
using namespace std;
int main()
{
const int Lx = 4; // Array dimensions
const int Lt = 2;
int var_x;
int var_T;
double F[Lx][Lt][2];
fftw_complex *in, *out;
fftw_plan p;
// Declare one-dimensional contiguous arrays of dimension Lx*Lt
in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Lx*Lt);
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Lx*Lt);
// create plan, an object containing necessary parameters
p = fftw_plan_dft_2d(Lx, Lt, in, out, FFTW_FORWARD, FFTW_MEASURE);
ifstream input;
input.open("2d_function");
// Put input data into F[x][T]
// Careful - the loop order depends on how the data is recorded
for(int x = 0; x < Lx; x++)
{
for(int T = 0; T < Lt; T++)
{
input >> var_T;
input >> var_x;
input >> F[x][T][0];
//input >> Fi[x][T][1]; Can get imaginary part if it exists
}
} // Now turn the array F[x][T] into C-order (row major) 1-d array
// It also happens to be an array of fftw_complex values
// However, we only need to store each value in the real part
// The order here is important. Rightmost index loops on the inside
for(int x = 0; x < Lx; x++)
{
for(int T = 0; T < Lt; T++)
{
in[T + Lt*x][0] = F[x][T][0];
// in[T + Lt*x][1] = F[x][T][1] : Can enter imaginary part if needed
}
}
input.close(); // Perform FFT
fftw_execute(p); // Get output data into 2-d form - use same array as for input
for(int x = 0; x < Lx; x++)
{
for(int T = 0; T < Lt; T++)
{
F[x][T][0] = out[T + Lt*x][0];
F[x][T][1] = out[T + Lt*x][1]; // Can get imaginary data too if it exists
}
}
// Print output to file
ofstream output;
output.open("2d_transform_out");
for(int x = 0; x < Lx; x++)
{
for(int T = 0; T < Lt; T++)
{
output << x << "\t" << T << "\t" << F[x][T][0] << "\t" << "+ i " << F[x][T][1] << "\n";
}
output << "\n";
}
fftw_destroy_plan(p);
fftw_free(in);
fftw_free(out);
output.close();
return 0;
}