CODEWITHSUNDEEP
opencvfftwspectrum
Narren
Narren

Reputation: 21

Drawing spectrum of an image in C++ (fftw, OpenCV)

I'm trying to create a program that will draw a 2d greyscale spectrum of a given image. I'm using OpenCV and FFTW libraries. By using tips and codes from the internet and modifying them I've managed to load an image, calculate fft of this image and recreate the image from the fft (it's the same). What I'm unable to do is to draw the fourier spectrum itself. Could you please help me? Here's the code (less important lines removed):

/* Copy input image */

/* Create output image */

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

/* Create plans */
plan_f = fftw_plan_dft_2d(w, h, in, dft, FFTW_FORWARD, 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.;
    }
}

/* forward DFT */
fftw_execute(plan_f);

/* spectrum */
for (i = 0, k = 0; i < h; i++)
{
    for (j = 0; j < w; j++, k++)
        ((uchar*)(img2->imageData + i * img2->widthStep))[j] = sqrt(pow(dft[k][0],2) + pow(dft[k][1],2));
}       

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

/* Free memory */

}

The problem is in the "Spectrum" section. Instead of a spectrum I get some noise. What am I doing wrong? I would be grateful for your help.

Upvotes: 2

Views: 4380

Answers (2)

mrgloom
mrgloom

Reputation: 21682

You need to draw magnitude of spectrum. here is the code.

void ForwardFFT(Mat &Src, Mat *FImg)
{
    int M = getOptimalDFTSize( Src.rows );
    int N = getOptimalDFTSize( Src.cols );
    Mat padded;    
    copyMakeBorder(Src, padded, 0, M - Src.rows, 0, N - Src.cols, BORDER_CONSTANT, Scalar::all(0));
    // Создаем комплексное представление изображения
    // planes[0] содержит само изображение, planes[1] его мнимую часть (заполнено нулями)
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexImg;
    merge(planes, 2, complexImg); 
    dft(complexImg, complexImg);    
    // После преобразования результат так-же состоит из действительной и мнимой части
    split(complexImg, planes);

    // обрежем спектр, если у него нечетное количество строк или столбцов
    planes[0] = planes[0](Rect(0, 0, planes[0].cols & -2, planes[0].rows & -2));
    planes[1] = planes[1](Rect(0, 0, planes[1].cols & -2, planes[1].rows & -2));

    Recomb(planes[0],planes[0]);
    Recomb(planes[1],planes[1]);
    // Нормализуем спектр
    planes[0]/=float(M*N);
    planes[1]/=float(M*N);
    FImg[0]=planes[0].clone();
    FImg[1]=planes[1].clone();
}
void ForwardFFT_Mag_Phase(Mat &src, Mat &Mag,Mat &Phase)
{
    Mat planes[2];
    ForwardFFT(src,planes);
    Mag.zeros(planes[0].rows,planes[0].cols,CV_32F);
    Phase.zeros(planes[0].rows,planes[0].cols,CV_32F);
    cv::cartToPolar(planes[0],planes[1],Mag,Phase);
}
Mat LogMag;
    LogMag.zeros(Mag.rows,Mag.cols,CV_32F);
    LogMag=(Mag+1);
    cv::log(LogMag,LogMag);
    //---------------------------------------------------
    imshow("Логарифм амплитуды", LogMag);
    imshow("Фаза", Phase);
    imshow("Результат фильтрации", img);

Upvotes: 2

Antonio
Antonio

Reputation: 579

Can you try to do the IFFT step and see if you recover the original image ? then , you can check step by step where is your problem. Another solution to find the problem is to do this process with a small matrix predefined by you ,and calculate it FFT in MATLAB, and check step by step, it worked for me!

Upvotes: 0

Related Questions