title: OpenCV + FFTW - magnitude image
body: <h3>Hello again.</h3>
  <p>Today I'm working on extending my simple OpenCV image processing application. I would like to calculate phase and magnitude of my loaded cv::Mat. I have to use FFTW c++ library for this purpose (i know about dft in OpenCV).</p>
  <p>My work is based on tutorial: <a href="http://www.admindojo.com/discrete-fourier-transform-in-c-with-fftw/" rel="nofollow noreferrer">http://www.admindojo.com/discrete-fourier-transform-in-c-with-fftw/</a></p>
  <h3>What's my problem</h3>
  <p>So according to tutorial my output magnitude should be like:
  <img src="https://i.stack.imgur.com/iZrUw.jpg" alt="enter image description here" /></p>
  <p>Unfortunately my output is quite different:
  <img src="https://i.stack.imgur.com/P6TU9.png" alt="enter image description here" /></p>
  <p>On the other hand, image of phase is almost the same as tutorial image so this part is good.</p>
  <h3>Code and My thoughts</h3>
  <p>Take a look on the most important code: (what am I doing there is trying to port tutorial as it is to work with OpenCV)</p>
  <p>EDITED: (Both posts merged)
  Ok. So I changed code a bit, but output is still different from the tutorial.
  Take a look at the code:</p>
  <pre><code>void Processing::fft_moc(cv::Mat &amp;pixels, cv::Mat &amp;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 to real parts (values between 0 and MaxRGB)
  for( int x = 0; x &lt; pixels.rows; x++ )
  {
      for( int y = 0; y &lt; pixels.cols; y++ )
      {
          double blue = pixels.at&lt;cv::Vec3b&gt;(x,y)[0];
          double green = pixels.at&lt;cv::Vec3b&gt;(x,y)[1];
          double red = pixels.at&lt;cv::Vec3b&gt;(x,y)[2];
  
          // save as real numbers
          inB[squareSize*x+y][0] = blue;
          inG[squareSize*x+y][0] = green;
          inR[squareSize*x+y][0] = red;
      }
  }
  
  // perform FORWARD fft
  fftw_execute(planB);
  fftw_execute(planG);
  fftw_execute(planR);
  
  
  double ***outMagF=new double**[pixels.rows];
  for(int i = 0 ; i &lt; pixels.rows ; i++)
  {
      outMagF[i]=new double *[pixels.cols];
      for(int j = 0 ; j &lt; pixels.cols ; j++)
      {
          outMagF[i][j]= new double[3];
      }
  }
  
  //calculate magnitude
  //find min and max for 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 &lt; pixels.rows; x++ )
  {
      for( int y = 0; y &lt; pixels.cols; y++ )
      {
  
          int i = 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 &gt; magB ? magB : n_minB;
          n_maxB = n_maxB &lt; magB ? magB : n_maxB;
  
          n_minG = n_minG &gt; magG ? magG : n_minG;
          n_maxG = n_maxG &lt; magG ? magG : n_maxG;
  
          n_minR = n_minR &gt; magR ? magR : n_minR;
          n_maxR = n_maxR &lt; magR ? magR : n_maxR;
  
          outMagF[x][y][0] = magB;
          outMagF[x][y][1] = magG;
          outMagF[x][y][2] = magR;
      }
  }
  
  for( int x = 0; x &lt; pixels.rows; x++ )
  {
      for( int y = 0; y &lt; pixels.cols; y++ )
      {
          int i = 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 to output = (value-min)/(max-min)
          outMag.at&lt;cv::Vec3f&gt;(x,y)[0] = (double)(outMagF[x][y][0]-n_minB)/(n_maxB-n_minB);
          outMag.at&lt;cv::Vec3f&gt;(x,y)[1] = (double)(outMagF[x][y][1]-n_minG)/(n_maxG-n_minG);
          outMag.at&lt;cv::Vec3f&gt;(x,y)[2] = (double)(outMagF[x][y][2]-n_minR)/(n_maxR-n_minR);
  
          // std::complex for arg()
          std::complex&lt;double&gt; cB(realB, imagB);
          std::complex&lt;double&gt; cG(realG, imagG);
          std::complex&lt;double&gt; cR(realR, imagR);
  
          // phase
          double phaseB = arg(cB) + M_PI;
          double phaseG = arg(cG) + M_PI;
          double phaseR = arg(cR) + M_PI;
  
          // scale and write to output
          outPhase.at&lt;cv::Vec3f&gt;(x,y)[0] = (phaseB / (double)(2 * M_PI)) * 1;
          outPhase.at&lt;cv::Vec3f&gt;(x,y)[1] = (phaseG / (double)(2 * M_PI)) * 1;
          outPhase.at&lt;cv::Vec3f&gt;(x,y)[2] = (phaseR / (double)(2 * M_PI)) * 1;
      }
  }
  
  // move zero frequency to (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);
  }
  </code></pre>
  <p>I store the final output in cv::Mat with type CV_32FC3. And Yes, the way I normalize magnitude is quite ugly but I just wanted to be sure that everything is working like I expect.</p>
  <p>Take a look at my output again:</p>
  <p><img src="https://i.stack.imgur.com/GcKdf.jpg" alt="enter image description here" /></p>
  <p>So as You can see I still need help with that.</p>