Logo ROOT   6.13/01
Reference Guide
Delaunay2D.h
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id: Delaunay2D.h,v 1.00
2 // Author: Daniel Funke, Lorenzo Moneta
3 
4 /*************************************************************************
5  * Copyright (C) 2015 ROOT Math Team *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 // Header file for class Delaunay2D
13 
14 #ifndef ROOT_Math_Delaunay2D
15 #define ROOT_Math_Delaunay2D
16 
17 //for testing purposes HAS_CGAL can be [un]defined here
18 //#define HAS_CGAL
19 
20 //for testing purposes THREAD_SAFE can [un]defined here
21 //#define THREAD_SAFE
22 
23 
24 #include "TNamed.h"
25 
26 #include <map>
27 #include <vector>
28 #include <set>
29 #include <functional>
30 
31 #ifdef HAS_CGAL
32  /* CGAL uses the name PTR as member name in its Handle class
33  * but its a macro defined in mmalloc.h of ROOT
34  * Safe it, disable it and then re-enable it later on*/
35  #pragma push_macro("PTR")
36  #undef PTR
37 
38  #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
39  #include <CGAL/Delaunay_triangulation_2.h>
40  #include <CGAL/Triangulation_vertex_base_with_info_2.h>
41  #include <CGAL/Interpolation_traits_2.h>
42  #include <CGAL/natural_neighbor_coordinates_2.h>
43  #include <CGAL/interpolation_functions.h>
44 
45  #pragma pop_macro("PTR")
46 #endif
47 
48 #ifdef THREAD_SAFE
49  #include<atomic> //atomic operations for thread safety
50 #endif
51 
52 
53 namespace ROOT {
54 
55 
56 
57  namespace Math {
58 
59 /**
60 
61  Class to generate a Delaunay triangulation of a 2D set of points.
62  Algorithm based on **Triangle**, a two-dimensional quality mesh generator and
63  Delaunay triangulator from Jonathan Richard Shewchuk.
64 
65  See [http://www.cs.cmu.edu/~quake/triangle.html]
66 
67  \ingroup MathCore
68  */
69 
70 
71 class Delaunay2D {
72 
73 public:
74 
75  struct Triangle {
76  double x[3];
77  double y[3];
78  UInt_t idx[3];
79 
80  #ifndef HAS_CGAL
81  double invDenom; //see comment below in CGAL fall back section
82  #endif
83  };
84 
85  typedef std::vector<Triangle> Triangles;
86 
87 public:
88 
89 
90  Delaunay2D(int n, const double *x, const double * y, const double * z, double xmin=0, double xmax=0, double ymin=0, double ymax=0);
91 
92  /// set the input points for building the graph
93  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);
94 
95  /// Return the Interpolated z value corresponding to the (x,y) point
96  double Interpolate(double x, double y);
97 
98  /// Find all triangles
99  void FindAllTriangles();
100 
101  /// return the number of triangles
102  Int_t NumberOfTriangles() const {return fNdt;}
103 
104  double XMin() const {return fXNmin;}
105  double XMax() const {return fXNmax;}
106  double YMin() const {return fYNmin;}
107  double YMax() const {return fYNmax;}
108 
109  /// set z value to be returned for points outside the region
110  void SetZOuterValue(double z=0.) { fZout = z; }
111 
112  /// return the user defined Z-outer value
113  double ZOuterValue() const { return fZout; }
114 
115  // iterators on the found triangles
116  Triangles::const_iterator begin() const { return fTriangles.begin(); }
117  Triangles::const_iterator end() const { return fTriangles.end(); }
118 
119 
120 private:
121 
122  // internal methods
123 
124 
125  inline double Linear_transform(double x, double offset, double factor){
126  return (x+offset)*factor;
127  }
128 
129  /// internal function to normalize the points
130  void DoNormalizePoints();
131 
132  /// internal function to find the triangle
133  /// use Triangle or CGAL if flag is set
134  void DoFindTriangles();
135 
136  /// internal method to compute the interpolation
137  double DoInterpolateNormalized(double x, double y);
138 
139 
140 
141 private:
142  // class is not copyable
143  Delaunay2D(const Delaunay2D&); // Not implemented
144  Delaunay2D& operator=(const Delaunay2D&); // Not implemented
145 
146 protected:
147 
148  //typedef std::function<double(double)> Transformer;
149 
150  Int_t fNdt; //!Number of Delaunay triangles found
151  Int_t fNpoints; //!Number of data points
152 
153  const double *fX; //!Pointer to X array (managed externally)
154  const double *fY; //!Pointer to Y array
155  const double *fZ; //!Pointer to Z array
156 
157  double fXNmin; //!Minimum value of fXN
158  double fXNmax; //!Maximum value of fXN
159  double fYNmin; //!Minimum value of fYN
160  double fYNmax; //!Maximum value of fYN
161 
162  //Transformer xTransformer; //!transform x values to mapped space
163  //Transformer yTransformer; //!transform y values to mapped space
164 
165  double fOffsetX; //!Normalization offset X
166  double fOffsetY; //!Normalization offset Y
167 
168  double fScaleFactorX; //!Normalization factor X
169  double fScaleFactorY; //!Normalization factor Y
170 
171  double fZout; //!Height for points lying outside the convex hull
172 
173 #ifdef THREAD_SAFE
174 
175  enum class Initialization : char {UNINITIALIZED, INITIALIZING, INITIALIZED};
176  std::atomic<Initialization> fInit; //!Indicate initialization state
177 
178 #else
179  Bool_t fInit; //!True if FindAllTriangels() has been performed
180 #endif
181 
182 
183  Triangles fTriangles; //!Triangles of Triangulation
184 
185 #ifdef HAS_CGAL
186 
187  //Functor class for accessing the function values/gradients
188  template< class PointWithInfoMap, typename ValueType >
189  struct Data_access : public std::unary_function< typename PointWithInfoMap::key_type,
190  std::pair<ValueType, bool> >
191  {
192 
193  Data_access(const PointWithInfoMap& points, const ValueType * values)
194  : _points(points), _values(values){};
195 
196  std::pair< ValueType, bool>
197  operator()(const typename PointWithInfoMap::key_type& p) const {
198  typename PointWithInfoMap::const_iterator mit = _points.find(p);
199  if(mit!= _points.end())
200  return std::make_pair(_values[mit->second], true);
201  return std::make_pair(ValueType(), false);
202  };
203 
204  const PointWithInfoMap& _points;
205  const ValueType * _values;
206  };
207 
208  typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
209  typedef CGAL::Triangulation_vertex_base_with_info_2<uint, K> Vb;
210  typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
211  typedef CGAL::Delaunay_triangulation_2<K, Tds> Delaunay;
212  typedef CGAL::Interpolation_traits_2<K> Traits;
213  typedef K::FT Coord_type;
214  typedef K::Point_2 Point;
215  typedef std::map<Point, Vb::Info, K::Less_xy_2> PointWithInfoMap;
216  typedef Data_access< PointWithInfoMap, double > Value_access;
217 
218  Delaunay fCGALdelaunay; //! CGAL delaunay triangulation object
219  PointWithInfoMap fNormalizedPoints; //! Normalized function values
220 
221 #else // HAS_CGAL
222  //fallback to triangle library
223 
224  /* Using barycentric coordinates for inTriangle test and interpolation
225  *
226  * Given triangle ABC and point P, P can be expressed by
227  *
228  * P.x = la * A.x + lb * B.x + lc * C.x
229  * P.y = la * A.y + lb * B.y + lc * C.y
230  *
231  * with lc = 1 - la - lb
232  *
233  * P.x = la * A.x + lb * B.x + (1-la-lb) * C.x
234  * P.y = la * A.y + lb * B.y + (1-la-lb) * C.y
235  *
236  * Rearranging yields
237  *
238  * la * (A.x - C.x) + lb * (B.x - C.x) = P.x - C.x
239  * la * (A.y - C.y) + lb * (B.y - C.y) = P.y - C.y
240  *
241  * Thus
242  *
243  * 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) )
244  * 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) )
245  * lc = 1 - la - lb
246  *
247  * We save the inverse denominator to speedup computation
248  *
249  * invDenom = 1 / ( (B.y - C.y)*(A.x - C.x) + (C.x - B.x)*(A.y - C.y) )
250  *
251  * P is in triangle (including edges if
252  *
253  * 0 <= [la, lb, lc] <= 1
254  *
255  * The interpolation of P.z is
256  *
257  * P.z = la * A.z + lb * B.z + lc * C.z
258  *
259  */
260 
261  std::vector<double> fXN; //! normalized X
262  std::vector<double> fYN; //! normalized Y
263 
264  /* To speed up localisation of points a grid is layed over normalized space
265  *
266  * A reference to triangle ABC is added to _all_ grid cells that include ABC's bounding box
267  */
268 
269  static const int fNCells = 25; //! number of cells to divide the normalized space
270  double fXCellStep; //! inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
271  double fYCellStep; //! inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
272  std::set<UInt_t> fCells[(fNCells+1)*(fNCells+1)]; //! grid cells with containing triangles
273 
274  inline unsigned int Cell(UInt_t x, UInt_t y) const {
275  return x*(fNCells+1) + y;
276  }
277 
278  inline int CellX(double x) const {
279  return (x - fXNmin) * fXCellStep;
280  }
281 
282  inline int CellY(double y) const {
283  return (y - fYNmin) * fYCellStep;
284  }
285 
286 #endif //HAS_CGAL
287 
288 
289 };
290 
291 
292 } // namespace Math
293 } // namespace ROOT
294 
295 
296 #endif
unsigned int Cell(UInt_t x, UInt_t y) const
grid cells with containing triangles
Definition: Delaunay2D.h:274
double fZout
Normalization factor Y.
Definition: Delaunay2D.h:171
void DoNormalizePoints()
internal function to normalize the points
Definition: Delaunay2D.cxx:244
std::vector< Triangle > Triangles
Definition: Delaunay2D.h:85
constexpr Double_t K()
Boltzmann&#39;s constant in .
Definition: TMath.h:247
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
Bool_t fInit
Height for points lying outside the convex hull.
Definition: Delaunay2D.h:179
double XMax() const
Definition: Delaunay2D.h:105
void FindAllTriangles()
Find all triangles.
Definition: Delaunay2D.cxx:128
int CellY(double y) const
Definition: Delaunay2D.h:282
void SetZOuterValue(double z=0.)
set z value to be returned for points outside the region
Definition: Delaunay2D.h:110
double fScaleFactorX
Normalization offset Y.
Definition: Delaunay2D.h:168
double DoInterpolateNormalized(double x, double y)
internal method to compute the interpolation
Definition: Delaunay2D.cxx:379
double fOffsetY
Normalization offset X.
Definition: Delaunay2D.h:166
double fScaleFactorY
Normalization factor X.
Definition: Delaunay2D.h:169
double fOffsetX
Maximum value of fYN.
Definition: Delaunay2D.h:165
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)
set the input points for building the graph
Definition: Delaunay2D.cxx:59
std::set< UInt_t > fCells[(fNCells+1) *(fNCells+1)]
inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
Definition: Delaunay2D.h:272
double fXNmin
Pointer to Z array.
Definition: Delaunay2D.h:157
Triangles::const_iterator begin() const
Definition: Delaunay2D.h:116
const double * fY
Pointer to X array (managed externally)
Definition: Delaunay2D.h:154
double fXNmax
Minimum value of fXN.
Definition: Delaunay2D.h:158
Delaunay2D & operator=(const Delaunay2D &)
std::vector< double > fYN
normalized X
Definition: Delaunay2D.h:262
double fYNmax
Minimum value of fYN.
Definition: Delaunay2D.h:160
double fXCellStep
number of cells to divide the normalized space
Definition: Delaunay2D.h:270
double YMax() const
Definition: Delaunay2D.h:107
Int_t fNpoints
Number of Delaunay triangles found.
Definition: Delaunay2D.h:151
double fYCellStep
inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
Definition: Delaunay2D.h:271
double ZOuterValue() const
return the user defined Z-outer value
Definition: Delaunay2D.h:113
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
double YMin() const
Definition: Delaunay2D.h:106
static const int fNCells
normalized Y
Definition: Delaunay2D.h:269
void DoFindTriangles()
internal function to find the triangle use Triangle or CGAL if flag is set
Definition: Delaunay2D.cxx:256
std::vector< double > fXN
Triangles of Triangulation.
Definition: Delaunay2D.h:261
Int_t NumberOfTriangles() const
return the number of triangles
Definition: Delaunay2D.h:102
double XMin() const
Definition: Delaunay2D.h:104
const double * fX
Number of data points.
Definition: Delaunay2D.h:153
Class to generate a Delaunay triangulation of a 2D set of points.
Definition: Delaunay2D.h:71
double fYNmin
Maximum value of fXN.
Definition: Delaunay2D.h:159
Triangles::const_iterator end() const
Definition: Delaunay2D.h:117
int CellX(double x) const
Definition: Delaunay2D.h:278
double Interpolate(double x, double y)
Return the Interpolated z value corresponding to the (x,y) point.
Definition: Delaunay2D.cxx:104
Delaunay2D(int n, const double *x, const double *y, const double *z, double xmin=0, double xmax=0, double ymin=0, double ymax=0)
class constructor from array of data points
Definition: Delaunay2D.cxx:34
double Linear_transform(double x, double offset, double factor)
Definition: Delaunay2D.h:125
const double * fZ
Pointer to Y array.
Definition: Delaunay2D.h:155
Triangles fTriangles
True if FindAllTriangels() has been performed.
Definition: Delaunay2D.h:183