#ifndef DATAFUNCTION_HH
#define DATAFUNCTION_HH

#include <dune/acfem/functions/localfunctionwrapper.hh>
#include <dune/acfem/functions/gridfunctionexpression.hh>
#include <dune/acfem/models/modelinterface.hh>


//include CIMg library for reading in images
//documentation; http://cimg.eu/reference/
#define cimg_display 0
#include "CImg.h"



namespace Dune {
  namespace ACFem {


    //Construct Image gridFunction from some function

    template<class DiscreteFunctionType>
    class LocalDataAdapter
    {
    public:
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;

      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename GridPartType::template Codim<0>::EntityType EntityType;

      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;

      /**Construct the right hand side function from a given PDE model
       * and a given -- potentially non-discrete -- function.
       */
      LocalDataAdapter(const std::string& name = "")
        : name_(name)
      {}

      /**@internal
       * Needed, otherwise solution_ and localSolution_ may "fall apart"
       */
      LocalDataAdapter(const LocalDataAdapter& other)
        : name_(other.name_)
      {}

      const std::string& name() const { return name_; }

      // Implement needed interface method for the adapter class
      void init(const EntityType& entity)
      {
        //return circle
        const auto center =  entity.geometry().center();
        if( (center[0] - 0.5) * (center[0] -0.5 ) +(center[1] - 0.5) * (center[1] -0.5 ) <= 0.09) result_ = 1; //circle center 0.5,0.5 radius 0.3
        else
        /* a set of squares to check effect of alphas.
        //a describes the size difference
        double a=1;
        //r is the initial width
        const double r=0.25;
        //h is the height
        double h = 0.9;
        for (int i = 0; i <6 ; ++i)
        {
          if(center[0]<h && center[1] < h && center[0]>h-a*r && center[1] > h-a*r)
          {
            result_ =1;
            return;
          }
          h -= a*r+a*0.05;
          a *= 1./1.5;
        }
        */

        result_ = 0;
      }

      //! evaluate local function
      template<class PointType>
      void evaluate(const PointType& x, RangeType& result) const
      {
        //is constant on entities
        result = result_;
      }

      //! jacobian of local function
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        DUNE_THROW(NotImplemented, "Jacobian of an image just cannot be computed.");
      }

      // hessian of local function
      template<class PointType>
      void hessian(const PointType& x, HessianRangeType& ret) const
      {
        DUNE_THROW(NotImplemented, "Hessian of an image just cannot be computed.");
      }

    protected:
      mutable RangeType result_;
      const std::string name_;
    };


    template<class DiscreteFunctionType>
    LocalFunctionWrapper<LocalDataAdapter<DiscreteFunctionType>, typename DiscreteFunctionType::GridPartType>
    data(const typename DiscreteFunctionType::GridPartType & gridPart, const std::string& name = "")
    {
      typedef LocalDataAdapter<DiscreteFunctionType> LocalFunctionType;
      typedef LocalFunctionWrapper<LocalFunctionType, typename DiscreteFunctionType::GridPartType> GridFunctionType;

      LocalFunctionType L_f_dot(name);

      return GridFunctionType(L_f_dot.name(), L_f_dot, gridPart, 0);
    }


    //Read in Image as a gridFunction

    template<class DiscreteFunctionType>
    class LocalImageAdapter
    {
    public:
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;

      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename GridPartType::template Codim<0>::EntityType EntityType;

      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;

      typedef cimg_library::CImg<double> ImageType;

      /**Construct the right hand side function from a given PDE model
       * and a given -- potentially non-discrete -- function.
       */
      LocalImageAdapter(const std::string filename, const std::string& name = "")
        : image_(ImageType(filename.c_str()).normalize(0,1)), width_(image_.width()), height_(image_.height()), name_(name)
      {}

      /**@internal
       * Needed, otherwise solution_ and localSolution_ may "fall apart"
       */
      LocalImageAdapter(const LocalImageAdapter& other)
        : image_(other.image_), width_(other.width_), height_(other.height_), name_(other.name_)
      {}

      const std::string& name() const { return name_; }

      // Implement needed interface method for the adapter class
      void init(const EntityType& entity)
      {
        //actually read in image here
        DomainType center = entity.geometry().center();
        result_ = RangeType(*image_.data(std::round(center[0] * width_), std::round( (1 - center[1]) * height_) ) );
      }

      //! evaluate local function
      template<class PointType>
      void evaluate(const PointType& x, RangeType& result) const
      {
        //is constant on entities
        result = result_;
      }

      //! jacobian of local function
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        DUNE_THROW(NotImplemented, "Jacobian of an image just cannot be computed.");
      }

      // hessian of local function
      template<class PointType>
      void hessian(const PointType& x, HessianRangeType& ret) const
      {
        DUNE_THROW(NotImplemented, "Hessian of an image just cannot be computed.");
      }

    protected:
      mutable RangeType result_;
      const ImageType image_;
      const unsigned width_, height_;
      const std::string name_;
    };


    template<class DiscreteFunctionType>
    LocalFunctionWrapper<LocalImageAdapter<DiscreteFunctionType>, typename DiscreteFunctionType::GridPartType>
    image(const typename DiscreteFunctionType::GridPartType & gridPart, const std::string filename, const std::string& name = "")
    {
      typedef LocalImageAdapter<DiscreteFunctionType> LocalFunctionType;
      typedef LocalFunctionWrapper<LocalFunctionType, typename DiscreteFunctionType::GridPartType> GridFunctionType;

      LocalFunctionType L_f_dot(filename, name);

      return GridFunctionType(L_f_dot.name(), L_f_dot, gridPart, 0);
    }


    //Read in Image and noise it - export as gridFunction


    template<class DiscreteFunctionType>
    class LocalNoiseImageAdapter
    {
    public:
      typedef typename DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
      typedef typename DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType;

      typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
      typedef typename GridPartType::template Codim<0>::EntityType EntityType;

      typedef typename FunctionSpaceType::RangeType RangeType;
      typedef typename FunctionSpaceType::DomainType DomainType;
      typedef typename FunctionSpaceType::JacobianRangeType JacobianRangeType;
      typedef typename FunctionSpaceType::HessianRangeType HessianRangeType;

      typedef cimg_library::CImg<double> ImageType;

      /**Construct the right hand side function from a given PDE model
       * and a given -- potentially non-discrete -- function.
       */
      LocalNoiseImageAdapter(const std::string filename, const double noiseLevel1, const unsigned noiseType1 = 0, const double noiseLevel2 = 0, const unsigned noiseType2 = 0, const std::string& name = "")
        : image_(ImageType(filename.c_str()).normalize(0,1).noise(noiseLevel1, noiseType1).noise(noiseLevel2, noiseType2)), width_(image_.width()), height_(image_.height()), name_(name)
      {}

      /**@internal
       * Needed, otherwise solution_ and localSolution_ may "fall apart"
       */
      LocalNoiseImageAdapter(const LocalNoiseImageAdapter& other)
        : image_(other.image_), width_(other.width_), height_(other.height_), name_(other.name_)
      {}

      const std::string& name() const { return name_; }

      // Implement needed interface method for the adapter class
      void init(const EntityType& entity)
      {
        //actually read in image here
        DomainType center = entity.geometry().center();
        result_ = RangeType(*image_.data(std::round(center[0] * width_), std::round( (1 - center[1]) * height_) ) );
      }

      //! evaluate local function
      template<class PointType>
      void evaluate(const PointType& x, RangeType& result) const
      {
        //is constant on entities
        result = result_;
      }

      //! jacobian of local function
      template<class PointType>
      void jacobian(const PointType& x, JacobianRangeType& ret) const
      {
        DUNE_THROW(NotImplemented, "Jacobian of an image just cannot be computed.");
      }

      // hessian of local function
      template<class PointType>
      void hessian(const PointType& x, HessianRangeType& ret) const
      {
        DUNE_THROW(NotImplemented, "Hessian of an image just cannot be computed.");
      }

    protected:
      mutable RangeType result_;
      const ImageType image_;
      const unsigned width_, height_;
      const std::string name_;
    };


    template<class DiscreteFunctionType>
    LocalFunctionWrapper<LocalNoiseImageAdapter<DiscreteFunctionType>, typename DiscreteFunctionType::GridPartType>
    noiseImage(const typename DiscreteFunctionType::GridPartType & gridPart, const std::string filename, const double noiseLevel1, const unsigned noiseType1 = 0, const double noiseLevel2 = 0., const unsigned noiseType2 = 2, const std::string& name = "")
    {
      typedef LocalNoiseImageAdapter<DiscreteFunctionType> LocalFunctionType;
      typedef LocalFunctionWrapper<LocalFunctionType, typename DiscreteFunctionType::GridPartType> GridFunctionType;

      LocalFunctionType L_f_dot(filename, noiseLevel1, noiseType1, noiseLevel2, noiseType2, name);

      return GridFunctionType(L_f_dot.name(), L_f_dot, gridPart, 0);
    }


  } // ACFem::

} // Dune::


#endif // DATAFUNCTION_HH
