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.

Answer 2

You did everything else right, but you need a logarithmic scale to see the effects of the typically small value changes of the fourier coefficients (like in your solution, there are small white domains in the output). The large values will then be white, the small changes will turn out to be what you are looking for. Even if you are not using OpenCV (here v. 2.4.2), here you can find a complete tutorial that describes it pretty cool in addition to theory: http://docs.opencv.org/doc/tutorials/core/discrete_fourier_transform/discrete_fourier_transform.html

I just saw your question. Maybe the answer comes too late, but it might help others in the future.. So could you please vote this answer as the right one to your question? I'm quite new in stackoverflow and need some reputation for complete participation ^^