画像にFFT

画像に2次元のフーリエ変換を行うには、縦、横の各ラインに対して1次元のFFTを行うことでできます。手順は、以下の図のようになります。

1次元のFFTの関数が与えられてるとして、画像にFFTを適用するには以下のようにします。

(FFTの関数にfftw3ライブラリを使用した場合）

#include <fftw3.h>

/* 反時計回りに行列を転置 */
void Transverse_matrix_L( double *a, double *b, int nX, int nY )
{
    int i, j;
    
    for ( i = 0; i < nY; i++ )
    {
        for ( j = 0; j < nX; j++ )
        {
            b[j][i] = a[i][j];
        }
    }
}

/* 時計回りに行列を転置 */
void Transverse_matrix_R( double *a, double *b, int nX, int nY )
{
    int i, j;
    
    for ( i = 0; i < nY; i++ )
    {
        for ( j = 0; j < nX; j++ )
        {
            b[i][j] = a[j][i];
        }
    }
}

/* 2次元フーリエ変換 */
void fft2( double *a_re, double *a_im, int nX, int nY )
{
    fftw_complex *in, *out;
    fftw_plan p;
    int i, j;
    double *w_re, *w_im;
    
    /* 領域の確保 */
    in = fftw_malloc( nX * sizeof(fftw_complex) );
    out = fftw_malloc( nX * sizeof(fftw_complex) );

    /* 横方向のFFT */
    for ( i = 0; i < nY; i++ )
    {
        /* inの実部に画像を詰め込む */
        for ( j = 0; j < nX; j++ )
        {
            in[j][0] = a_re[i*nX+j];
        }
        
        /* プランの作成 */
        p = fftw_plan_dft_1d(nX, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
        
        /* FFT実行 */
        fftw_execute(p);
        
        /* プランの削除 */
        fftw_destroy_plan(p);
        
        /* 結果をコピー */
        for ( j = 0; j < nX; j++ )
        {
            a_re[i*nX+j] = out[j][0] / nX;            /* 正規化が必要 */
            a_im[i*nX+j] = out[j][1] / nX;
        }
    }
    
    /* FFTW用領域の開放 */
    fftw_free(in); fftw_free(out);
    
    /* 作業領域確保 */
    w_re = new double[nY*nX];
    w_im = new double[nY*nX];

    /* 実部、虚部両方のデータを転置 */
    Transverse_matrix( a_re, w_re, nX, nY );
    Transverse_matrix( a_im, w_im, nX, nY );
    
    /* 再び領域確保、領域サイズに注意 */
    in = fftw_malloc( nY * sizeof(fftw_complex) );
    out = fftw_malloc( nY * sizeof(fftw_complex) );
    
    /* 縦方向のFFT */
    for ( i = 0; i < nX; i++ )
    {
        for ( j = 0; j < nY; j++ )
        {
            in[j][0] = w_re[i*nY+j];
            in[j][1] = w_im[i*nY+j];
        }
        
        /* プランの作成 */
        p = fftw_plan_dft_1d(nY, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
        
        /* FFT実行 */
        fftw_execute(p);
        
        /* プランの削除 */
        fftw_destroy_plan(p);
        
        /* 結果をコピー */
        for ( j = 0; j < nY; j++ )
        {
            w_re[i*nY+j] = out[j][0] / nY;            /* 正規化が必要 */
            w_im[i*nY+j] = out[j][1] / nY;
        }
    }

    fftw_free(in); fftw_free(out);

    /* 転置して元に戻した後、結果をコピー */
    Transverse_matrix_R( w_re, a_re, nY, nX );
    Transverse_matrix_R( w_im, a_im, nY, nX );
    
    delete w_re;
    delete w_im;
}

参考文献