Logo ROOT  
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 "RtypesCore.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
53namespace 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
72
73public:
74
75 struct Triangle {
76 double x[3];
77 double y[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
87public:
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
120private:
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
141private:
142 // class is not copyable
143 Delaunay2D(const Delaunay2D&); // Not implemented
144 Delaunay2D& operator=(const Delaunay2D&); // Not implemented
145
146protected:
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
int Int_t
Definition: RtypesCore.h:45
unsigned int UInt_t
Definition: RtypesCore.h:46
bool Bool_t
Definition: RtypesCore.h:63
float xmin
Definition: THbookFile.cxx:95
float ymin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
float ymax
Definition: THbookFile.cxx:95
TRObject operator()(const T1 &t1) const
point * points
Definition: X3DBuffer.c:22
Class to generate a Delaunay triangulation of a 2D set of points.
Definition: Delaunay2D.h:71
Triangles::const_iterator end() const
Definition: Delaunay2D.h:117
double fXNmin
Pointer to Z array.
Definition: Delaunay2D.h:157
double DoInterpolateNormalized(double x, double y)
internal method to compute the interpolation
Definition: Delaunay2D.cxx:379
int CellY(double y) const
Definition: Delaunay2D.h:282
Triangles::const_iterator begin() const
Definition: Delaunay2D.h:116
Delaunay2D(const Delaunay2D &)
double Interpolate(double x, double y)
Return the Interpolated z value corresponding to the (x,y) point.
Definition: Delaunay2D.cxx:104
const double * fZ
Pointer to Y array.
Definition: Delaunay2D.h:155
double fXNmax
Minimum value of fXN.
Definition: Delaunay2D.h:158
void SetZOuterValue(double z=0.)
set z value to be returned for points outside the region
Definition: Delaunay2D.h:110
double fZout
Normalization factor Y.
Definition: Delaunay2D.h:171
Int_t NumberOfTriangles() const
return the number of triangles
Definition: Delaunay2D.h:102
double YMin() const
Definition: Delaunay2D.h:106
double fOffsetY
Normalization offset X.
Definition: Delaunay2D.h:166
Delaunay2D & operator=(const Delaunay2D &)
double XMax() const
Definition: Delaunay2D.h:105
double ZOuterValue() const
return the user defined Z-outer value
Definition: Delaunay2D.h:113
Triangles fTriangles
True if FindAllTriangels() has been performed.
Definition: Delaunay2D.h:183
unsigned int Cell(UInt_t x, UInt_t y) const
grid cells with containing triangles
Definition: Delaunay2D.h:274
double Linear_transform(double x, double offset, double factor)
Definition: Delaunay2D.h:125
double fYCellStep
inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
Definition: Delaunay2D.h:271
std::vector< double > fYN
normalized X
Definition: Delaunay2D.h:262
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
void FindAllTriangles()
Find all triangles.
Definition: Delaunay2D.cxx:128
static const int fNCells
normalized Y
Definition: Delaunay2D.h:269
double YMax() const
Definition: Delaunay2D.h:107
std::vector< double > fXN
Triangles of Triangulation.
Definition: Delaunay2D.h:261
double fXCellStep
number of cells to divide the normalized space
Definition: Delaunay2D.h:270
std::set< UInt_t > fCells[(fNCells+1) *(fNCells+1)]
inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
Definition: Delaunay2D.h:272
double fYNmin
Maximum value of fXN.
Definition: Delaunay2D.h:159
double fOffsetX
Maximum value of fYN.
Definition: Delaunay2D.h:165
Int_t fNpoints
Number of Delaunay triangles found.
Definition: Delaunay2D.h:151
double XMin() const
Definition: Delaunay2D.h:104
std::vector< Triangle > Triangles
Definition: Delaunay2D.h:85
Bool_t fInit
Height for points lying outside the convex hull.
Definition: Delaunay2D.h:179
const double * fY
Pointer to X array (managed externally)
Definition: Delaunay2D.h:154
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
const double * fX
Number of data points.
Definition: Delaunay2D.h:153
void DoNormalizePoints()
internal function to normalize the points
Definition: Delaunay2D.cxx:244
double fYNmax
Minimum value of fYN.
Definition: Delaunay2D.h:160
double fScaleFactorY
Normalization factor X.
Definition: Delaunay2D.h:169
double fScaleFactorX
Normalization offset Y.
Definition: Delaunay2D.h:168
void DoFindTriangles()
internal function to find the triangle use Triangle or CGAL if flag is set
Definition: Delaunay2D.cxx:256
int CellX(double x) const
Definition: Delaunay2D.h:278
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
Namespace for new Math classes and functions.
static const uint32_t K[64]
Definition: RSha256.hxx:148
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
static const Int_t UNINITIALIZED
Definition: TNeuron.cxx:41