向前FFT图像和向后FFT图像得到相同的结果

Forward FFT an image and backward FFT an image to get the same result

本文关键字:图像 FFT 结果 向前      更新时间:2023-10-16

我正在尝试使用来自http://www.fftw.org/的库来FFT图像,以便我可以在频域中进行卷积。但我不知道该怎么做。为了理解如何做到这一点,我试图将图像作为pixelcolors数组转发FFT,然后向后FFT以获得相同的pixelcolors数组。我是这样做的:

fftw_plan planR, planG, planB;
fftw_complex *inR, *inG, *inB, *outR, *outG, *outB, *resultR, *resultG, *resultB;
//Allocate arrays.
inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        inR[y * width + x][0] = pixelColors[currentIndex];
        inG[y * width + x][0] = pixelColors[currentIndex + 1];
        inB[y * width + x][0] = pixelColors[currentIndex + 2];
    }
}
//Forward plans.
planR = fftw_plan_dft_2d(width, width, inR, outR, FFTW_FORWARD, FFTW_MEASURE);
planG = fftw_plan_dft_2d(width, width, inG, outG, FFTW_FORWARD, FFTW_MEASURE);
planB = fftw_plan_dft_2d(width, width, inB, outB, FFTW_FORWARD, FFTW_MEASURE);
//Forward FFT.
fftw_execute(planR);
fftw_execute(planG);
fftw_execute(planB);
//Backward plans.
planR = fftw_plan_dft_2d(width, width, outR, resultR, FFTW_BACKWARD, FFTW_MEASURE);
planG = fftw_plan_dft_2d(width, width, outG, resultG, FFTW_BACKWARD, FFTW_MEASURE);
planB = fftw_plan_dft_2d(width, width, outB, resultB, FFTW_BACKWARD, FFTW_MEASURE);
//Backward fft
fftw_execute(planR);
fftw_execute(planG);
fftw_execute(planB);
//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        pixelColors[currentIndex] = resultR[y * width + x][0];
        pixelColors[currentIndex + 1] = resultG[y * width + x][0];
        pixelColors[currentIndex + 2] = resultB[y * width + x][0];
    }
}

有人能告诉我一个例子,如何转发FFT图像,然后向后FFT使用FFTW图像得到相同的结果?我一直在看很多展示如何使用FFTW到FFT的例子,但我无法弄清楚它如何适用于我的情况,我有一个代表图像的pixelcolors数组。

需要注意的一件重要的事情是,当你在反向FFT之后进行正向FFT时,这通常会导致N的缩放因子被应用到最终结果中,即得到的图像像素值需要除以N才能匹配原始像素值。(N是FFT的大小。)所以你的输出循环应该看起来像这样:

//Overwrite the pixelcolors with the result.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        pixelColors[currentIndex] = resultR[y * width + x][0] / (width * height);
        pixelColors[currentIndex + 1] = resultG[y * width + x][0] / (width * height);
        pixelColors[currentIndex + 2] = resultB[y * width + x][0] / (width * height);
    }
}

还要注意,您可能想要做一个从真到复杂的FFT,然后是一个从复杂到真的IFFT(在内存和性能方面都更有效)。现在看起来你在两个方向上都是复数到复数,这很好,但是你没有正确地填充输入数组。如果你坚持使用复杂到复杂的格式,那么你可能需要把你的输入循环改成这样:

//Fill in arrays with the pixelcolors.
for (int y = 0; y < height; y++) {
    for (int x = 0; x < width; x++) {
        int currentIndex = ((y * width) + (x)) * 3;
        inR[y * width + x][0] = (double)pixelColors[currentIndex];
        inR[y * width + x][1] = 0.0;
        inG[y * width + x][0] = (double)pixelColors[currentIndex + 1];
        inG[y * width + x][1] = 0.0;
        inB[y * width + x][0] = (double)pixelColors[currentIndex + 2];
        inB[y * width + x][1] = 0.0;
    }
}

。像素值进入复数输入值的实部,虚部需要归零。

还有一件事要注意:当您最终使其工作时,您会发现性能非常糟糕-相对于实际FFT所花费的时间,创建计划需要花费很长时间。其思想是,您只创建一次计划,但使用它来执行许多fft。因此,您需要将计划创建从实际的FFT代码中分离出来,并将其放入初始化例程或构造函数或其他任何东西中。

但是如果你使用realToComplex或ComplexToRealFunction,请注意,图像将存储在维度的矩阵中[height x (width/2 +1)],如果你想在频域中做一些中间计算,它们将变得有点困难…

它不起作用的原因是fftw_plan_dft_2d()做了一些基准测试以找到最佳算法并在过程中更改输入数据,因此您必须在fftw_plan_dft_2d()之后填充输入数据,而不是在它之前。