// @(#)root/mathcore:$Id: Delaunay2D.h,v 1.00
// Author: Daniel Funke, Lorenzo Moneta

/*************************************************************************
 * Copyright (C) 2015 ROOT Math Team                                     *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

// Header file for class Delaunay2D

#ifndef ROOT_Math_Delaunay2D
#define ROOT_Math_Delaunay2D

//for testing purposes HAS_CGAL can be [un]defined here
//#define HAS_CGAL

//for testing purposes THREAD_SAFE can [un]defined here
//#define THREAD_SAFE


//#include "RtypesCore.h"

#include <map>
#include <vector>
#include <set>
#include <functional>

#ifdef HAS_CGAL
   /* CGAL uses the name PTR as member name in its Handle class
    * but its a macro defined in mmalloc.h of ROOT
    * Safe it, disable it and then re-enable it later on*/
   #pragma push_macro("PTR")
   #undef PTR

   #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
   #include <CGAL/Delaunay_triangulation_2.h>
   #include <CGAL/Triangulation_vertex_base_with_info_2.h>
   #include <CGAL/Interpolation_traits_2.h>
   #include <CGAL/natural_neighbor_coordinates_2.h>
   #include <CGAL/interpolation_functions.h>

   #pragma pop_macro("PTR")
#endif

#ifdef THREAD_SAFE
   #include<atomic> //atomic operations for thread safety
#endif


namespace ROOT {



   namespace Math {

/**

   Class to generate a Delaunay triangulation of a 2D set of points.
   Algorithm based on **Triangle**, a two-dimensional quality mesh generator and
   Delaunay triangulator from Jonathan Richard Shewchuk.

   See [http://www.cs.cmu.edu/~quake/triangle.html]

   After having found the triangles using the above library,  barycentric coordinates are used
   to test whether a point is inside a triangle (inTriangle test) and for interpolation.
   All this below is implemented in the DoInterpolateNormalized function.

   Given triangle ABC and point P, P can be expressed by

     P.x = la * A.x + lb * B.x + lc * C.x
     P.y = la * A.y + lb * B.y + lc * C.y

   with lc = 1 - la - lb

     P.x = la * A.x + lb * B.x + (1-la-lb) * C.x
     P.y = la * A.y + lb * B.y + (1-la-lb) * C.y

   Rearranging yields

     la * (A.x - C.x) + lb * (B.x - C.x) = P.x - C.x
     la * (A.y - C.y) + lb * (B.y - C.y) = P.y - C.y

   Thus

     la = ( (B.y - C.y)*(P.x - C.x) + (C.x - B.x)*(P.y - C.y) ) / ( (B.y - C.y)*(A.x - C.x) + (C.x - B.x)*(A.y - C.y) )
     lb = ( (C.y - A.y)*(P.x - C.x) + (A.x - C.x)*(P.y - C.y) ) / ( (B.y - C.y)*(A.x - C.x) + (C.x - B.x)*(A.y - C.y) )
     lc = 1 - la - lb

   We save the inverse denominator to speedup computation

     invDenom = 1 / ( (B.y - C.y)*(A.x - C.x) + (C.x - B.x)*(A.y - C.y) )

     P is in triangle (including edges if

     0 <= [la, lb, lc] <= 1

     The interpolation of P.z is

     P.z = la * A.z + lb * B.z + lc * C.z

   To speed up localisation of points (to see to which triangle belong) a grid is laid over the internal coordinate space.
   A reference to triangle ABC is added to _all_ grid cells that include ABC's bounding box.
   The size of the grid is defined to be 25x25

   Optionally (if the compiler macro `HAS_GCAL` is defined ) the triangle findings and interpolation can be computed
   using the GCAL library. This is however not supported when using the class within ROOT

   \ingroup MathCore
 */


class Delaunay2D  {

public:

   struct Triangle {
      double x[3];           // x of triangle vertices
      double y[3];           // y of triangle vertices
      unsigned int idx[3];   // point corresponding to vertices

      #ifndef HAS_CGAL
      double invDenom; // cached inv denominator for computing barycentric coordinates (see above)
      #endif
   };

   typedef std::vector<Triangle> Triangles;

public:


   Delaunay2D(int n, const double *x, const double * y, const double * z, double xmin=0, double xmax=0, double ymin=0, double ymax=0);

   /// set the input points for building the graph
   void SetInputPoints(int n, const double *x, const double * y, const double * z, double xmin=0, double xmax=0, double ymin=0, double ymax=0);

   /// Return the Interpolated z value corresponding to the given (x,y) point.
   /// Note that in case no Delaunay triangles are found, for example when the
   /// points are aligned, then a default value of zero is always return.
   /// See the class documentation for  how the interpolation is computed.
   double  Interpolate(double x, double y);

   /// Find all triangles
   void      FindAllTriangles();

   /// return the number of triangles
   int    NumberOfTriangles() const {return fNdt;}

   double  XMin() const {return fXNmin;}
   double  XMax() const {return fXNmax;}
   double  YMin() const {return fYNmin;}
   double  YMax() const {return fYNmax;}

   /// set z value to be returned for points  outside the region
   void      SetZOuterValue(double z=0.) { fZout = z; }

   /// return the user defined Z-outer value
   double ZOuterValue() const { return fZout; }

   // iterators on the found triangles
   Triangles::const_iterator begin() const { return fTriangles.begin(); }
   Triangles::const_iterator end()  const { return fTriangles.end(); }


private:

   // internal methods


   inline double Linear_transform(double x, double offset, double factor){
      return (x+offset)*factor;
   }

   /// internal function to normalize the points.
   /// See the class documentation for the details on how it is computed.
   void DoNormalizePoints();

   /// internal function to find the triangle
   /// use Triangle or CGAL if flag is set
   void DoFindTriangles();

   /// internal method to compute the interpolation
   double  DoInterpolateNormalized(double x, double y);



private:
   // class is not copyable
   Delaunay2D(const Delaunay2D&); // Not implemented
   Delaunay2D& operator=(const Delaunay2D&); // Not implemented

protected:

   int         fNdt;           ///<! Number of Delaunay triangles found
   int         fNpoints;       ///<! Number of data points

   const double   *fX;         ///<! Pointer to X array (managed externally)
   const double   *fY;         ///<! Pointer to Y array
   const double   *fZ;         ///<! Pointer to Z array

   double    fXNmin;           ///<! Minimum value of fXN
   double    fXNmax;           ///<! Maximum value of fXN
   double    fYNmin;           ///<! Minimum value of fYN
   double    fYNmax;           ///<! Maximum value of fYN

   double    fOffsetX;         ///<! Normalization offset X
   double    fOffsetY;         ///<! Normalization offset Y

   double    fScaleFactorX;    ///<! Normalization factor X
   double    fScaleFactorY;    ///<! Normalization factor Y

   double    fZout;            ///<! Height for points lying outside the convex hull

   bool      fInit;            ///<! True if FindAllTriangles() has been performed


   Triangles   fTriangles;     ///<! Triangles of Triangulation

#ifndef HAS_CGAL

   //using triangle library

   std::vector<double> fXN; ///<! normalized X
   std::vector<double> fYN; ///<! normalized Y

   static const int fNCells = 25; ///<! number of cells to divide the normalized space
   double fXCellStep; ///<! inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
   double fYCellStep; ///<! inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
   std::set<unsigned int> fCells[(fNCells+1)*(fNCells+1)]; ///<! grid cells with containing triangles

   inline unsigned int Cell(unsigned int x, unsigned int y) const {
      return x*(fNCells+1) + y;
   }

   inline int CellX(double x) const {
      return (x - fXNmin) * fXCellStep;
   }

   inline int CellY(double y) const {
      return (y - fYNmin) * fYCellStep;
   }

#else // HAS_CGAL
      // case of using GCAL
      //Functor class for accessing the function values/gradients
      template< class PointWithInfoMap, typename ValueType >
      struct Data_access : public std::unary_function< typename PointWithInfoMap::key_type,
                std::pair<ValueType, bool> >
      {

        Data_access(const PointWithInfoMap& points, const ValueType * values)
              : _points(points), _values(values){};

        std::pair< ValueType, bool>
        operator()(const typename PointWithInfoMap::key_type& p) const {
         typename PointWithInfoMap::const_iterator mit = _points.find(p);
         if(mit!= _points.end())
           return std::make_pair(_values[mit->second], true);
         return std::make_pair(ValueType(), false);
        };

        const PointWithInfoMap& _points;
        const ValueType * _values;
      };

      typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
      typedef CGAL::Triangulation_vertex_base_with_info_2<uint, K> Vb;
      typedef CGAL::Triangulation_data_structure_2<Vb>             Tds;
      typedef CGAL::Delaunay_triangulation_2<K, Tds>               Delaunay;
      typedef CGAL::Interpolation_traits_2<K>                      Traits;
      typedef K::FT                                                Coord_type;
      typedef K::Point_2                                           Point;
      typedef std::map<Point, Vb::Info, K::Less_xy_2>              PointWithInfoMap;
      typedef Data_access< PointWithInfoMap, double >              Value_access;

   Delaunay fCGALdelaunay; //! CGAL delaunay triangulation object
   PointWithInfoMap fNormalizedPoints; //! Normalized function values

#endif //HAS_CGAL


};


} // namespace Math
} // namespace ROOT


#endif