在C++中进行图像的奇怪但紧密的FFT和IFFT

4
我写了一个程序,可以加载、保存和对黑白png图像执行fft和ifft操作。经过长时间的调试,我终于得到了一些连贯的输出,但发现它扭曲了原始图像。
输入: input fft: fft ifft: enter image description here 据我测试,每个数组中的像素数据都被正确存储和转换。像素以两个数组的形式存储:'data'包含每个像素的黑白值,而'complex_data'则是'data'长度的两倍,并在交替索引中存储每个像素的实部和虚部的黑白值。我的fft算法在一个类似'complex_data'的数组上运行。在从用户读取命令的代码之后,这里是相关代码:
if (cmd == "fft")
        {              
            if (height > width) size = height;
            else size = width;

            N = (int)pow(2.0, ceil(log((double)size)/log(2.0)));

            temp_data = (double*) malloc(sizeof(double) * width * 2); //array to hold each row of the image for processing in FFT()

            for (i = 0; i < (int) height; i++)
            {
                for (j = 0; j < (int) width; j++)
                {
                    temp_data[j*2] = complex_data[(i*width*2)+(j*2)];
                    temp_data[j*2+1] = complex_data[(i*width*2)+(j*2)+1];
                }
                FFT(temp_data, N, 1);
                for (j = 0; j < (int) width; j++)
                {
                    complex_data[(i*width*2)+(j*2)] = temp_data[j*2];
                    complex_data[(i*width*2)+(j*2)+1] = temp_data[j*2+1];
                }
            }
            transpose(complex_data, width, height); //tested
            free(temp_data);
            temp_data = (double*) malloc(sizeof(double) * height * 2);
            for (i = 0; i < (int) width; i++)
            {
                for (j = 0; j < (int) height; j++)
                {
                    temp_data[j*2] = complex_data[(i*height*2)+(j*2)];
                    temp_data[j*2+1] = complex_data[(i*height*2)+(j*2)+1];
                }
                FFT(temp_data, N, 1);
                for (j = 0; j < (int) height; j++)
                {
                    complex_data[(i*height*2)+(j*2)] = temp_data[j*2];
                    complex_data[(i*height*2)+(j*2)+1] = temp_data[j*2+1];
                }
            }
            transpose(complex_data, height, width);

            free(temp_data);
            free(data);

            data = complex_to_real(complex_data, image.size()/4); //tested
            image = bw_data_to_vector(data, image.size()/4); //tested
            cout << "*** fft success ***" << endl << endl;



void FFT(double* data, unsigned long nn, int f_or_b){ // f_or_b is 1 for fft, -1 for ifft

unsigned long n, mmax, m, j, istep, i;
double wtemp, w_real, wp_real, wp_imaginary, w_imaginary, theta;
double temp_real, temp_imaginary;

// reverse-binary reindexing to separate even and odd indices
// and to allow us to compute the FFT in place

n = nn<<1;
j = 1;
for (i = 1; i < n; i += 2) {
    if (j > i) {
        swap(data[j-1], data[i-1]);
        swap(data[j], data[i]);
    }
    m = nn;
    while (m >= 2 && j > m) {
        j -= m;
        m >>= 1;
    }
    j += m;
};

// here begins the Danielson-Lanczos section

mmax = 2;
while (n > mmax) {
    istep = mmax<<1;
    theta = f_or_b * (2 * M_PI/mmax);
    wtemp = sin(0.5 * theta);
    wp_real = -2.0 * wtemp * wtemp;
    wp_imaginary = sin(theta);
    w_real = 1.0;
    w_imaginary = 0.0;
    for (m = 1; m < mmax; m += 2) {
        for (i = m; i <= n; i += istep) {
            j = i + mmax;
            temp_real = w_real * data[j-1] - w_imaginary * data[j];
            temp_imaginary = w_real * data[j] + w_imaginary * data[j-1];

            data[j-1] = data[i-1] - temp_real;
            data[j] = data[i] - temp_imaginary;
            data[i-1] += temp_real;
            data[i] += temp_imaginary;
        }
        wtemp = w_real;
        w_real += w_real * wp_real - w_imaginary * wp_imaginary;
        w_imaginary += w_imaginary * wp_real + wtemp * wp_imaginary;
    }
    mmax=istep;
}}

我的ifft与f_or_b设置为-1而不是1相同。我的程序在每行上调用FFT(),转置图像,再次在每行上调用FFT(),然后转置回来。也许我的索引有错误吗?


3
与您的问题无关,但您应该避免在C++中使用malloc/free,而应该使用new/delete。 - Christophe
从编写C代码的旧习惯。 - Justa Ghost
1
当您使用非 POD 数据时,危险的习惯就会出现。 - Christophe
请查看以下链接,了解如何计算DFFT:https://dev59.com/questions/toTba4cB1Zd3GeqP3T49#26355569。在子链接中,您将找到我编写的C++源代码,用于(I)DFT/DFFT/DCT/DFCT,因此您可以将其用于调试(与您的结果进行比较)。 - Spektre
3个回答

2

由于这个问题只是调试,所以以下是一些提示而不是实际答案:

你的结果非常糟糕

应该像这样:

good results

  • 第一行是实际的DFFT结果
  • Re,Im,Power被一个常数放大,否则你会看到一个黑色图像
  • 最后一张图片是原始未放大的Re,IM结果的IDFFT
  • 第二行是相同的,但是DFFT结果在x,y中都被半尺寸包裹起来,以匹配大多数DIP/CV文本中的常见结果

如您所见,如果您反向IDFFT包装结果,则结果不正确(棋盘掩模)

您只有单个图像作为DFFT结果

  • 这是功率谱吗?
  • 或者您忘记了包含虚部?只是用于显示还是也在计算中使用了它?

您的1D **DFFT是否有效

  • 对于实际数据,结果应该是对称的
  • 检查我的评论中的链接,并比较一些样本1D数组的结果
  • 先调试/修复您的1D FFT,然后再进入下一个级别
  • 不要忘记测试实数和复杂数据...

您的IDFFT看起来饱和(没有灰度)

  • 所以你是否放大了DFFT结果以查看图像,并将其用于IDFFT而不是原始的DFFT结果?
  • 还要检查是否在计算过程中四舍五入到整数

注意(I)DFFT的溢出/下溢

如果您的图像像素强度很大且图像分辨率也很大,则您的计算可能会失去精度。虽然从未在图像中看到过这个问题,但如果您的图像是HDR,则有可能发生。这是使用DFFT计算大多项式进行卷积的常见问题。


2
谢谢大家提供的意见。虽然关于内存损坏的问题很有道理,但它并不是问题的根源。我分配内存的数据大小并不是过度大的,并且我也在正确的位置释放内存。在学习C语言时,我练习了很多次这个问题。问题也不是FFT算法,甚至不是2D实现问题。
我错过的只是在我的IFFT代码末尾进行 1/(M*N) 的比例缩放。因为图像是512x512,所以我需要将IFFT输出按1 /(512 * 512)进行比例缩放。另外,我的FFT看起来像白噪声,因为像素数据没有重新缩放以适合0到255之间的范围。

你应该接受你的答案(点击左侧投票计数器旁边的检查图标,以便其他人可以清楚快速地看到你的问题已解决(以及如何解决)。 - Spektre

0
建议您查看文章http://www.yolinux.com/TUTORIALS/C++MemoryCorruptionAndMemoryLeaks.html Christophe提出了一个很好的观点,但他错了,因为在现代时代中使用malloc而不是new()/free()不会初始化内存或选择最佳数据类型,这将导致以下所有问题:
可能的原因是:
  1. 如果数字的符号在某处改变,我曾经看到过类似的问题,当使用平台调用调用dll并且传递了一个按值而非参考值时。这是由于内存不一定为空,因此当您的图像数据进入时,它的值将进行布尔数学运算。我建议您在将图像数据放置在内存之前确保内存为空。

  2. 内存向右旋转(汇编语言中的ROR)或向左旋转(ROL)。如果正在使用不一定匹配的数据类型,例如,有符号值进入无符号数据类型,或者如果一个变量的位数与另一个变量不同,则会发生这种情况。

  3. 由于无符号值进入有符号变量而导致数据丢失。结果是1个比特被丢失,因为它将用于确定负或正,或者在极端情况下,如果进行二补码计算,则该数字的含义将被反转,请查看维基百科上的二补码。

另外,请查看如何在使用之前清除/分配内存。http://www.cprogramming.com/tutorial/memory_debugging_parallel_inspector.html


将来,请“编辑”您的答案以添加更多信息,而不是在同一个问题上发布多个答案。 - Matt

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接