1

でFFTベースのサブピクセルシフト(変換)アルゴリズムを実装しようとしていPythonます。フーリエシフト定理により、アレイをサブピクセル量で次のように変換できます。1.順方向FFTアレイ2.フーリエ空間での線形位相ランプによるアレイの乗算3.逆方向FFTアレイ

このアルゴリズムは、numpy / scipyを使用してPythonで簡単に実装できますが、256 ** 2アレイの場合、シフトごとに非常に低速(〜10msec)になります。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_complex配列「i」にコピーする必要があることがわかります。そして最後に、結果「i」を出力のnumpy配列「b」にコピーする必要があります。fftwでnumpy配列「a」と「b」を直接使用する方がはるかに効率的ですが、これを機能させることはできません。

fftwが複雑なnumpy配列を直接使用する方法について誰かが考えていますscipy.weave.inlineか?

ありがとう

4

1 に答える 1

0

fftwのマニュアルによると、complex.h前にインポートできます。これにより、ネイティブCデータ型に対応するfftw.hことが保証されます。fftw_complexnumpyデータ型もネイティブCデータ型と互換性があることが保証されている(または実際には互換性がある可能性が高い)と確信しています。

この場合、配列データへのポインタにとしてアクセスできますa.data_as(ctypes.c_void_p)。残念ながら、ctypesは複雑な型を認識しませんが、うまくいけば、voidポインターにキャストすることでうまくいくでしょう。

これを行うときは、配列がC連続で格納され、配列の作成時aにパラメーターで指定されるように注意する必要があります。order='C'

于 2012-10-25T07:11:49.230 に答える