Hello StackOverflow community, i have a problem with the dft algorithm of the fftw library. All i want to do is to transform a certain pattern forward and backward to receive the input pattern again, of course there will be some sort of filtering in between the transformations later on.
So, what my program does atm is:
My problem is that when is use no filtering my backward transformed image is wrapped somehow and i can't calculate the correct magnitude, which should be indentical to my input image / test signal. When i set the fitler/"window" to a value of 0.5 the backward transformation works fine, but my input image is just half as bright as it should be.
The following image illustrates my problem: (from top left to bottom right) 1. Input signal, 2. Real part of backward transformation, 3. From backward transformated data calculated magnitude, 4. Input signal multiplied with 0.5, 5. Real part of backward transformation, 6. From backward transformated data calculated magnitude. http://imageshack.com/a/img538/5426/nbL9YZ.png
Does anybody have an idea why the dft performs in that way?! It's kind of strange...
My code looks like this atm:
/***** parameters **************************************************************************/
int imSize = 256;
int imN = imSize * imSize;
char* interferogram = new char[imN];
double* spectrumReal = new double[imN];
double* spectrumImaginary = new double[imN];
double* outputReal = new double[imN];
double* outputImaginary = new double[imN];
double* phase = new double[imN];
char* spectrumRealChar = new char[imN];
char* spectrumImaginaryChar = new char[imN];
char* outputRealChar = new char[imN];
char* outputImaginaryChar = new char[imN];
char* phaseChar = new char[imN];
Mat interferogramMat = Mat(imSize, imSize, CV_8U, interferogram);
Mat spectrumRealCharMat = Mat(imSize, imSize, CV_8U, spectrumRealChar);
Mat spectrumImaginaryCharMat = Mat(imSize, imSize, CV_8U, spectrumImaginaryChar);
Mat outputRealCharMat = Mat(imSize, imSize, CV_8U, outputRealChar);
Mat outputImaginaryCharMat = Mat(imSize, imSize, CV_8U, outputImaginaryChar);
Mat phaseCharMat = Mat(imSize, imSize, CV_8U, phaseChar);
/***** compute interferogram ****************************************************************/
fill_n(interferogram, imN, 0);
double value = 0;
double window = 0;
for (int y = 0; y < imSize; y++)
{
for (int x = 0; x < imSize; x++)
{
value = 127.5 + 127.5 * cos((2*PI) / 10000 * (pow(double(x - imSize/2), 2) + pow(double(y - imSize/2), 2)));
window = 1;
value *= window;
interferogram[y * imSize + x] = (unsigned char)value;
}
}
/***** create fftw arays and plans **********************************************************/
fftw_complex* input;
fftw_complex* spectrum;
fftw_complex* output;
fftw_plan p_fw;
fftw_plan p_bw;
input = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
spectrum = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
output = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * imN);
p_fw = fftw_plan_dft_2d(imSize, imSize, input, spectrum, FFTW_FORWARD, FFTW_ESTIMATE);
p_bw = fftw_plan_dft_2d(imSize, imSize, spectrum, output, FFTW_BACKWARD, FFTW_ESTIMATE);
/***** copy data ****************************************************************************/
for (int i = 0; i < imN; i++)
{
input[i][0] = double(interferogram[i]) / 255.;
input[i][1] = 0.;
spectrum[i][0] = 0.;
spectrum[i][1] = 0.;
output[i][0] = 0.;
output[i][1] = 0.;
}
/***** FPS algorithm ************************************************************************/
fftw_execute(p_fw);
fftw_execute(p_bw);
for (int i = 0; i < imN; i++)
{
phase[i] = sqrt(pow(output[i][0], 2) + pow(output[i][1], 2));
}
/***** copy data ****************************************************************************/
for (int i = 0; i < imN; i++)
{
spectrumReal[i] = spectrum[i][0];
spectrumImaginary[i] = spectrum[i][1];
outputReal[i] = output[i][0] / imN;
outputImaginary[i] = output[i][1];
}
SaveCharImage(interferogram, imN, "01_interferogram_512px_8bit.raw");
SaveDoubleImage(spectrumReal, imN, "02_spectrum_real_512px_64bit.raw");
SaveDoubleImage(spectrumImaginary, imN, "03_spectrum_imaginary_512px_64bit.raw");
SaveDoubleImage(outputReal, imN, "03_output_real_512px_64bit.raw");
DoubleToCharArray(spectrumReal, spectrumRealChar, imSize);
DoubleToCharArray(spectrumImaginary, spectrumImaginaryChar, imSize);
DoubleToCharArray(outputReal, outputRealChar, imSize);
DoubleToCharArray(outputImaginary, outputImaginaryChar, imSize);
DoubleToCharArray(phase, phaseChar, imSize);
/***** show images **************************************************************************/
imshow("interferogram", interferogramMat);
imshow("spectrum real", spectrumRealCharMat);
imshow("spectrum imaginary", spectrumImaginaryCharMat);
imshow("out real", outputRealCharMat);
imshow("out imaginary", outputImaginaryCharMat);
imshow("phase", phaseCharMat);
int key = waitKey(0);
Here are some lines of your code :
char* interferogram = new char[imN];
...
double value = 0;
double window = 0;
for (int y = 0; y < imSize; y++)
{
for (int x = 0; x < imSize; x++)
{
value = 127.5 + 127.5 * cos((2*PI) / 10000 * (pow(double(x - imSize/2), 2) + pow(double(y - imSize/2), 2)));
window = 1;
value *= window;
interferogram[y * imSize + x] = (unsigned char)value;
}
}
The problem is that a char is between -128 and 127, while unsigned char ranges from 0 to 255. In interferogram[y * imSize + x] = (unsigned char)value;, there is an implicit cast to char.
It does not affect the output if window=0.5, but it triggers a change if window=1 as value becomes higher than 127. This is exactly the problem that you noticed in your question !
It does not affect the first displayed image since CV_8U corresponds to unsigned char : interferogram is therefore cast back into a unsigned char*. Take a look at Can I turn unsigned char into char and vice versa? to know more about char to unsigned char cast.
The problem occurs at input[i][0] = double(interferogram[i]) / 255.; : if window=1, interferogram[i] may be negative and input[i][0] becomes negative.
Change all char to unsigned char and it should solve the problem.
You may also change
outputReal[i] = output[i][0] / imN;
outputImaginary[i] = output[i][1];
for
outputReal[i] = output[i][0];
outputImaginary[i] = output[i][1];
Calls to fftw seems to be fine.