Weird but close fft and ifft of image in c++

Question

I wrote a program that loads, saves, and performs the fft and ifft on black and white png images. After much debugging headache, I finally got some coherent output only to find that it distorted the original image.

input:

fft:

ifft:

As far as I have tested, the pixel data in each array is stored and converted correctly. Pixels are stored in two arrays, 'data' which contains the b/w value of each pixel and 'complex_data' which is twice as long as 'data' and stores real b/w value and imaginary parts of each pixel in alternating indices. My fft algorithm operates on an array structured like 'complex_data'. After code to read commands from the user, here's the code in question:

if (cmd == "fft")
        {              
            if (height > width) size = height;
            else size = width;

            N = (int)pow(2.0, ceil(log((double)size)/log(2.0)));

            temp_data = (double*) malloc(sizeof(double) * width * 2); //array to hold each row of the image for processing in FFT()

            for (i = 0; i < (int) height; i++)
            {
                for (j = 0; j < (int) width; j++)
                {
                    temp_data[j*2] = complex_data[(i*width*2)+(j*2)];
                    temp_data[j*2+1] = complex_data[(i*width*2)+(j*2)+1];
                }
                FFT(temp_data, N, 1);
                for (j = 0; j < (int) width; j++)
                {
                    complex_data[(i*width*2)+(j*2)] = temp_data[j*2];
                    complex_data[(i*width*2)+(j*2)+1] = temp_data[j*2+1];
                }
            }
            transpose(complex_data, width, height); //tested
            free(temp_data);
            temp_data = (double*) malloc(sizeof(double) * height * 2);
            for (i = 0; i < (int) width; i++)
            {
                for (j = 0; j < (int) height; j++)
                {
                    temp_data[j*2] = complex_data[(i*height*2)+(j*2)];
                    temp_data[j*2+1] = complex_data[(i*height*2)+(j*2)+1];
                }
                FFT(temp_data, N, 1);
                for (j = 0; j < (int) height; j++)
                {
                    complex_data[(i*height*2)+(j*2)] = temp_data[j*2];
                    complex_data[(i*height*2)+(j*2)+1] = temp_data[j*2+1];
                }
            }
            transpose(complex_data, height, width);

            free(temp_data);
            free(data);

            data = complex_to_real(complex_data, image.size()/4); //tested
            image = bw_data_to_vector(data, image.size()/4); //tested
            cout << "*** fft success ***" << endl << endl;



void FFT(double* data, unsigned long nn, int f_or_b){ // f_or_b is 1 for fft, -1 for ifft

unsigned long n, mmax, m, j, istep, i;
double wtemp, w_real, wp_real, wp_imaginary, w_imaginary, theta;
double temp_real, temp_imaginary;

// reverse-binary reindexing to separate even and odd indices
// and to allow us to compute the FFT in place

n = nn<<1;
j = 1;
for (i = 1; i < n; i += 2) {
    if (j > i) {
        swap(data[j-1], data[i-1]);
        swap(data[j], data[i]);
    }
    m = nn;
    while (m >= 2 && j > m) {
        j -= m;
        m >>= 1;
    }
    j += m;
};

// here begins the Danielson-Lanczos section

mmax = 2;
while (n > mmax) {
    istep = mmax<<1;
    theta = f_or_b * (2 * M_PI/mmax);
    wtemp = sin(0.5 * theta);
    wp_real = -2.0 * wtemp * wtemp;
    wp_imaginary = sin(theta);
    w_real = 1.0;
    w_imaginary = 0.0;
    for (m = 1; m < mmax; m += 2) {
        for (i = m; i <= n; i += istep) {
            j = i + mmax;
            temp_real = w_real * data[j-1] - w_imaginary * data[j];
            temp_imaginary = w_real * data[j] + w_imaginary * data[j-1];

            data[j-1] = data[i-1] - temp_real;
            data[j] = data[i] - temp_imaginary;
            data[i-1] += temp_real;
            data[i] += temp_imaginary;
        }
        wtemp = w_real;
        w_real += w_real * wp_real - w_imaginary * wp_imaginary;
        w_imaginary += w_imaginary * wp_real + wtemp * wp_imaginary;
    }
    mmax=istep;
}}

My ifft is the same only with the f_or_b set to -1 instead of 1. My program calls FFT() on each row, transposes the image, calls FFT() on each row again, then transposes back. Is there maybe an error with my indexing?

Answer 1

Not an actual answer as this question is Debug only so some hints instead:

your results are really bad

it should look like this:

• first line is the actual DFFT result
• Re,Im,Power is amplified by a constant otherwise you would see a black image
• the last image is IDFFT of the original not amplified Re,IM result
• the second line is the same but the DFFT result is wrapped by half size of image in booth x,y to match the common results in most DIP/CV texts

As you can see if you IDFFT back the wrapped results the result is not correct (checker board mask)

You have just single image as DFFT result

• is it power spectrum?
• or you forget to include imaginary part? to view only or perhaps also to computation somewhere as well?

is your 1D **DFFT working?**

• for real data the result should be symmetric
• check the links from my comment and compare the results for some sample 1D array
• debug/repair your 1D FFT first and only then move to the next level
• do not forget to test Real and complex data ...

your IDFFT looks BW (no gray) saturated

• so did you amplify the DFFT results to see the image and used that for IDFFT instead of the original DFFT result?
• also check if you do not round to integers somewhere along the computation

beware of (I)DFFT overflows/underflows

If your image pixel intensities are big and the resolution of image too then your computation could loss precision. Newer saw this in images but if your image is HDR then it is possible. This is a common problem with convolution computed by DFFT for big polynomials.

Answer 2

Suggest you look at the article http://www.yolinux.com/TUTORIALS/C++MemoryCorruptionAndMemoryLeaks.html

Christophe has a good point but he is wrong about it not being related to the problem because it seems that in modern times using malloc instead of new()/free() does not initialise memory or select best data type which would result in all problems listed below:-

Possibly causes are:

1. Sign of a number changing somewhere, I have seen similar issues when a platform invoke has been used on a dll and a value is passed by value instead of reference. It is caused by memory not necessarily being empty so when your image data enters it will have boolean maths performed on its values. I would suggest that you make sure memory is empty before you put your image data there.

2. Memory rotating right (ROR in assembly langauge) or left (ROL) . This will occur if data types are being used which do not necessarily match, eg. a signed value entering an unsigned data type or if the number of bits is different in one variable to another.

3. Data being lost due to an unsigned value entering a signed variable. Outcomes are 1 bit being lost because it will be used to determine negative or positive, or at extremes if twos complement takes place the number will become inverted in meaning, look for twos complement on wikipedia.

Also see how memory should be cleared/assigned before use. http://www.cprogramming.com/tutorial/memory_debugging_parallel_inspector.html

Answer 3

Thank you everyone for your opinions. All that stuff about memory corruption, while it makes a point, is not the root of the problem. The sizes of data I'm mallocing are not overly large, and I am freeing them in the right places. I had a lot of practice with this while learning c. The problem was not the fft algorithm either, nor even my 2D implementation of it.

All I missed was the scaling by 1/(M*N) at the very end of my ifft code. Because the image is 512x512, I needed to scale my ifft output by 1/(512*512). Also, my fft looks like white noise because the pixel data was not rescaled to fit between 0 and 255.