#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

#include <stdio.h>

using namespace cv;
using namespace std;

static void help()
{
    printf("\nThis program demonstrated the use of the discrete Fourier transform (dft)\n"
           "The dft of an image is taken and it's power spectrum is displayed.\n"
           "Usage:\n"
            "./dft [image_name -- default lena.jpg)]\n");
}

const char* keys =
{
    "{1| |lena.jpg)|input image file}"
};

int main(int argc, const char ** argv)
{
    help();
    CommandLineParser parser(argc, argv, keys);
    string filename = parser.get<string>("1");

    //画像をグレースケールで読み込む+エラー処理
    Mat img = imread(filename.c_str(), CV_LOAD_IMAGE_GRAYSCALE);
    if( img.empty() )
    {
        help();
        printf("Cannot read image file: %s\n", filename.c_str());
        return -1;
    }
    //DFTに最適なサイズを取得
    int M = getOptimalDFTSize( img.rows );
    int N = getOptimalDFTSize( img.cols );
    //最適なサイズはimgより大きいので周辺を0で埋めたMatを作成
    Mat padded;
    copyMakeBorder(img, padded, 0, M - img.rows, 0, N - img.cols, BORDER_CONSTANT, Scalar::all(0));
    //複素数画像(実部チャネルと虚部チャネルの2チャネル)の準備
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    //2枚のMatをmergeしてcomplexImgへ収める
    Mat complexImg;
    merge(planes, 2, complexImg);
    //DFT処理
    dft(complexImg, complexImg);

    // DFTの結果は値の幅が大きいのでlogにするが，値が0だとエラーになるので全体に1を加えてから
    // DFTの結果は√((実部)^2+(虚部)^2) とする
    // compute log(1 + sqrt(Re(DFT(img))**2 + Im(DFT(img))**2))
    split(complexImg, planes);
    //要素を2乗したものをplanes[0]に上書きで収める
    magnitude(planes[0], planes[1], planes[0]);
    Mat mag = planes[0];
    // 1足してlog
    mag += Scalar::all(1);
    log(mag, mag);

    // crop the spectrum, if it has an odd number of rows or columns
    // 大きさが奇数の場合削る
    mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));

    int cx = mag.cols/2;
    int cy = mag.rows/2;

    // rearrange the quadrants of Fourier image
    // so that the origin is at the image center
    //低周波が四隅になっているので入れ替えて中心を直流に
    // MATLABならfftshiftで済むが
    // OpenCVの場合はMatにRectを指定してROI(Region of Interest)
    // をq0-3の4つ作成し入れ替える
    Mat tmp;
    Mat q0(mag, Rect(0, 0, cx, cy));
    Mat q1(mag, Rect(cx, 0, cx, cy));
    Mat q2(mag, Rect(0, cy, cx, cy));
    Mat q3(mag, Rect(cx, cy, cx, cy));

    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

    //最小最大を0-1に割当てるよう正規化
    normalize(mag, mag, 0, 1, CV_MINMAX);

    imshow("spectrum magnitude", mag);
    waitKey();
    return 0;
}

    IplImage *myDFT(IplImage *src)
    {
      IplImage* fftImg=0;
      IplImage* realImg=0;
      IplImage* imageImg=0;
      IplImage* prealImg=0;
      IplImage* pimageImg=0;
      IplImage* modulo=0;
    
      //for only 1 channel image
      if (src->nChannels > 1)
        return NULL;
      //for FFT
      realImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
      imageImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
      fftImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,2);
      prealImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
      pimageImg=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
      modulo=cvCreateImage(cvGetSize(srcImg),IPL_DEPTH_64F,1);
    
      cvConvertScale(src,realImg,1,0); // uchar image to double image
      cvZero(imageImg);
      cvMerge(realImg,imageImg,NULL,NULL,fftImg); 
      cvDFT(fftImg,fftImg,CV_DXT_FORWARD);
      cvSplit(fftImg,realImg,imageImg,NULL,NULL);
      cvPow(realImg,prealImg,2); //prealImg=realImg^2
      cvPow(imageImg,pimageImg,2);//pimageImg=imageImg^2
      cvAdd(prealImg,pimageImg,modulo); //modulo = prealImg+pimageImg
      cvPow(modulo,modulo,0.5); //modulo = sqrt(modulo)
      
      //rearrange quadrant
      double tmp13,tmp24;
      int nRow = srcImg->height;
      int nCol = srcImg->width;
      int cy = nRow/2; // image center
      int cx = nCol/2;
      int i, j;
    
      for( j = 0; j < cy; j++ ){
        for( i = 0; i < cx; i++ ){
          tmp13 = CV_IMAGE_ELEM( modulo, double, j, i);
          CV_IMAGE_ELEM( modulo, double, j, i) = CV_IMAGE_ELEM(modulo, double, j+cy, i+cx);
          CV_IMAGE_ELEM( modulo, double, j+cy, i+cx) = tmp13;
    
          tmp24 = CV_IMAGE_ELEM( modulo, double, j, i+cx);
          CV_IMAGE_ELEM( modulo, double, j, i+cx) = CV_IMAGE_ELEM(modulo, double, j+cy, i);
          CV_IMAGE_ELEM( modulo, double, j+cy, i) = tmp24;
        }
      }
    
      //normalize
      CvPoint maxPoint,minPoint;
      double max,min;
      cvMinMaxLoc(modulo,&min,&max,&minPoint,&maxPoint);
      double scale=255.0/(max-min);
      double shift=min*scale;
      cvConvertScale(modulo,modulo,scale,shift);
    
      //clean up
      cvReleaseImage(&fftImg);
      cvReleaseImage(&realImg);
      cvReleaseImage(&imageImg);
      cvReleaseImage(&prealImg=0);
      cvReleaseImage(&pimageImg=0);
      return modulo;
    }

    #define IGNORE_FREQ 3
    
    void myDisplayDFT(int width, double *f64scanedLine)
    {
      //2チャンネルの1次元マトリクス
      CvMat *dftData = cvCreateMat(width, 1, CV_64FC2);
      //RealとIm用の1チャンネルの1次元マトリクス
      CvMat *dftReData = cvCreateMat(width, 1, CV_64FC1);
      CvMat *dftImData = cvCreateMat(width, 1, CV_64FC1);
      //表示用ダミーデータ
      //(この関数を繰り返し呼ぶ場合は下のmalloc,freeを繰り返すのは得策ではない。
      //mainなどの呼び出し元で確保してポインタを渡すのがスマート)
      double *f64DispDFTLine;
      //グラフ用画像
      //グラフは扱い易いように横軸をwidth[pixel]の幅にし、縦は255。上下左右のマージンは20とする。
      IplImage *graphImage=cvCreateImage(cvSize(width+40,296),IPL_DEPTH_8U, 3);
      int x;
      double max=0.0;
      f64DispDFTLine = (double*)malloc(sizeof(double)*width);
    
      //貰ったデータをdftDataの1チャネルに格納。アクセスの仕方はOpenCV wikiのFAQ参照
      //2チャネルは0で初期化(追加)
      for (x=0; x<width; x++){
        ((double*)(dftData->data.ptr+dftData->step*x))[0] = f64scanedLine[x];
        ((double*)(dftData->data.ptr+dftData->step*x))[1] = 0.0;
      }
      cvDFT(dftData,dftData,CV_DXT_FORWARD);
      //1,2チャンネルにそれぞれReal,Imパートのデータが入っているので分ける
      cvSplit(dftData,dftReData,dftImData,NULL,NULL);
      //二乗和のルート。マトリクスは使い回し。
      cvPow(dftReData,dftReData,2);
      cvPow(dftImData,dftImData,2);
      cvAdd(dftReData,dftImData,dftReData);
      cvPow(dftReData,dftReData,0.5);
      //マトリクスから扱い易いただの配列にデータを移す。同時に描画用に最大値を調べる
      for (x=0; x<width; x++){
        f64DFTLine[x]= *((double*)(dftReData->data.ptr+dftReData->step*x));
        if(x>=IGNORE_FREQ && x<width/2 && f64DFTLine[x]>max){
          max=f64DFTLine[x];
        }
      }
      //表示用配列に正規化した値を代入
      for (x=0; x<width; x++){
        f64DispDFTLine[x]=(255.0*f64DFTLine[x]/max)>255.0?255.0:(255.0*f64DFTLine[x]/max);
      }
      //グラフ作成
      cvSet(graphImage,CV_RGB(255,255,255)); //白塗り
      cvLine(graphImage,cvPoint(20,20),cvPoint(20,276),CV_RGB(0,0,0)); //軸
      cvLine(graphImage,cvPoint(20,276),cvPoint(width+20,276),CV_RGB(0,0,0)); //軸
      //スペクトルの範囲はwidth/2になるのでグラフを描画する際は横軸を倍にしている
      for (x=0; x<width/2; x++){  
        cvLine(graphImage,cvPoint(20+x*2,276-(int)f64DispDFTLine[x]),
               cvPoint(20+(x+1)*2,276-(int)f64DispDFTLine[x+1]),CV_RGB(255,0,0));
      }
      cvNamedWindow("Graph",1);
      cvShowImage("Graph",graphImage);
      cvWaitKey(0);
      cvSaveImage("graph.png)",graphImage);
      cvDestroyWindow("Graph");
      cvReleaseImage(&graphImage);
      free(f64DispDFTLine);
      cvReleaseMat(&dftData);
      cvReleaseMat(&dftReData);
      cvReleaseMat(&dftImData);
    }