0

我想要做的过程是对图像进行 FFT（存储在“imagen”中），然后将其与过滤器“H”相乘，之后也会进行逆 FFT。代码如下所示：

int ancho;
int alto;
ancho=ui.imageframe->imagereader->GetBufferedRegion().GetSize()[0];     //ancho=widht of the image
alto=ui.imageframe->imagereader->GetBufferedRegion().GetSize()[1];      //alto=height of the image

double *H ;
H =matrix2D_H(ancho,alto,eta,sigma); // H is calculated

// We want to get: F= fft(f) ; H*F ; f'=ifft(H*F)
// Inicialization of the neccesary elements for the calculation of the fft
fftw_complex *out;
fftw_plan p;

int N= (ancho/2+1)*alto; //number of points of the image
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*N);

double *in = (double*) imagen.GetPointer(); // conversion of itk.smartpointer --> double*
p = fftw_plan_dft_r2c_2d(ancho, alto, in, out, FFTW_ESTIMATE); // FFT planning
fftw_execute(p); // FFT calculation

/* Multiplication of the Output of the FFT with the Filter H*/ 
int a = alto;
int b = ancho/2 +1; // The reason for the second dimension to have this value is that when the FFT calculation of a real image is performed only the non-redundants outputs are calculated, that’s the reason for the output of the FFT and the filter ‘H’ to be equal. 

// Matrix point-by-point multiplicaction: [axb]*[axb]
fftw_complex* res ; // result will be stored here
res = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*a*b);
res = multiply_matrix_2D(out,H, a, b);

问题出在此处，在函数“multiply_matrix_2D”内的循环中：

 fftw_complex*  prueba_r01::multiply_matrix_2D(fftw_complex* out, double* H, int M ,int N){
/* The matrix out[MxN] or [n0x(n1/2)+1] is the image after the FFT , and the out_H[MxN] is the filter in the frequency domain,
both are multiplied POINT TO POINT, it has to be called  twice, one for the imaginary part and another for the normal part
*/
fftw_complex *H_cast;
H_cast = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
H_cast= reinterpret_cast<fftw_complex*> (H); // casting from double* to fftw_complex*

fftw_complex *res; // the result of the multiplication will be stored here
res = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);

//Loop for calculating the matrix point-to-point multiplication
for (int x = 0; x<M ; x++){
    for (int y = 0; y<N ; y++){
            res[x*N+y][0] = out[x*N+y][0]*(H_cast[x*N+y][0]+H_cast[x*N+y][1]); 
            res[x*N+y][1] = out[x*N+y][1]*(H_cast[x*N+y][0]+H_cast[x*N+y][1]); 
        }
}
fftw_free(H_cast);
return res;
}

x = 95 和 y = 93 的值是 M = 191 和 N = 96；prueba_r01.exe 中 0x004273ab 处的不受控制的异常：0xC0000005 访问违规读取 0x01274000。

图片 http://img846.imageshack.us/img846/4585/accessviolationproblem.png

很多变量的值都是红色的，并且对于翻译问题：H_cast[][1] 在值框中：“Error30CXX0000：不可能评估表达式”。

我将非常感谢任何形式的帮助！

安东尼奥

4

2 回答 2

1

ITK 内部有一个 FFT 类，可以使用 cmake 中的 fftw (USE_FFTW) 进行配置。此类描述如何从 fftw 引用 ITK 原始缓冲内存。

PS：即将推出的 ITKv4 对 fftw 的兼容性有了很大的提升。

于 2011-09-23T02:26:24.903 回答
1

这部分代码

H_cast = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*M*N);
H_cast= reinterpret_cast<fftw_complex*> (H); // casting from double* to fftw_complex*

首先为其分配一个新缓冲区H_cast，然后立即将其设置为指向原始缓冲区H。它不复制数据，只是复制指针。

在函数结束时，一些缓冲区被释放

fftw_free(H_cast);

这似乎释放了指向的数据，H而不是函数中分配的缓冲区。

当回到来电者H那里时，那里迷路了！

于 2011-07-18T16:06:40.053 回答