用 C++ 绘制图像的光谱（fftw，OpenCV）

Posted 2023-02-25

【中文标题】用 C++ 绘制图像的光谱（fftw，OpenCV）

【英文标题】：Drawing spectrum of an image in C++ (fftw, OpenCV)

【发布时间】：2011-07-31 15:45:15

【问题描述】：

我正在尝试创建一个程序来绘制给定图像的二维灰度光谱。我正在使用 OpenCV 和 FFTW 库。通过使用来自互联网的提示和代码并对其进行修改，我设法加载图像，计算该图像的 fft 并从 fft 重新创建图像（它是相同的）。我无法做的是绘制傅立叶光谱本身。请你帮助我好吗？ 这是代码（删除了不太重要的行）：

/* 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 */

问题出在“频谱”部分。我得到的不是频谱，而是一些噪音。我究竟做错了什么？非常感谢您的帮助。

【问题讨论】：

听起来像是一个缩放问题 - 您应该检查 FFT 输出幅度值的范围。

你能建议我该怎么做吗？从我在不同论坛上读到的内容来看，如果它是关于大小的，我会得到一个黑色图像（中间有一个白点）。谢谢你的回答。

由于您没有对幅度进行任何范围检查，因此如果它很大，当您将其分配给 8 位像素时，它将模 2 包装（它看起来像噪声）。这就是为什么我说你应该检查范围 - 例如添加另一个循环以找到最大和最小幅度，然后相应地缩放您的值，以便您知道它们适合 8 位范围。

对功率谱应用对数缩放，然后重新缩放数据以适应 [0,255]。

我检查了最高值（和最低值），然后将输出数据除以该值并乘以 255。我现在有一个带有一个白点的黑色图像。我读到这样的问题是当我们使用整数而不是浮点数时。但就我而言，如您所见，我将输出数据直接保存到图像 img2。中间没有变数。我现在做错了什么？

【参考方案1】：

您需要绘制频谱的大小。这是代码。

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);  

【参考方案2】：

您可以尝试执行 IFFT 步骤，看看您是否恢复了原始图像？然后，您可以逐步检查您的问题在哪里。另一个解决问题的方法是用你预定义的一个小矩阵来做这个过程，然后在 MATLAB 中对其进行 FFT 计算，然后逐步检查，它对我有用！