1. 程式人生 > >FFTW3庫對影象進行FFT的例項

FFTW3庫對影象進行FFT的例項

正常使用fftw_plan_dft_2d對影象進行FFT

要使用fftw_plan_dft_2d函式首先得下載並安裝官方的庫,安裝的時候要注意64位系統與32位系統的不同的執行步驟,我安裝的時候弄了挺久的

程式碼塊

這裡以圖pRefImg為例

    int nFFTRow=0;
    int nFFTCol=0;

    nFFTRow=nRefRow+nMaskRow-1;
    nFFTCol=nRefCol+nMaskCol-1;

    fftw_complex *pRefFFTIn= (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * nFFTRow*nFFTCol);//定義FFT的輸入
fftw_complex *pRefFFTOut= (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * nFFTRow*nFFTCol);//定義FFT的輸出 fftw_plan p_for_ref;//生成plan p_for_ref = fftw_plan_dft_2d(nFFTRow, nFFTCol, pRefFFTIn, pRefFFTOut, FFTW_FORWARD, FFTW_ESTIMATE);//只要把FFTW_FORWARD改為FFTW_BACKWARD就是IFFT了 for (i=0; i<nFFTRow; i++)//對拓展的影象進行補零
{ for (j=0; j<nFFTCol; j++) { if ((i < nRefRow) && (j < nRefCol)) { pRefFFTIn[i*nFFTCol+j][0]= (double)pRefImg[i][j];//實部,pRefImg為參考圖片 pRefFFTIn[i*nFFTCol+j][1] = 0;//虛部 } else { pRefFFTIn[i*nFFTCol+j][0
] = 0; pRefFFTIn[i*nFFTCol+j][1] = 0; } } } fftw_execute(p_for_ref); //執行plan fftw_destroy_plan(p_for_ref); if(pRefFFTIn!=NULL) fftw_free(pRefFFTIn);//釋放記憶體

執行此段程式碼後陣列pRefFFTOut就存放了變換後的結果,pRefFFTOut[…][0]儲存的為實部、pRefFFTOut[…][1]儲存的為虛部。

可以注意到,2dfft的輸入輸出都是兩個很長很長的陣列,實際上在函式的內部也是進行的1dfft,我們告訴函式有多少行多少列,目的就是讓函式知道要進行多少點的fft進行多少次。

使用fftw_plan_dft_1d函式對影象進行FFT

如何由1DFFT對影象進行2D的FFT呢,很簡單,只需要對影象的每一行進行FFT,然後轉置,然後再對每一行進行FFT,然後再轉置就行了

程式碼塊

以下為擷取的一部分


    fftw_complex *pRefFFTIn = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nFFTRow*nFFTCol);
    for (i=0; i<nFFTRow; i++)//拓展補零
    {
        for (j=0; j<nFFTCol; j++) 
        {
            if ((i < nRefRow) && (j < nRefCol))
            {
                pRefFFTIn[i*nFFTCol+j][0] = (double)pRefImg[i][j];//輸入影象
                pRefFFTIn[i*nFFTCol+j][1] = 0;
            }
            else
            {
                pRefFFTIn[i*nFFTCol+j][0] = 0;
                pRefFFTIn[i*nFFTCol+j][1] = 0;
            }
        }   
    }

    double **pRefFFTResult1_re = (double**)malloc(sizeof(double*)*nFFTRow);//開闢儲存陣列 實部
    double **pRefFFTResult1_im = (double**)malloc(sizeof(double*)*nFFTRow);//開闢儲存陣列 虛部
    for(int i=0;i<nFFTRow;i++)                     
    {
        pRefFFTResult1_re[i] = (double*)malloc(sizeof(double)*nFFTCol);  
        pRefFFTResult1_im[i] = (double*)malloc(sizeof(double)*nFFTCol);
    } 
    //第一次fft
    for (i=0; i<nFFTRow; i++)
    {
        fftw_complex *pRefFFTRowIn;
        fftw_complex *pRefFFTRowOut;
        pRefFFTRowIn  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nFFTCol);
        pRefFFTRowOut = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nFFTCol);
        fftw_plan p_for_ref1;
        p_for_ref1 = fftw_plan_dft_1d(nFFTCol, pRefFFTRowIn, pRefFFTRowOut, FFTW_FORWARD,FFTW_ESTIMATE);

        for (j=0;j<nFFTCol;j++)//行提取
        {
            pRefFFTRowIn[j][0] = pRefFFTIn[i*nFFTCol+j][0];
            pRefFFTRowIn[j][1] = pRefFFTIn[i*nFFTCol+j][1];
        }

        fftw_execute(p_for_ref1); 

        for (j=0;j<nFFTCol;j++)
        {
            pRefFFTResult1_re[i][j] = pRefFFTRowOut[j][0];
            pRefFFTResult1_im[i][j] = pRefFFTRowOut[j][1];
        }
        fftw_destroy_plan(p_for_ref1);
        //釋放記憶體
        if(pRefFFTRowIn!=NULL) fftw_free(pRefFFTRowIn);
        if(pRefFFTRowOut!=NULL) fftw_free(pRefFFTRowOut);
    }

    //轉置
    double **pRefFFTResult1T_re = (double**)malloc(sizeof(double*)*nFFTRow);
    double **pRefFFTResult1T_im = (double**)malloc(sizeof(double*)*nFFTRow);
    for(int i=0;i<nFFTRow;i++)                     
    {
        pRefFFTResult1T_re[i] = (double*)malloc(sizeof(double)*nFFTCol);  
        pRefFFTResult1T_im[i] = (double*)malloc(sizeof(double)*nFFTCol);
    } 
    for (i=0; i<nFFTRow; i++)
    {
        for (j=0; j<nFFTCol;j++)
        {
            pRefFFTResult1T_re[i][j] = pRefFFTResult1_re[j][i];
            pRefFFTResult1T_im[i][j] = pRefFFTResult1_im[j][i];
        }
    }

    //第二次fft
    double **pRefFFTResult_re = (double**)malloc(sizeof(double*)*nFFTRow);
    double **pRefFFTResult_im = (double**)malloc(sizeof(double*)*nFFTRow);
    for(int i=0;i<nFFTRow;i++)                 
    {
        pRefFFTResult_re[i] = (double*)malloc(sizeof(double)*nFFTCol);  
        pRefFFTResult_im[i] = (double*)malloc(sizeof(double)*nFFTCol);
    } 

    for (i=0; i<nFFTRow; i++)
    {

        fftw_complex *pRefFFTColIn ;
        fftw_complex *pRefFFTColOut ;
        pRefFFTColIn  = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * nFFTRow);
        pRefFFTColOut = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * nFFTRow);
        fftw_plan p_for_ref2;
        p_for_ref2 = fftw_plan_dft_1d(nFFTRow, pRefFFTColIn, pRefFFTColOut, FFTW_FORWARD,FFTW_ESTIMATE);

        for(j=0;j<nFFTCol;j++)//行提取
        {
            pRefFFTColIn[j][0] = pRefFFTResult1T_re[i][j];
            pRefFFTColIn[j][1] = pRefFFTResult1T_im[i][j];
        }

        fftw_execute(p_for_ref2); 

        for(j=0;j<nFFTCol;j++)
        {
            pRefFFTResult_re[i][j] = pRefFFTColOut[j][0];
            pRefFFTResult_im[i][j] = pRefFFTColOut[j][1];
        }
        fftw_destroy_plan(p_for_ref2);
        if(pRefFFTColIn!=NULL) fftw_free(pRefFFTColIn);
        if(pRefFFTColOut!=NULL) fftw_free(pRefFFTColOut);
    }

    //轉置
    double **pRefFFTResultT_re = (double**)malloc(sizeof(double*)*nFFTRow);
    double **pRefFFTResultT_im = (double**)malloc(sizeof(double*)*nFFTRow);
    for(int i=0;i<nFFTRow;i++)                 
    {
        pRefFFTResultT_re[i] = (double*)malloc(sizeof(double)*nFFTCol);
        pRefFFTResultT_im[i] = (double*)malloc(sizeof(double)*nFFTCol);
    } 
    for (i=0; i<nFFTRow; i++)
    {
        for (j=0; j<nFFTCol;j++)
        {
            pRefFFTResultT_re[i][j] = pRefFFTResult_re[j][i];//最終結果的實部
            pRefFFTResultT_im[i][j] = pRefFFTResult_im[j][i];//最終結果的虛部
        }
    }

程式解讀: 1.首先將需要進行變換的圖片根據輸入一幅圖片還是兩幅圖片進行按需要補零; 2.定義陣列來儲存中間結果並且做轉置變換,注意此時需要自己定義double型的二維陣列,不能直接用庫裡的fftw_double來定義陣列; 3.進行第一次fft 4.進行儲存並轉置 5.進行第二次fft 6.進行儲存並轉置,得最終結果 有一個問題,以上兩種方法轉換後進行模板匹配都是成功得,但是比較輸出的資料,兩種方法得到的結果並不是一樣的,這一點有待討論