// =================
// || Worksheet 11
// || In this exercise, we are going be solving the 2d wave equation using fourier derivatives.
//    d^2/dt^2 f(x,y,t) = - k (d^2/dx^2+d^2/dy^2) f(x,y,t)
// =================

// see also, e.g. http://www.mtnmath.com/whatrh/node66.html

#include "program.h"

using namespace std;
using namespace cv;

int worksheetthirteen(Mat& src)
{
    int o_sz = getOptimalDFTSize(200);
    Mat input = Mat::zeros(o_sz, o_sz, CV_32F);
    
    // draw a little circle in the middle of our image.
    circle(input, Point(input.rows/2, input.cols/2), 25.0f, Scalar( 1.0, 1.0, 1.0 ),CV_FILLED);
    
    // First, lets embed our source image in a complex image, so that we can easily take its fourier transform.
    Mat complexinput = embedInComplexImage(src);
    
    // a window to show our outputs
    cv::namedWindow( "current frame", WINDOW_AUTOSIZE );
    
    int MAX_NUMBER_OF_FRAMES = 5000;
    
    // dissipation factor.
    float k = 0.01f;
    
    // this is of the initial data
    Mat prev_frame = input.clone() ;
    Mat prev_prev_frame = input.clone();
    
    Mat colorized_frame;
    
    // repeat this for each of the frames.
    for (int n = 0; n < MAX_NUMBER_OF_FRAMES; n++)
    {
        // embed the previous image in a complex image.
        Mat cmplx_prev = embedInComplexImage(prev_frame);
        // and now take the DFT.
        dft(cmplx_prev, cmplx_prev);
        
        // Now, to compute the derivative of an image in fourier space, we use the formula
        // d/dx^2 sin( w x) = - w^2 sin (w x)
        // i.e. when working with sine waves, taking two derivatives is the same
        // as multiplying by the negative of the frequency squared (-w^2)  !!!
        // Now, in fourier space, our frequency^2 is the same as the distance^2 from the corners.
        // but if we re arrange our image...
        rearrange(cmplx_prev);
        // then our frequency squared is the same as the distance^2 from the center
        
        // so. loop over all the pixels in the complex_image.
        
        for(int i =0; i < o_sz; i ++)
        {
            for( int j = 0; j < o_sz; j++)
            {
                // multiply the pixel at (i,j) by its distance from the center pixel.
                // the center pixel is (o_sz/2, o_sz/2)
                
                // cmplx_prev.at<complex<float>>(i,j) *= ????
                
                // note we are using the operator *= here
                // this means 'times equals'
                // myVariable *= 5;
                // is exactly the same as:
                // myVariable = 5* myVariable;
            }
            
        }
        
        // rearrange our image back
        rearrange(cmplx_prev);
        
        // perform the IDFT of our complex image.
        cv::Mat Lframe; // this is the laplacian
        cv::dft(cmplx_prev, Lframe, cv::DFT_INVERSE | cv::DFT_REAL_OUTPUT);
        
        // the image Lframe is now the Laplacian of our image!
        
        // f(t+1) =   -f(t-1)        +   2 * f(t)      - k * L f(t)
        Mat frame = -prev_prev_frame + 2.0f*prev_frame - k * Lframe;
        
        // lets colorize our image too, so that it looks better.
        frame.convertTo(colorized_frame,CV_8UC1, 255);
        applyColorMap(colorized_frame, colorized_frame,COLORMAP_AUTUMN);
        // lets also make it 4 times the size, so we dont strain our eyes.
        resize(colorized_frame,colorized_frame,Size(colorized_frame.cols*4,colorized_frame.rows*4 ) );
        
        // display the current frame.
        cv::imshow( "current frame", colorized_frame);
        
        //wait for 'esc' key press for 300ms. If 'esc' key is pressed, break loop
        if (waitKey(10) == 27)
        {
            cout << "esc key is pressed by user" << endl;
            break;
        }
        
        // need to update PREVIOUS two frames. !!!!!
        prev_prev_frame = prev_frame.clone();
        prev_frame = frame.clone();
    }
    
    return 0;
}