I am trying to use the 1-D Fourier Transform algorithm in chapter 12.2 to do 2-D image processing.

What we are trying to do is take an 8-bit grayscale image, input the image into a 2-D array. The original image is a PGM file (raw format) and the 2-D array used to store the original image data is a array of unsigned characters (0-255). The data is then cast into an array of floating point numbers.

Then this array is prepared by transferring it into an array that is twice the width of the original and the same height. The data is then transferred to this new array. The original image data becomes the real parts of the complex array, and the imaginary parts are initialized to zero. The result is an array with real and complex parts, with the imaginary components set to zero initially.

The next step in my program is to take each row of the complex array and copy it into a 1-D array of the same size. This 1-D array is then tranferred as an argument into the four1 function from numerical recipes. After the fourier transform is done on the array it is then copied back into the complex array. This process continues for all the rows.

Then all the columns are transformed in a similar manner. The only change with the columns is to take 2 columns (column 1 being real, and column 2 being imaginary) from every row and transferred them to the 1-D array.

The result is then copied into 2 separate arrays. The real and imaginary components each go to separate arrays. The magnitude is taken and is then output to an image file.

I have used a simple test image of a white box on a black background. The image of the frequency is almost correct, except for some distortion in the middle. This little bit of distortion is causing problems when a filter is applied. I can't figure out what the problem is.

Here is the code from my program


void fft(float **datafft, int xsize, int ysize, int dir)
/************************************************** *******/
/* Function fft */
/* Function performs an fft on the prepares array. Produces an array */
/* with real and imaginary parts. dir is used for the direction to */
/* perform the fourier transform. 1 for forward fourier, and -1 for */
/* inverse fourier transform. */
/************************************************** *******/
{
int i, j, k, l;
float *datainput;
datainput = (float *) calloc (2*xsize, sizeof (float));

for (i = 0; i < ysize; i++)
{
for (j = 0; j < 2*xsize; j++)
{
tempRI[i][j] = 0.0;
}
}


/* Do transform on all rows */
for (i = 0; i < ysize; i++)
{
for (j = 0; j < 2*xsize; j++)
{
/* Copy each row of prepared array into temp. 1-D array */
datainput[j] = datafft[i][j];
}
/* Input temporary array into 1-D fft function (four1) */
four1(datainput-1, xsize, dir);

/* if (dir == 1)
{
for (l = 0; l < 2*xsize; l++)
datainput[l]= datainput[l]/(xsize);
}
*/
for (j = 0; j < 2*xsize; j++)
{
/* Copy each row from temp. array back to temp. array */
imagefltin[i][j] = datainput[j];
}
}


/* Do transform on all columns */
k = 0;
for (j = 0; j < 2*xsize; j+=2)
{
for (i = 0; i < ysize; i++)
{
/* Copy each column into temporary 1-D array */
datainput[k] = imagefltin[i][j];
datainput[k+1] = imagefltin[i][j+1];
k+=2;
}


/* Input temp. array into 1-D fft function (four1) */
four1(datainput-1, xsize, dir);
/*
if (dir == 1)
{
for (l = 0; l < 2*xsize; l++)
datainput[l]= datainput[l]/(xsize);
}
*/
k = 0;
/* Copy each column from temp. array back to temp. array */
for (i = 0; i < ysize; i++)
{
tempRI[i][j] = datainput[k];
tempRI[i][j+1] = datainput[k+1];
k+=2;
}
k = 0;
}


k =0;
for (i = 0; i < ysize; i++)
{
for (j = 0; j < 2*xsize; j+=2)
{
realarray[i][k] = tempRI[i][j];
imaginaryarray[i][k] = tempRI[i][j+1];
k++;
}
k = 0;
}


}

/* End of Function fft */

Could you clarify a little bit what you regard as distortion? In other words, what is there about the result that strikes you as being wrong? The center of the transformed image corresponds to zero-frequency in both x and y directions. The image that you described (a white box on a black field) has a non-zero average value, and will therefore lead to a large DC component in the center of the image. If this is part of what you object to, it can be removed by adding a constant to the image to adjust its average to be zero.

Other than that, a simple box in a black field will have a fundamental frequency that is very low, and a lot of the image power spectrum will likely fall close to the center of the transformed image. That's as it should be.