#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/highgui.hpp>
using namespace std;
/*
* dft函数的flags参数说明：
* DFT_INVERSE           反傅里叶变换
* DFT_SCALE与DFT_INVERSE结合使用, 输出结果除以元素数进行缩放
* DFT_ROWS 对输入的矩阵每一行进行傅里叶变换或反傅里叶变换, 同时对多个向量进行变换, 可以减小三维或多维变换操作的开销
* DFT_COMPLEX_OUTPUT 傅里叶变换的结果关于原点共轭对称，对称点的值实部相等虚部互为相反数。所以在一个象限中存储实部,
*                    对称象限存储虚部。这样, 两个通道存储的复数就只需要一个通道, 节省内存。默认输入矩阵为单通道,输出
*                    矩阵也为单通道(复数共轭对称压缩存储复数); 如果输入的矩阵双通道, 输出矩阵也为双通道; 如果指定此参数
*                    总是使用双通道输出结果
* DFT_REAL_OUTPUT    对共轭对称矩阵(如正向傅里叶变换生成的矩阵)进行反傅里叶变换时, 结果是一个实数矩阵。但函数不会判断共轭
*                    可以通过这个参数告诉函数输入的矩阵是共轭对称的, 从而输出实数矩阵。如果输入只有一个通道(复数共轭对称压缩存储复数)
*                    函数认为是一个经过压缩的共轭对称矩阵, 从而输出实数矩阵
* DFT_COMPLEX_INPUT 输入必须有两个通道分别表示实部和虚部, 如果输入有两个通道, 默认也是实部和虚部通道
*/
int main()
{
    cv::Mat gray = cv::imread("test.jpg", cv::IMREAD_GRAYSCALE);
    if (gray.empty())
    {
        cout << "failed to read image" << endl;
        return -1;
    }
    cv::namedWindow("src", cv::WINDOW_AUTOSIZE);
    cv::imshow("src", gray);

    cv::Mat matPadded;
    const int optimalRows = cv::getOptimalDFTSize(gray.rows);
    const int optimalCols = cv::getOptimalDFTSize(gray.cols);
    cv::copyMakeBorder(gray, matPadded, 0, optimalRows - gray.rows, 0, optimalCols - gray.cols, cv::BORDER_CONSTANT, cv::Scalar::all(0));

    cv::Mat matPaddedDouble;
    matPadded.convertTo(matPaddedDouble, CV_64FC1);

    cv::Mat dftResult;
    cv::dft(matPaddedDouble, dftResult, cv::DFT_COMPLEX_OUTPUT);  // need DFT_COMPLEX_OUTPUT, otherwise output is one channel
    if (dftResult.empty())
    {
        cout << "failed to dft" << endl;
        return -1;
    }

    cv::Mat dftResultChannels[] = { cv::Mat::zeros(dftResult.size(), CV_64FC1), cv::Mat::zeros(dftResult.size(), CV_64FC1) };
    cv::split(dftResult, dftResultChannels);                                  // split real and imag

    cv::Mat matMagnitude;                                                     // calculate magnitude
    cv::magnitude(dftResultChannels[0], dftResultChannels[1], matMagnitude);  // magnitude = sqrt(real ^ 2 + imag ^ 2);
    cv::imshow("raw magnitude", matMagnitude);                                // all white or all black
    matMagnitude += cv::Scalar::all(1);
    cv::log(matMagnitude, matMagnitude);
    cv::imshow("log magnitude", matMagnitude);
    cv::normalize(matMagnitude, matMagnitude, 0, 1, cv::NORM_MINMAX);
    cv::imshow("normalized log magnitude", matMagnitude);

    int rows = ((matMagnitude.rows & (2 - 1)) == 1) ? (matMagnitude.rows - 1) : matMagnitude.rows;
    int cols = ((matMagnitude.cols & (2 - 1)) == 1) ? (matMagnitude.cols - 1) : matMagnitude.cols;
    matMagnitude = matMagnitude(cv::Rect(0, 0, cols, rows));
    int centerX = cols >> 1, centerY = rows >> 1;
    cv::Mat topLeft(matMagnitude, cv::Rect(0, 0, centerX, centerY)); 
    cv::Mat topRight(matMagnitude, cv::Rect(centerX, 0, centerX, centerY));
    cv::Mat bottomLeft(matMagnitude, cv::Rect(0, centerY, centerX, centerY));
    cv::Mat bottomRight(matMagnitude, cv::Rect(centerX, centerY, centerX, centerY));
    cv::Mat tmp;
    topLeft.copyTo(tmp),  bottomRight.copyTo(topLeft), tmp.copyTo(bottomRight);
    topRight.copyTo(tmp), bottomLeft.copyTo(topRight), tmp.copyTo(bottomLeft);
    cv::imshow("centered magnitude", matMagnitude);
    cv::waitKey(0);
    return 0;
}