我试图在Python
中实现一种基于FFT的亚像素移位(平移)算法。Fourier移位定理允许将数组按亚像素量转换为:
一。前向FFT阵列
2。Fourier空间中的线性相位斜率乘阵列
三。逆FFT阵列
这个算法很容易在python中使用numpy/scipy实现,但是对于256**2数组来说,它的每次移位速度非常慢(大约10毫秒)。我正试图通过使用scipy.weave.inline直接从python调用c代码来加速这一过程。
不过,我在将复杂的numpy数组传递给FFTW时遇到了问题。c代码看起来像:
#include <fftw3.h>
#include <stdlib.h>
#define INVERSE +1
#define FORWARD -1
fftw_complex *i, *o;
int n, m;
fftw_plan pf, pi;
#line 22 "test_scipy_weave.py"
i = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * xdim*ydim);
o = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * xdim*ydim);
pf = fftw_plan_dft_2d(xdim, ydim, i, o, -1, FFTW_PATIENT);
pi = fftw_plan_dft_2d(xdim, ydim, o, i, 1, FFTW_PATIENT);
# Copy data to fftw_complex array. How to use python arrays directly
for (n=0; n<xdim;n++){
for (m=0; m<ydim; m++){
i[n*xdim+m][0]=a[n*xdim+m].real();
i[n*xdim+m][1]=a[n*xdim+m].imag();
}
}
fftw_execute(pf);
/* Mult by linear phase ramp here */
fftw_execute(pi);
for (n=0; n<xdim;n++){
for (m=0; m<ydim; m++){
b[n*xdim+m] = std::complex<double>([in*xdim+m][0], i[n*xdim+m][1]);
}
}
fftw_destroy_plan(p);
所以你可以看到我必须将存储在numpy数组“a”中的数据复制到fftw_复杂数组“I”中。最后,我还要将结果“I”复制到输出numpy数组“b”中。直接在fftw中使用numpy数组“a”和“b”会更有效率,但我不能让它工作。
有人知道如何让fftw直接在
scipy.weave.inline
中使用复杂的numpy数组吗?谢谢
最佳答案
根据fftw manual,您可以在complex.h
之前导入fftw.h
,这将保证fftw_complex
与本机C数据类型相对应。我确信numpy数据类型也保证(或者实际上很可能)与本机C数据类型兼容。
在这种情况下,可以将指向数组数据的指针访问为a.data_as(ctypes.c_void_p)
。不幸的是,ctypes不能识别复杂的类型,但希望将其转换为一个空指针可以做到这一点。
执行此操作时,必须注意数组a
以C-连续方式存储,在创建数组时由参数order='C'
指定。