root/Math/Delaunay2D.h

0001 // @(#)root/mathcore:$Id: Delaunay2D.h,v 1.00
0002 // Authors: Daniel Funke, Lorenzo Moneta, Olivier Couet
0003
0004 /*************************************************************************
0005  * Copyright (C) 2015 ROOT Math Team                                     *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011
0012 // Header file for class Delaunay2D
0013
0014 #ifndef ROOT_Math_Delaunay2D
0015 #define ROOT_Math_Delaunay2D
0016
0017 //for testing purposes HAS_CGAL can be [un]defined here
0018 //#define HAS_CGAL
0019
0020 //for testing purposes THREAD_SAFE can [un]defined here
0021 //#define THREAD_SAFE
0022
0023
0024 //#include "RtypesCore.h"
0025
0026 #include <map>
0027 #include <vector>
0028 #include <set>
0029 #include <functional>
0030
0031 #ifdef HAS_CGAL
0032    /* CGAL uses the name PTR as member name in its Handle class
0033     * but its a macro defined in mmalloc.h of ROOT
0034     * Safe it, disable it and then re-enable it later on*/
0035    #pragma push_macro("PTR")
0036    #undef PTR
0037
0038    #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
0039    #include <CGAL/Delaunay_triangulation_2.h>
0040    #include <CGAL/Triangulation_vertex_base_with_info_2.h>
0041    #include <CGAL/Interpolation_traits_2.h>
0042    #include <CGAL/natural_neighbor_coordinates_2.h>
0043    #include <CGAL/interpolation_functions.h>
0044
0045    #pragma pop_macro("PTR")
0046 #endif
0047
0048 #ifdef THREAD_SAFE
0049    #include<atomic> //atomic operations for thread safety
0050 #endif
0051
0052
0053 namespace ROOT {
0054
0055
0056
0057    namespace Math {
0058
0059 /**
0060
0061    Class to generate a Delaunay triangulation of a 2D set of points.
0062    Algorithm based on [CDT](https://github.com/artem-ogre/CDT), a C++ library for
0063    generating constraint or conforming Delaunay triangulations.
0064
0065    After having found the triangles using the above library,  barycentric coordinates are used
0066    to test whether a point is inside a triangle (inTriangle test) and for interpolation.
0067    All this below is implemented in the DoInterpolateNormalized function.
0068
0069    Given triangle ABC and point P, P can be expressed by
0070
0071      P.x = la * A.x + lb * B.x + lc * C.x
0072      P.y = la * A.y + lb * B.y + lc * C.y
0073
0074    with lc = 1 - la - lb
0075
0076      P.x = la * A.x + lb * B.x + (1-la-lb) * C.x
0077      P.y = la * A.y + lb * B.y + (1-la-lb) * C.y
0078
0079    Rearranging yields
0080
0081      la * (A.x - C.x) + lb * (B.x - C.x) = P.x - C.x
0082      la * (A.y - C.y) + lb * (B.y - C.y) = P.y - C.y
0083
0084    Thus
0085
0086      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) )
0087      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) )
0088      lc = 1 - la - lb
0089
0090    We save the inverse denominator to speedup computation
0091
0092      invDenom = 1 / ( (B.y - C.y)*(A.x - C.x) + (C.x - B.x)*(A.y - C.y) )
0093
0094      P is in triangle (including edges if
0095
0096      0 <= [la, lb, lc] <= 1
0097
0098      The interpolation of P.z is
0099
0100      P.z = la * A.z + lb * B.z + lc * C.z
0101
0102    To speed up localisation of points (to see to which triangle belong) a grid is laid over the internal coordinate space.
0103    A reference to triangle ABC is added to _all_ grid cells that include ABC's bounding box.
0104    The size of the grid is defined to be 25x25
0105
0106    Optionally (if the compiler macro `HAS_GCAL` is defined ) the triangle findings and interpolation can be computed
0107    using the GCAL library. This is however not supported when using the class within ROOT
0108
0109    \ingroup MathCore
0110  */
0111
0112
0113 class Delaunay2D  {
0114
0115 public:
0116
0117    struct Triangle {
0118       double x[3];           // x of triangle vertices
0119       double y[3];           // y of triangle vertices
0120       unsigned int idx[3];   // point corresponding to vertices
0121
0122       #ifndef HAS_CGAL
0123       double invDenom; // cached inv denominator for computing barycentric coordinates (see above)
0124       #endif
0125    };
0126
0127    typedef std::vector<Triangle> Triangles;
0128
0129 public:
0130
0131
0132    Delaunay2D(int n, const double *x, const double * y, const double * z, double xmin=0, double xmax=0, double ymin=0, double ymax=0);
0133
0134    /// set the input points for building the graph
0135    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);
0136
0137    /// Return the Interpolated z value corresponding to the given (x,y) point.
0138    /// Note that in case no Delaunay triangles are found, for example when the
0139    /// points are aligned, then a default value of zero is always return.
0140    /// See the class documentation for  how the interpolation is computed.
0141    double  Interpolate(double x, double y);
0142
0143    /// Find all triangles
0144    void      FindAllTriangles();
0145
0146    /// return the number of triangles
0147    int    NumberOfTriangles() const {return fNdt;}
0148
0149    double  XMin() const {return fXNmin;}
0150    double  XMax() const {return fXNmax;}
0151    double  YMin() const {return fYNmin;}
0152    double  YMax() const {return fYNmax;}
0153
0154    /// set z value to be returned for points  outside the region
0155    void      SetZOuterValue(double z=0.) { fZout = z; }
0156
0157    /// return the user defined Z-outer value
0158    double ZOuterValue() const { return fZout; }
0159
0160    // iterators on the found triangles
0161    Triangles::const_iterator begin() const { return fTriangles.begin(); }
0162    Triangles::const_iterator end()  const { return fTriangles.end(); }
0163
0164
0165 private:
0166
0167    // internal methods
0168
0169
0170    inline double Linear_transform(double x, double offset, double factor){
0171       return (x+offset)*factor;
0172    }
0173
0174    /// internal function to normalize the points.
0175    /// See the class documentation for the details on how it is computed.
0176    void DoNormalizePoints();
0177
0178    /// internal function to find the triangle
0179    /// use Triangle or CGAL if flag is set
0180    void DoFindTriangles();
0181
0182    /// internal method to compute the interpolation
0183    double  DoInterpolateNormalized(double x, double y);
0184
0185
0186
0187 private:
0188    // class is not copyable
0189    Delaunay2D(const Delaunay2D&); // Not implemented
0190    Delaunay2D& operator=(const Delaunay2D&); // Not implemented
0191
0192 protected:
0193
0194    int         fNdt;           ///<! Number of Delaunay triangles found
0195    int         fNpoints;       ///<! Number of data points
0196
0197    const double   *fX;         ///<! Pointer to X array (managed externally)
0198    const double   *fY;         ///<! Pointer to Y array
0199    const double   *fZ;         ///<! Pointer to Z array
0200
0201    double    fXNmin;           ///<! Minimum value of fXN
0202    double    fXNmax;           ///<! Maximum value of fXN
0203    double    fYNmin;           ///<! Minimum value of fYN
0204    double    fYNmax;           ///<! Maximum value of fYN
0205
0206    double    fOffsetX;         ///<! Normalization offset X
0207    double    fOffsetY;         ///<! Normalization offset Y
0208
0209    double    fScaleFactorX;    ///<! Normalization factor X
0210    double    fScaleFactorY;    ///<! Normalization factor Y
0211
0212    double    fZout;            ///<! Height for points lying outside the convex hull
0213
0214    bool      fInit;            ///<! True if FindAllTriangles() has been performed
0215
0216
0217    Triangles   fTriangles;     ///<! Triangles of Triangulation
0218
0219 #ifndef HAS_CGAL
0220
0221    //using triangle library
0222
0223    std::vector<double> fXN; ///<! normalized X
0224    std::vector<double> fYN; ///<! normalized Y
0225
0226    static const int fNCells = 25; ///<! number of cells to divide the normalized space
0227    double fXCellStep; ///<! inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
0228    double fYCellStep; ///<! inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
0229    std::set<unsigned int> fCells[(fNCells+1)*(fNCells+1)]; ///<! grid cells with containing triangles
0230
0231    inline unsigned int Cell(unsigned int x, unsigned int y) const {
0232       return x*(fNCells+1) + y;
0233    }
0234
0235    inline int CellX(double x) const {
0236       return (x - fXNmin) * fXCellStep;
0237    }
0238
0239    inline int CellY(double y) const {
0240       return (y - fYNmin) * fYCellStep;
0241    }
0242
0243 #else // HAS_CGAL
0244       // case of using GCAL
0245       //Functor class for accessing the function values/gradients
0246       template< class PointWithInfoMap, typename ValueType >
0247       struct Data_access : public std::unary_function< typename PointWithInfoMap::key_type,
0248                 std::pair<ValueType, bool> >
0249       {
0250
0251         Data_access(const PointWithInfoMap& points, const ValueType * values)
0252               : _points(points), _values(values){};
0253
0254         std::pair< ValueType, bool>
0255         operator()(const typename PointWithInfoMap::key_type& p) const {
0256          typename PointWithInfoMap::const_iterator mit = _points.find(p);
0257          if(mit!= _points.end())
0258            return std::make_pair(_values[mit->second], true);
0259          return std::make_pair(ValueType(), false);
0260         };
0261
0262         const PointWithInfoMap& _points;
0263         const ValueType * _values;
0264       };
0265
0266       typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
0267       typedef CGAL::Triangulation_vertex_base_with_info_2<uint, K> Vb;
0268       typedef CGAL::Triangulation_data_structure_2<Vb>             Tds;
0269       typedef CGAL::Delaunay_triangulation_2<K, Tds>               Delaunay;
0270       typedef CGAL::Interpolation_traits_2<K>                      Traits;
0271       typedef K::FT                                                Coord_type;
0272       typedef K::Point_2                                           Point;
0273       typedef std::map<Point, Vb::Info, K::Less_xy_2>              PointWithInfoMap;
0274       typedef Data_access< PointWithInfoMap, double >              Value_access;
0275
0276    Delaunay fCGALdelaunay; //! CGAL delaunay triangulation object
0277    PointWithInfoMap fNormalizedPoints; //! Normalized function values
0278
0279 #endif //HAS_CGAL
0280
0281
0282 };
0283
0284
0285 } // namespace Math
0286 } // namespace ROOT
0287
0288
0289 #endif