//discrete fourier tranform, 離散傅立葉變換

//標頭檔案
#include "opencv2/core.hpp"     //Core functionality，核心函式相關
#include "opencv2/imgproc.hpp"  //Image processing， 影象處理相關
#include "opencv2/imgcodecs.hpp"//Image file reading and writing， 影象的載入和寫出相關
#include "opencv2/highgui.hpp"  //High-level GUI，圖形介面GUI相關

#include <iostream>
 


/**
 * 程式流程
 * 1、載入影象，格式為灰度圖
 * 2、獲取圖片dft變換的最佳大小
 * 3、邊框加0的方式填充圖片，即非0部分為dft變換的最佳大小 
 * 4、建立陣列儲存影象實部虛部，且合併到complexI
 * 5、傅立葉變換 dft(complexI, complexI)
 * 6、重新分離實部虛部，並且計算幅度
 * 7、將幅度對映到對數域
 * 8、以影象中心為原點劃分象限，每個象限建立一個ROI
 * 9、對角象限互換
 * 10、顯示結果
 */

//名稱空間
using namespace cv;
using namespace std;

//幫助函式，輸出程式的資訊
 

static void help(void)
{
    cout << endl
        <<  "This program demonstrated the use of the discrete Fourier transform (DFT). " << endl   //離散傅立葉變換示例
        <<  "The dft of an image is taken and it's power spectrum is displayed."          << endl   //離散傅立葉變換後顯示功率譜
        <<  "Usage:"
 
                                                                      << endl
        <<  "./discrete_fourier_transform [image_name -- default ../data/lena.jpg]"       << endl;  //預設載入圖片路徑
}

int main(int argc, char ** argv)
{
    help();

    //獲取影象路徑（檔名），命令列輸入否則預設
    const char* filename = argc >=2 ? argv[1] : "../data/lena.jpg";

    //載入影象，方式為載入灰度圖
    Mat I = imread(filename, IMREAD_GRAYSCALE);
    //檢查是否成功載入
    if( I.empty()){
        cout << "Error opening image" << endl;
        return -1;
    }

//! [expand]
    Mat padded;                            
    //expand input image to optimal size， 將輸入影象擴充套件到最佳大小
    int m = getOptimalDFTSize( I.rows );
    int n = getOptimalDFTSize( I.cols ); 

    // on the border add zero values，在邊框上新增零值，使用copyMakeBorder函式
    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));
//! [expand]

//! [complex_and_real] 實部和虛部
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; //Mat 陣列儲存影象的實部和虛部
    Mat complexI; 
    merge(planes, 2, complexI);         // Add to the expanded another plane with zeros，用零新增到擴充套件的另一平面

//! [complex_and_real]

//! [dft]
    //離散傅立葉變換
    dft(complexI, complexI);            // this way the result may fit in the source matrix，這種方式的結果可能適合在源矩陣

//! [dft]

    // compute the magnitude and switch to logarithmic scale，計算幅度並對映到對數刻度
    //公式 => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)) 
//! [magnitude] 幅度
    split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I)),分離實部和虛部
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude，計算幅度且存放到planes[0]
    Mat magI = planes[0]; //幅度 
//! [magnitude]

//! [log]
    magI += Scalar::all(1);                    // switch to logarithmic scale，對映到對數刻度
    log(magI, magI);
//! [log]

//! [crop_rearrange]裁剪重新排列
    // crop the spectrum, if it has an odd number of rows or columns， 裁剪頻譜, 如果它有奇數行或列數
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

    // rearrange the quadrants of Fourier image  so that the origin is at the image center
    //重新排列傅立葉影象的象限, 使原點位於影象中心
    int cx = magI.cols/2;
    int cy = magI.rows/2;

    //每個象限新建一個ROI
    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant， 左上，第二象限
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right， 右上，第一象限
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left，左下，第三象限
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right， 右下，第四象限

    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)，交換左上和右下象限
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)，交換右上和左下象限
    q2.copyTo(q1);
    tmp.copyTo(q2);
//! [crop_rearrange]

//! [normalize]
    //歸一化，畫素值都對映到[0,1]之間
    normalize(magI, magI, 0, 1, NORM_MINMAX); // Transform the matrix with float values into a
                                            // viewable image form (float between values 0 and 1).
//! [normalize]

    //顯示結果
    imshow("Input Image"       , I   );    // Show the result
    imshow("spectrum magnitude", magI);
    waitKey();

    return 0;
}

/**
 * 要點總結：
 * 載入圖片格式為灰度圖
 * getOptimalDFTSize()函式獲取最佳大小
 * copyMakeBorder()加框函式
 * 實部虛部
 * merge()合併函式
 * dft()函式
 * 幅度公式sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)
 * split()分離函式
 * magnitude()計算幅度
 * log()對數函式
 * normalize()歸一化函式
 * /