c++ - OpenCV + FFTW - magnitude image -

c++ - OpenCV + FFTW - magnitude image -

hello again.

today i'm working on extending simple opencv image processing application. calculate phase , magnitude of loaded cv::mat. have utilize fftw c++ library purpose (i know dft in opencv).

my work based on tutorial: http://www.admindojo.com/discrete-fourier-transform-in-c-with-fftw/

what's problem

so according tutorial output magnitude should like:

unfortunately output quite different:

on other hand, image of phase same tutorial image part good.

code , thoughts

take on of import code: (what doing there trying port tutorial work opencv)

edited: (both posts merged) ok. changed code bit, output still different tutorial. take @ code:

void processing::fft_moc(cv::mat &pixels, cv::mat &outmag, cv::mat outphase, int mode) { int squaresize = pixels.cols; fftw_plan planr, plang, planb; fftw_complex *inr, *ing, *inb, *outr, *outg, *outb; // allocate input arrays inb = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squaresize * squaresize); ing = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squaresize * squaresize); inr = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squaresize * squaresize); // allocate output arrays outb = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squaresize * squaresize); outg = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squaresize * squaresize); outr = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * squaresize * squaresize); if (mode == fft) { // create plans planb = fftw_plan_dft_2d(squaresize, squaresize, inr, outb, fftw_forward, fftw_estimate); plang = fftw_plan_dft_2d(squaresize, squaresize, ing, outg, fftw_forward, fftw_estimate); planr = fftw_plan_dft_2d(squaresize, squaresize, inb, outr, fftw_forward, fftw_estimate); } // assig1n values real parts (values between 0 , maxrgb) for( int x = 0; x < pixels.rows; x++ ) { for( int y = 0; y < pixels.cols; y++ ) { double bluish = pixels.at<cv::vec3b>(x,y)[0]; double greenish = pixels.at<cv::vec3b>(x,y)[1]; double reddish = pixels.at<cv::vec3b>(x,y)[2]; // save real numbers inb[squaresize*x+y][0] = blue; ing[squaresize*x+y][0] = green; inr[squaresize*x+y][0] = red; } } // perform forwards fft fftw_execute(planb); fftw_execute(plang); fftw_execute(planr); double ***outmagf=new double**[pixels.rows]; for(int = 0 ; < pixels.rows ; i++) { outmagf[i]=new double *[pixels.cols]; for(int j = 0 ; j < pixels.cols ; j++) { outmagf[i][j]= new double[3]; } } //calculate magnitude //find min , max each channel double n_ming = 0.0; double n_maxg = 0.0; double n_minb = 0.0; double n_maxb = 0.0; double n_minr = 0.0; double n_maxr = 0.0; for( int x = 0; x < pixels.rows; x++ ) { for( int y = 0; y < pixels.cols; y++ ) { int = squaresize*x+y; // normalize values double realb = outb[i][0] / (double)(squaresize * squaresize); double imagb = outb[i][1] / (double)(squaresize * squaresize); double realg = outg[i][0] / (double)(squaresize * squaresize); double imagg = outg[i][1] / (double)(squaresize * squaresize); double realr = outr[i][0] / (double)(squaresize * squaresize); double imagr = outr[i][1] / (double)(squaresize * squaresize); // magnitude double magb = log(1+sqrt((realb * realb) + (imagb * imagb))); double magg = log(1+sqrt((realg * realg) + (imagg * imagg))); double magr = log(1+sqrt((realr * realr) + (imagr * imagr))); n_minb = n_minb > magb ? magb : n_minb; n_maxb = n_maxb < magb ? magb : n_maxb; n_ming = n_ming > magg ? magg : n_ming; n_maxg = n_maxg < magg ? magg : n_maxg; n_minr = n_minr > magr ? magr : n_minr; n_maxr = n_maxr < magr ? magr : n_maxr; outmagf[x][y][0] = magb; outmagf[x][y][1] = magg; outmagf[x][y][2] = magr; } } for( int x = 0; x < pixels.rows; x++ ) { for( int y = 0; y < pixels.cols; y++ ) { int = squaresize*x+y; double realb = outb[i][0] / (double)(squaresize * squaresize); double imagb = outb[i][1] / (double)(squaresize * squaresize); double realg = outg[i][0] / (double)(squaresize * squaresize); double imagg = outg[i][1] / (double)(squaresize * squaresize); double realr = outr[i][0] / (double)(squaresize * squaresize); double imagr = outr[i][1] / (double)(squaresize * squaresize); // write normalized output = (value-min)/(max-min) outmag.at<cv::vec3f>(x,y)[0] = (double)(outmagf[x][y][0]-n_minb)/(n_maxb-n_minb); outmag.at<cv::vec3f>(x,y)[1] = (double)(outmagf[x][y][1]-n_ming)/(n_maxg-n_ming); outmag.at<cv::vec3f>(x,y)[2] = (double)(outmagf[x][y][2]-n_minr)/(n_maxr-n_minr); // std::complex arg() std::complex<double> cb(realb, imagb); std::complex<double> cg(realg, imagg); std::complex<double> cr(realr, imagr); // phase double phaseb = arg(cb) + m_pi; double phaseg = arg(cg) + m_pi; double phaser = arg(cr) + m_pi; // scale , write output outphase.at<cv::vec3f>(x,y)[0] = (phaseb / (double)(2 * m_pi)) * 1; outphase.at<cv::vec3f>(x,y)[1] = (phaseg / (double)(2 * m_pi)) * 1; outphase.at<cv::vec3f>(x,y)[2] = (phaser / (double)(2 * m_pi)) * 1; } } // move 0 frequency (squaresize/2, squaresize/2) swapquadrants(squaresize, outmag); swapquadrants(squaresize, outphase); // free memory fftw_destroy_plan(planr); fftw_destroy_plan(plang); fftw_destroy_plan(planb); fftw_free(inr); fftw_free(outr); fftw_free(ing); fftw_free(outg); fftw_free(inb); fftw_free(outb); }

i store final output in cv::mat type cv_32fc3. , yes, way normalize magnitude quite ugly wanted sure working expect.

take @ output again:

so can see still need help that.

fft planes contain big difference between 0th elements (the dc) large, , rest of elements close zero.

when displaying magnitude mutual practice show log of magnitude big values cut down more little ones. tutorial states explicitly: "the magnitude appears black isn’t. create info visible scale image logarithmically."

you need display log of values see similar image.

c++ image opencv fftw