#ifndef _FINITE_DIFFERENCE_METHOD_H_
#define _FINITE_DIFFERENCE_METHOD_H_

#include <opencv2\opencv.hpp>  
#include <opencv2\core\core.hpp>
#include <opencv2\core\mat.hpp>
#include <iostream>  
using namespace std;
using namespace cv;

#ifndef M_PI
#define M_PI 3.1415926535
#endif

/* 
 * Central Finite Difference Operators 
 * 
 * using OpenCV matrix data structure
 */
class CFiniteDifferenceOperator {

public:
    //constructor
    CFiniteDifferenceOperator() {}
    //destructor
    ~CFiniteDifferenceOperator() {};

    /*
     *  set the ghost cells to satisfy the 
     *  Neumann boundary condition
     *  u - the input function
     *  U - the extension of u with ghost cells
     */
    void set_ghost_cells(cv::Mat & U, const cv::Mat & u) {

        int width  = U.cols;
        int height = U.rows;

        assert(u.cols == width - 2 && u.rows == height - 2);

        /*
        *   set the central regular cells
        */
        for (int i = 1; i < width - 1; i++)
            for (int j = 1; j < height - 1; j++)
                U.at<double>(i, j) = u.at<double>(i - 1, j - 1);

        /* set the ghost cells on the boundary edges
        *  by copying the function value
        *  from the closest regular cells
        */
        for (int j = 1; j < height - 1; j++)
            U.at<double>(0, j) = U.at<double>(1, j);

        for (int j = 1; j < height - 1; j++)
            U.at<double>(width - 1, j) = U.at<double>(width - 2, j);

        for (int i = 1; i < width - 1; i++)
            U.at<double>(i, 0) = U.at<double>(i, 1);

        for (int i = 1; i < width - 1; i++)
            U.at<double>(i, height - 1) = U.at<double>(i, height - 2);

        /* set the ghost cells at the corners
        *  by copying the function value
        *  from the closest regular cells
        */
        U.at<double>(0, 0) = U.at<double>(1, 1);
        U.at<double>(width - 1, height - 1) = U.at<double>(width - 2, height - 2);
        U.at<double>(width - 1, 0)  = U.at<double>(width - 2, 1);
        U.at<double>(0, height - 1) = U.at<double>(1, height - 2);
    }

    /*
     *  Compute Dxx(U), namely u_{xx} using 
     *  the central difference operator
     * 
     *  U   - the input matrix with ghost cells
     *  Dxx - the output u_{xx}
     *  dx  - the horizontal spacial step length
     */
    void _dxx(const cv::Mat& U, cv::Mat& Dxx, const double hx)
    {
        int width = U.cols;
        int height = U.rows;

        assert(Dxx.width == width - 2 && Dxx.height == height - 2);
        
        for (int i = 1; i < width - 1; i++) {
            for (int j = 1; j < height - 1; j++) {
                //central difference operator
                Dxx.at<double>(i - 1, j - 1) = (U.at<double>(i + 1, j) + U.at<double>(i - 1, j) - 2.0 * U.at<double>(i, j)) / (hx * hx);
            }
        }
    }

    /*
     *  Compute Dxy(U), namely u_{yy} using 
     *  the central difference operator
     * 
     *  U   - the input matrix with ghost cells
     *  Dyy - the output u_{yy}
     *  dy  - the vertical spacial step length
     */
    void _dyy(const cv::Mat& U, cv::Mat& Dyy, const double hy)
    {
        int width = U.cols;
        int height = U.rows;

        assert(Dyy.width == width - 2 && Dyy.height == height - 2);

        for (int i = 1; i < width - 1; i++) {
            for (int j = 1; j < height - 1; j++) {
                //central difference operator
                Dyy.at<double>(i - 1, j - 1) = (U.at<double>(i, j + 1) + U.at<double>(i, j - 1) - 2.0 * U.at<double>(i, j)) / (hy * hy);
            }
        }
    }

    /*
     *  Compute Dxy(U), namely u_{xy} using
     *  the central difference operator
     *
     *  U   - the input matrix with ghost cells
     *  Dxy - the output u_{xy}
     *  dx  - the horizontal spacial step length
     *  dy  - the vertical spacial step length
     */
    void _dxy(const cv::Mat& U, cv::Mat& Dxy, const double hx, const double hy) {

        int width = U.cols;
        int height = U.rows;

        assert(Dxy.width == width - 2 && Dxy.height == height - 2);

        for (int i = 1; i < width - 1; i++) {
            for (int j = 1; j < height - 1; j++) {
                //central difference operator
                Dxy.at<double>(i - 1, j - 1) = (
                    U.at<double>(i + 1, j + 1) +
                    U.at<double>(i - 1, j - 1) -
                    U.at<double>(i - 1, j + 1) -
                    U.at<double>(i + 1, j - 1)) / (4 * hx * hy);
            }
        }
    }

    /*
     *  Compute Dx(U), namely u_{x} using
     *  the central difference operator
     *
     *  U   - the input matrix with ghost cells
     *  Dx - the output u_{x}
     *  dx  - the horizontal spacial step length
     */
    void _dx(cv::Mat & U, cv::Mat & Dx, double hx)
    {
        int width = U.cols;
        int height = U.rows;

        assert(Dx.cols == width - 2 && Dx.rows == height - 2);

        for (int i = 1; i < width - 1; i++) {
            for (int j = 1; j < height - 1; j++) {
                //central difference operator
                Dx.at<double>(i - 1, j - 1) = (U.at<double>(i + 1, j) - U.at<double>(i-1, j)) / (2*hx);
                //Dx.at<double>(i - 1,j - 1) = (U.at<double>(i + 1,j) - U.at<double>(i,j)) / hx;
            }
        }
    }

    /*
     *  Compute Dx(U), namely u_{y} using
     *  the central difference operator
     *
     *  U   - the input matrix with ghost cells
     *  Dx - the output u_{y}
     *  dy  - the vertical spacial step length
     */
    void _dy(cv::Mat & U, cv::Mat & Dy, double hy)
    {
        int width = U.cols;
        int height = U.rows;

        assert(Dy.cols == width - 2 && Dy.rows == height - 2);

        for (int i = 1; i < width - 1; i++) {
            for (int j = 1; j < height - 1; j++) {
                //central difference operator
                Dy.at<double>(i - 1,j - 1) = (U.at<double>(i,j + 1) - U.at<double>(i,j - 1)) / (2*hy);
            }
        }
    }
    /* convert matrix index (i,j) to the planar coordinates (x,y)
    *  [0,n1-1] X [0,n2-1] corresponds to [-1,1] X [-1, 1]
     * (0,0) corresponds to (-1 + h1/2, -1 + h2/2) 
     * (i,j)   - matrix element indices
     * (n1,n2) - dimension of the matrix
     * output is the correponding (x,y)
     */
    cv::Vec2d coord(const int i, const int j, const int n1, const int n2)
    {
        double h1 = 2.0 / n1;
        double h2 = 2.0 / n2;

        double x = -1.0 + h1 / 2 + i * h1;
        double y = -1.0 + h2 / 2 + j * h2;

        return cv::Vec2d(x, y);
    }

    /*
     *  compute the square 
     */
    double _sqr(const double x) {
        return x * x;
    }

    /*
     *  visualize the matrix by showing the normalized
     *  matrix as an image
     */
    void visualize(const std::string  title, const Mat& u)
    {
        Mat nu;
        u.copyTo(nu);
        
        int width  = u.cols;
        int height = u.rows;
        double min, max;
        //find the minimal and maximal values of the matrix
        minMaxIdx(nu, &min, &max);
        //normalize the matrix
        for (int i = 0; i < width; i++)
            for (int j = 0; j < height; j++) {
                double& v = nu.at<double>(i, j);
                v = (v - min) / (max - min);
            }
        //display the image
        cv::imshow(title, nu);
        std::cout << "Press any key to continue " << std::endl;
        cv::waitKey();
    }
};


#endif
