Digital Image Processing#3 - 離散傅立葉轉換(Discrete Fourier Transform) 與 Gaussian Filter

文章流程：

第一次傅立葉就上手
一行心得

第一次傅立葉就上手，實作流程詳述！

以下為助教範例提供的ex2.bmp，解析度64 x 64：

===這是原圖，雖然沒特色，但方便Debug===

第一步，前餐，在傅立葉之前

//前置準備工作, 讀取圖片 與 定義參數, 以下數值針對範例舉例
Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE); // 讀取灰階圖
print( I.type() ); // 確認型態為 CV_8UC1, uchar 8 / channel 1
int M = I.cols; // M(width), 對傅立葉loop重要, 值為64
int N = I.rows; // N(height), 值為64

...

// zero-padding, I -> padded, 變兩倍大小, 為了不混疊
Mat padded;
copyMakeBorder(I, padded, 0, N, 0, M, BORDER_CONSTANT, Scalar::all(0)); 
Mat plane = Mat_(padded); // 為了運算方便, 轉成double, 配合之後的CV_64F型態的Mat
int P = plane.cols; // P = M / 2, 對傅立葉loop重要, 值為128, 
int Q = plane.rows; // Q = N / 2, 值為128, 

// 奇數座標, 值乘以-1, 將頻域shift to center
// 參考 : Spatial domain*-1^(x+y) <=> Frequency domain F(u-P/2, v-Q/2)
for (int y = 0; y < N; ++y)  
　　　for (int x = 0; x < M; ++x)  
　　　　　　if( (x+y)&1 )
　　　　　　　　　plane.ptr(y)[x] = plane.ptr(y)[x] * -1.0;

第二步，魔王的離散傅立葉 D F T

    // 參數快取, 加速計算
    double PI_2 = 2.0 * CV_PI;
    double angleTerm = 0.0;
    double invP = 1.0 / P; // 注意
    double invQ = 1.0 / Q; // 注意

    // 實數. 虛數. Spectrum, Phase一次搞定, CV_64F型態
    Mat Real = Mat::zeros(plane.size(), CV_64F); // 注意, size = P*Q
    Mat Imag = Mat::zeros(plane.size(), CV_64F); // P*Q
    Mat Spectrum = Mat::zeros(plane.size(), CV_64F); 
    Mat Phase = Mat::zeros(plane.size(), CV_64F);
    
    // 傅立葉, 外圈循環計算F(u,v), 內圈循環f(x,y), 注意loop範圍, P.Q.M.N看清楚
    for (int v = 0; v < Q; ++v)  
    {  
        double v_invQ = (double)(v)*invQ; // 加速計算
        for (int u = 0; u < P; ++u)  
        {
            double u_invP = (double)(u)*invP;
            double rd = 0.0; // 實部加總
            double id = 0.0; // 虛部加總

            for (int y = 0; y < N; ++y)  
            {  
                double vy_invQ = (double)(y)*v_invQ;
                for (int x = 0; x < M; ++x)  
                {  
                    // 參照傅立葉. 尤拉公式, 參數聯想Cos與Sin函數的震幅.頻率.位移.週期
                    // e^(j*angleTerm) = cos(angleTerm) + j*sin(angleTerm)
                    angleTerm = PI_2*( (double)(x)*u_invP + vy_invQ ); 
                    double target = plane.ptr(y)[x]; // f(x,y)

                    rd += (double)(target*cos(angleTerm));

                    id -= (double)(target*sin(angleTerm));
                }
            }  
            Real.ptr(v)[u] = rd;
            Imag.ptr(v)[u] = id;
            
            // 參見公式
            Spectrum.ptr(v)[u] = (double)(sqrt(rd*rd + id*id)); 
            Phase.ptr(v)[u] = (double)(atan2(id, rd));
        }  
    }

DFT卡關很久，原因在於公式上的M. N 與 P. Q搞混，反覆嘗試後，才得出正確結果

轉換到，順便得出兩個目標，Spectrum與Phase angle，顯示出來！

//imshow可接受範圍為[0-255]或[0-1], 兩個Mat型態都是double, 選擇正規化[0,1]後顯示
normalize(Spectrum, Spectrum, 0, 1, CV_MINMAX);
normalize(Phase, Phase, 0, 1, CV_MINMAX);
// viewable image form (float between values 0 and 1).

// Show the result
imshow("Input Image", I);    // Show the result
imshow("Spectrum Image", Spectrum);    // Show the result
imshow("Phase angle Image", Phase);

===是狙擊手，快趴下！ Spectrum (沒做Log)===

===質 Phase angle 感===

第2.5步，如果顯示結果非預期...

提供點Debug技巧，拆解傅立葉公式，從觀念思考起：

===傅立葉公式，關鍵的angleTerm===

===尤拉公式，將 e 變成我們熟悉的cos與sin===

===根據cos與sin函數的特質，推出可能的問題點，嘗試看看===

第三步，EX關卡反離散傅立葉I D F T

想說解決傅立葉，剩下都小咖了，沒想到又被 IDFT 的小細節卡到，挑出細節講述

    //double invP = 1.0 / P; // 注意, 仍是P與Q
    //double invQ = 1.0 / Q; // 注意
    double invPQ = invP * invQ;
    double cosineA, sineA;
    Mat Recover = Mat::zeros(I.size(), CV_64F); // Size=M*N

    for (int y = 0; y < N; ++y)  
    {  
        double y_invQ = (double)(y)*invQ;
        for (int x = 0; x < M; ++x)  
        { 
            double x_invP = (double)(x)*invP;
            double inv_rd = 0.0; // 加總實部
            double inv_id = 0.0; // 加總虛部

            for (int v = 0; v < Q; ++v)  
            {  
                double vy_invQ = (double)(v)*y_invQ;
                for (int u = 0; u < P; ++u)  
                {
                    angleTerm = PI_2*( (double)(u)*x_invP + vy_invQ );
                    cosineA = cos(angleTerm);
                    sineA = sin(angleTerm);
                    double rd = Real.ptr(v)[u]; // Re(F(u,v))
                    double id = Imag.ptr(v)[u]; // Im(F(u,v))

                    // 實部係數為( real * cos x - imag * sin x )
                    // 虛數係數為( real * sin x + imag * cos x )
                    inv_rd += (rd*cosineA - id*sineA); // 實*實 - 虛*虛(變號)
                    inv_id += (rd*sineA + id*cosineA); // 實*虛 + 虛*實
                }
            }
            inv_rd *= invPQ; // 對應公式上的 1/MN
            inv_id *= invPQ;
            //double odd_mul = ( (x+y)&1 ) ? -1 : 1; // 不知為何不用乘-1, 待釐清
            // 實部平方 + 虛部平方, 再開根號
            Recover.ptr(y)[x] = (double)(sqrt(inv_rd*inv_rd + inv_id*inv_id));
        }
    }

   // 第二個顯示技巧, 用convertTo函數, 將CV_64F轉成CV_8U
   Mat rec_out;
   Recover.convertTo(rec_out, CV_8U);
   imshow("Recover", rec_out);

　　// 小技巧, Image Compare, 確認產生結果與原圖一致, errorL2與similarity越小越相似
   double errorL2 = norm( I, rec_out, CV_L2 );
   double similarity = errorL2 / (double)( I.rows * I.cols );

===IDFT還原，一模一樣就對了===

第四步，Gaussian Filter

著重介紹頻域的，空間域的則找維基，公式帶入

  // Gaussian Frequency Filter, 此段放入DFT後(轉完頻域), IDFT前
  double pow_D02 = 2.0*D0*D0; // D0 可以給使用者輸入
  for (int v = 0; v < Q; ++v)
  {  
     for (int u = 0; u <; P; ++u)
     {
        // M = P/2, N = Q/2, 公式一樣, 自行代換
        double Duv = (double)(sqrt( (u-M)*(u-M)+(v-N)*(v-N) ));
        double Huv = exp(-(Duv*Duv/pow_D02)); // 注意, 不是尤拉, 而是exp

        Real.ptr(v)[u] *= Huv;
        Imag.ptr(v)[u] *= Huv;
     }
  }

===Gaussian Frequency Filter，D0=10.0===

===Gaussian Spatial Filter 3x3===

===Gaussian Spatial Filter 5x5===

一行心得

資質不好，耗費三天，寫得好崩潰XD，收穫挺多的，有種把上課學到的，應用起來的FU！

991412 影像處理學習部落格

2014年5月12日星期一