OpenCV + FFTW -震级图像

Question

又见面了。

今天，我正在扩展我的简单OpenCV图像处理应用程序。我想计算我加载的cv::Mat的相位和大小。为此，我必须使用FFTW库(我知道OpenCV中的c++ )。

我的工作基于教程：http://www.admindojo.com/discrete-fourier-transform-in-c-with-fftw/

我的问题是什么

所以根据教程，我的输出幅度应该是这样的：

​

不幸的是，我的输出是完全不同的：

​

另一方面，相位图像与教程图像几乎相同，所以这部分很好。

代码和我的想法

看看最重要的代码：(我在那里做的是尝试移植教程，因为它是与OpenCV一起工作的)

编辑：(两篇文章合并) Ok。所以我稍微修改了一下代码，但输出结果仍然与教程有所不同。看一下代码：

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 to real parts (values between 0 and MaxRGB)
for( int x = 0; x < pixels.rows; x++ )
{
    for( int y = 0; y < pixels.cols; y++ )
    {
        double blue = pixels.at<cv::Vec3b>(x,y)[0];
        double green = pixels.at<cv::Vec3b>(x,y)[1];
        double red = pixels.at<cv::Vec3b>(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 < 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 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 < pixels.rows; x++ )
{
    for( int y = 0; y < 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 > 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 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<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 for 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 and write to 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 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);
}

我将最终输出存储在cv::Mat中，类型为CV_32FC3。是的，我归一化幅度的方式是相当丑陋的，但我只是想确保一切都像我预期的那样工作。

再看看我的输出：

​

​

因此，正如你所看到的，我仍然需要帮助。

Answer 1

您将计算值赋给uchar变量，您将失去精度，并且所有负值和大于255的值也会丢失。尝试在实值变量中进行计算，然后将最终结果规范化到0-255范围内，然后将其赋给CV_8U类型的结果图像。

Answer 2

FFT平面通常包含非常大的第0个元素( DC)和通常接近于零的其余元素之间的非常大的差异。

在显示震级时，通常的做法是实际显示震级的对数，以便较大的值比较小的值减少得更强烈。

本教程明确说明了这一点：“图像的大小看起来是黑色的，但实际上不是。为了使信息可见，我们对图像进行对数缩放。”

您需要显示值的对数才能看到相似的图像。