OpenCV + FFTW - 幅度图像

Question

你好。

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

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

我有什么问题

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

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

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

代码和我的想法

看看最重要的代码：（我在做什么是试图移植教程，因为它可以与 OpenCV 一起使用）

已编辑：（两个帖子合并） 行。所以我稍微修改了代码，但输出仍然与教程不同。 看一下代码：

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

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

在显示幅度时，通常的做法是实际显示幅度的对数，以便大值比小值更强烈地减少。
教程明确说明了这一点：“幅度看起来是黑色的，但不是。为了使信息可见，我们以对数方式缩放图像。”

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

Answer 2

您将计算值分配给 uchar 变量，您会失去精度，所有负值和高于 255 的值也会丢失。 尝试在实值变量中进行计算，然后将最终结果标准化为 0-255 范围，然后将其分配给 CV_8U 类型的结果图像。