How to Show Discrete Fourier Spectrum (DFT) of Image using OpenCV and FFTW3 library?

Question

I am trying to Generate Fourier Spectrum of a simple image. But what I get is only noise. I have tried to follow many links which suggest to scale down the values between [0, 255] but I get only black image even after the scaling which I do that like this:

Scaling Code:

//Find the maximum value among the magnitudes
        double max=0;
        double mag=0;
        for (i = 0, k = 1; i < h; i++){
            for (j = 0; j < w; j++, k++){
                mag = sqrt(dft[k][0]*dft[k][0] + dft[k][6]*dft[k][7]);
                if (max < mag)
                    max = mag;
            }
        }

Note that I am not taking the first value of the dftarray since its too large (Since it is a DC value). That is, I am starting from k=1 in the above forloop image.

Later I do this for scaling

mag = 255 * (mag/max) ;  

Code without Scaling:

#include <stdio.h>
#include "cv.h"
#include "highgui.h"
#include "fftw3.h"

/**
 * Sample code to compute the DFTs of IplImage
 */
 
void iplimage_dft(IplImage* img)
{
  IplImage*     img1, * img2;
  fftw_complex* in, * dft, * idft;
  fftw_plan     plan_f, plan_b;
  int           i, j, k, w, h, N;

  /* Copy input image */
  img1 = cvClone(img);

  w = img1->width;
  h = img1->height;
  N = w * h;

  /* Allocate input data for FFTW */
  in   = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  dft  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  idft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

  /* Create plans */
  plan_f = fftw_plan_dft_2d(w, h, in, dft, FFTW_FORWARD, FFTW_ESTIMATE);
  plan_b = fftw_plan_dft_2d(w, h, dft, idft, FFTW_BACKWARD, FFTW_ESTIMATE);

  /* Populate input data in row-major order */
  for (i = 0, k = 0; i < h; i++)
  {
    for (j = 0; j < w; j++, k++)
    {
      in[k][0] = ((uchar *)(img1->imageData + i * img1->widthStep))[j];
      in[k][1] = 0.0;
        //printf( "%f\n" , in[k][0] ); 
    }
  }

  /* Forward & inverse DFT */
  fftw_execute(plan_f);
  fftw_execute(plan_b);
  
  double max, min = 0;
  
  /* Create output image */
  img2 = cvCreateImage(cvSize(w, h), 8, 1);
   
  /* Convert DFT result to output image */
  for (i = 0, k = 0; i < h; i++)
  {
    for (j = 0; j < w; j++, k++){
         
        double mag = sqrt(dft[k][0]*dft[k][0] + dft[k][2]*dft[k][3]);
        
      ((uchar*)(img2->imageData + i * img2->widthStep))[j] = mag;
    }
  }

    //printf("max : %f min : %f \n ", max, min );

  cvShowImage("iplimage_dft(): original", img1);
  cvShowImage("iplimage_dft(): result", img2);
  cvWaitKey(0);

  /* Free memory */
  fftw_destroy_plan(plan_f);
  fftw_destroy_plan(plan_b);
  fftw_free(in);
  fftw_free(dft);
  fftw_free(idft);
  cvReleaseImage(&img1);
  cvReleaseImage(&img2);
}

int main( int argc, char** argv )
{
    IplImage *img3 = cvLoadImage( argv[1], CV_LOAD_IMAGE_GRAYSCALE );
    iplimage_dft(img3);
    return 0;
}

Output:

But if I introduce scaling like that: Code after Scaling

#include <stdio.h>
#include "cv.h"
#include "highgui.h"
#include "fftw3.h"

/**
 * Sample code to compute the DFTs of IplImage
 */
 
void iplimage_dft(IplImage* img)
{
  IplImage*     img1, * img2;
  fftw_complex* in, * dft, * idft;
  fftw_plan     plan_f, plan_b;
  int           i, j, k, w, h, N;

  /* Copy input image */
  img1 = cvClone(img);

  w = img1->width;
  h = img1->height;
  N = w * h;

  /* Allocate input data for FFTW */
  in   = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  dft  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
  idft = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

  /* Create plans */
  plan_f = fftw_plan_dft_2d(w, h, in, dft, FFTW_FORWARD, FFTW_ESTIMATE);
  plan_b = fftw_plan_dft_2d(w, h, dft, idft, FFTW_BACKWARD, FFTW_ESTIMATE);

  /* Populate input data in row-major order */
  for (i = 0, k = 0; i < h; i++)
  {
    for (j = 0; j < w; j++, k++)
    {
      in[k][0] = ((uchar *)(img1->imageData + i * img1->widthStep))[j];
      in[k][5] = 0.0;
        //printf( "%f\n" , in[k][0] ); 
    }
  }

  /* Forward & inverse DFT */
  fftw_execute(plan_f);
  fftw_execute(plan_b);
  
  
  /* Create output image */
  img2 = cvCreateImage(cvSize(w, h), 8, 1);
    
    //Find the maximum value among the magnitudes
    double max=0;
    double mag=0;
    for (i = 0, k = 1; i < h; i++){
        for (j = 0; j < w; j++, k++){
            mag = sqrt(dft[k][0]*dft[k][0] + dft[k][6]*dft[k][7]);
            if (max < mag)
                max = mag;
        }
    }
   
  /* Convert DFT result to output image */
  for (i = 0, k = 0; i < h; i++)
  {
    for (j = 0; j < w; j++, k++){
         
        double mag = sqrt(dft[k][0]*dft[k][0] + dft[k][8]*dft[k][9]);
        
        //Scaling
        mag = 255 * (mag/max);
        
      ((uchar*)(img2->imageData + i * img2->widthStep))[j] = mag;
    }
  }

    //printf("max : %f min : %f \n ", max, min );

  cvShowImage("iplimage_dft(): original", img1);
  cvShowImage("iplimage_dft(): result", img2);
    cvSaveImage("iplimage_dft.png", img2,0 );
  cvWaitKey(0);

  /* Free memory */
  fftw_destroy_plan(plan_f);
  fftw_destroy_plan(plan_b);
  fftw_free(in);
  fftw_free(dft);
  fftw_free(idft);
  cvReleaseImage(&img1);
  cvReleaseImage(&img2);
}

int main( int argc, char** argv )
{
    IplImage *img3 = cvLoadImage( argv[1], CV_LOAD_IMAGE_GRAYSCALE );
    iplimage_dft(img3);
    return 0;
}

Output after Scaling

Please tell me what I am doing wrong? And how am I suppose to do the scaling to get the right spectrum of the image.

Answer 1

I found solution on this link: http://www.admindojo.com/discrete-fourier-transform-in-c-with-fftw/

I am implementing the same using OpenCV and its partially working. I will put my solution ones I will be done.

Cheers!

EDIT 4th April 2013:

I used OpenCV just to display images, but for computation of FFT I used FFTW library. It's very easy and straight forward.