Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
Delaunay2D.h
Go to the documentation of this file.
1// @(#)root/mathcore:$Id: Delaunay2D.h,v 1.00
2// Authors: Daniel Funke, Lorenzo Moneta, Olivier Couet
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 [CDT](https://github.com/artem-ogre/CDT), a C++ library for
63 generating constraint or conforming Delaunay triangulations.
64
65 After having found the triangles using the above library, barycentric coordinates are used
66 to test whether a point is inside a triangle (inTriangle test) and for interpolation.
67 All this below is implemented in the DoInterpolateNormalized function.
68
69 Given triangle ABC and point P, P can be expressed by
70
71 P.x = la * A.x + lb * B.x + lc * C.x
72 P.y = la * A.y + lb * B.y + lc * C.y
73
74 with lc = 1 - la - lb
75
76 P.x = la * A.x + lb * B.x + (1-la-lb) * C.x
77 P.y = la * A.y + lb * B.y + (1-la-lb) * C.y
78
79 Rearranging yields
80
81 la * (A.x - C.x) + lb * (B.x - C.x) = P.x - C.x
82 la * (A.y - C.y) + lb * (B.y - C.y) = P.y - C.y
83
84 Thus
85
86 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) )
87 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) )
88 lc = 1 - la - lb
89
90 We save the inverse denominator to speedup computation
91
92 invDenom = 1 / ( (B.y - C.y)*(A.x - C.x) + (C.x - B.x)*(A.y - C.y) )
93
94 P is in triangle (including edges if
95
96 0 <= [la, lb, lc] <= 1
97
98 The interpolation of P.z is
99
100 P.z = la * A.z + lb * B.z + lc * C.z
101
102 To speed up localisation of points (to see to which triangle belong) a grid is laid over the internal coordinate space.
103 A reference to triangle ABC is added to _all_ grid cells that include ABC's bounding box.
104 The size of the grid is defined to be 25x25
105
106 Optionally (if the compiler macro `HAS_GCAL` is defined ) the triangle findings and interpolation can be computed
107 using the GCAL library. This is however not supported when using the class within ROOT
108
109 \ingroup MathCore
110 */
111
112
114
115public:
116
117 struct Triangle {
118 double x[3]; // x of triangle vertices
119 double y[3]; // y of triangle vertices
120 unsigned int idx[3]; // point corresponding to vertices
121
122 #ifndef HAS_CGAL
123 double invDenom; // cached inv denominator for computing barycentric coordinates (see above)
124 #endif
125 };
126
127 typedef std::vector<Triangle> Triangles;
128
129public:
130
131
132 Delaunay2D(int n, const double *x, const double * y, const double * z, double xmin=0, double xmax=0, double ymin=0, double ymax=0);
133
134 /// set the input points for building the graph
135 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);
136
137 /// Return the Interpolated z value corresponding to the given (x,y) point.
138 /// Note that in case no Delaunay triangles are found, for example when the
139 /// points are aligned, then a default value of zero is always return.
140 /// See the class documentation for how the interpolation is computed.
141 double Interpolate(double x, double y);
142
143 /// Find all triangles
144 void FindAllTriangles();
145
146 /// return the number of triangles
147 int NumberOfTriangles() const {return fNdt;}
148
149 double XMin() const {return fXNmin;}
150 double XMax() const {return fXNmax;}
151 double YMin() const {return fYNmin;}
152 double YMax() const {return fYNmax;}
153
154 /// set z value to be returned for points outside the region
155 void SetZOuterValue(double z=0.) { fZout = z; }
156
157 /// return the user defined Z-outer value
158 double ZOuterValue() const { return fZout; }
159
160 // iterators on the found triangles
161 Triangles::const_iterator begin() const { return fTriangles.begin(); }
162 Triangles::const_iterator end() const { return fTriangles.end(); }
163
164
165private:
166
167 // internal methods
168
169
170 inline double Linear_transform(double x, double offset, double factor){
171 return (x+offset)*factor;
172 }
173
174 /// internal function to normalize the points.
175 /// See the class documentation for the details on how it is computed.
176 void DoNormalizePoints();
177
178 /// internal function to find the triangle
179 /// use Triangle or CGAL if flag is set
180 void DoFindTriangles();
181
182 /// internal method to compute the interpolation
183 double DoInterpolateNormalized(double x, double y);
184
185
186
187private:
188 // class is not copyable
189 Delaunay2D(const Delaunay2D&); // Not implemented
190 Delaunay2D& operator=(const Delaunay2D&); // Not implemented
191
192protected:
193
194 int fNdt; ///<! Number of Delaunay triangles found
195 int fNpoints; ///<! Number of data points
196
197 const double *fX; ///<! Pointer to X array (managed externally)
198 const double *fY; ///<! Pointer to Y array
199 const double *fZ; ///<! Pointer to Z array
200
201 double fXNmin; ///<! Minimum value of fXN
202 double fXNmax; ///<! Maximum value of fXN
203 double fYNmin; ///<! Minimum value of fYN
204 double fYNmax; ///<! Maximum value of fYN
205
206 double fOffsetX; ///<! Normalization offset X
207 double fOffsetY; ///<! Normalization offset Y
208
209 double fScaleFactorX; ///<! Normalization factor X
210 double fScaleFactorY; ///<! Normalization factor Y
211
212 double fZout; ///<! Height for points lying outside the convex hull
213
214 bool fInit; ///<! True if FindAllTriangles() has been performed
215
216
217 Triangles fTriangles; ///<! Triangles of Triangulation
218
219#ifndef HAS_CGAL
220
221 //using triangle library
222
223 std::vector<double> fXN; ///<! normalized X
224 std::vector<double> fYN; ///<! normalized Y
225
226 static const int fNCells = 25; ///<! number of cells to divide the normalized space
227 double fXCellStep; ///<! inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
228 double fYCellStep; ///<! inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
229 std::set<unsigned int> fCells[(fNCells+1)*(fNCells+1)]; ///<! grid cells with containing triangles
230
231 inline unsigned int Cell(unsigned int x, unsigned int y) const {
232 return x*(fNCells+1) + y;
233 }
234
235 inline int CellX(double x) const {
236 return (x - fXNmin) * fXCellStep;
237 }
238
239 inline int CellY(double y) const {
240 return (y - fYNmin) * fYCellStep;
241 }
242
243#else // HAS_CGAL
244 // case of using GCAL
245 //Functor class for accessing the function values/gradients
246 template< class PointWithInfoMap, typename ValueType >
247 struct Data_access : public std::unary_function< typename PointWithInfoMap::key_type,
248 std::pair<ValueType, bool> >
249 {
250
251 Data_access(const PointWithInfoMap& points, const ValueType * values)
252 : _points(points), _values(values){};
253
254 std::pair< ValueType, bool>
255 operator()(const typename PointWithInfoMap::key_type& p) const {
256 typename PointWithInfoMap::const_iterator mit = _points.find(p);
257 if(mit!= _points.end())
258 return std::make_pair(_values[mit->second], true);
259 return std::make_pair(ValueType(), false);
260 };
261
262 const PointWithInfoMap& _points;
263 const ValueType * _values;
264 };
265
266 typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
267 typedef CGAL::Triangulation_vertex_base_with_info_2<uint, K> Vb;
268 typedef CGAL::Triangulation_data_structure_2<Vb> Tds;
269 typedef CGAL::Delaunay_triangulation_2<K, Tds> Delaunay;
270 typedef CGAL::Interpolation_traits_2<K> Traits;
271 typedef K::FT Coord_type;
272 typedef K::Point_2 Point;
273 typedef std::map<Point, Vb::Info, K::Less_xy_2> PointWithInfoMap;
274 typedef Data_access< PointWithInfoMap, double > Value_access;
275
276 Delaunay fCGALdelaunay; //! CGAL delaunay triangulation object
277 PointWithInfoMap fNormalizedPoints; //! Normalized function values
278
279#endif //HAS_CGAL
280
281
282};
283
284
285} // namespace Math
286} // namespace ROOT
287
288
289#endif
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h offset
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
float xmin
float ymin
float xmax
float ymax
TRObject operator()(const T1 &t1) const
Class to generate a Delaunay triangulation of a 2D set of points.
Definition Delaunay2D.h:113
Triangles::const_iterator end() const
Definition Delaunay2D.h:162
int fNpoints
! Number of data points
Definition Delaunay2D.h:195
double fXNmin
! Minimum value of fXN
Definition Delaunay2D.h:201
std::set< unsigned int > fCells[(fNCells+1) *(fNCells+1)]
! grid cells with containing triangles
Definition Delaunay2D.h:229
int fNdt
! Number of Delaunay triangles found
Definition Delaunay2D.h:194
double DoInterpolateNormalized(double x, double y)
internal method to compute the interpolation
int CellY(double y) const
Definition Delaunay2D.h:239
Triangles::const_iterator begin() const
Definition Delaunay2D.h:161
Delaunay2D(const Delaunay2D &)
const double * fZ
! Pointer to Z array
Definition Delaunay2D.h:199
double fXNmax
! Maximum value of fXN
Definition Delaunay2D.h:202
void SetZOuterValue(double z=0.)
set z value to be returned for points outside the region
Definition Delaunay2D.h:155
double fZout
! Height for points lying outside the convex hull
Definition Delaunay2D.h:212
double YMin() const
Definition Delaunay2D.h:151
double fOffsetY
! Normalization offset Y
Definition Delaunay2D.h:207
Delaunay2D & operator=(const Delaunay2D &)
double XMax() const
Definition Delaunay2D.h:150
double ZOuterValue() const
return the user defined Z-outer value
Definition Delaunay2D.h:158
Triangles fTriangles
! Triangles of Triangulation
Definition Delaunay2D.h:217
double Linear_transform(double x, double offset, double factor)
Definition Delaunay2D.h:170
double fYCellStep
! inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
Definition Delaunay2D.h:228
std::vector< double > fYN
! normalized Y
Definition Delaunay2D.h:224
void FindAllTriangles()
Find all triangles.
static const int fNCells
! number of cells to divide the normalized space
Definition Delaunay2D.h:226
double YMax() const
Definition Delaunay2D.h:152
bool fInit
! True if FindAllTriangles() has been performed
Definition Delaunay2D.h:214
unsigned int Cell(unsigned int x, unsigned int y) const
Definition Delaunay2D.h:231
std::vector< double > fXN
! normalized X
Definition Delaunay2D.h:223
double fXCellStep
! inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
Definition Delaunay2D.h:227
double fYNmin
! Minimum value of fYN
Definition Delaunay2D.h:203
double fOffsetX
! Normalization offset X
Definition Delaunay2D.h:206
double XMin() const
Definition Delaunay2D.h:149
std::vector< Triangle > Triangles
Definition Delaunay2D.h:127
const double * fY
! Pointer to Y array
Definition Delaunay2D.h:198
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
const double * fX
! Pointer to X array (managed externally)
Definition Delaunay2D.h:197
void DoNormalizePoints()
internal function to normalize the points.
double fYNmax
! Maximum value of fYN
Definition Delaunay2D.h:204
int NumberOfTriangles() const
return the number of triangles
Definition Delaunay2D.h:147
double fScaleFactorY
! Normalization factor Y
Definition Delaunay2D.h:210
double fScaleFactorX
! Normalization factor X
Definition Delaunay2D.h:209
void DoFindTriangles()
internal function to find the triangle use Triangle or CGAL if flag is set
int CellX(double x) const
Definition Delaunay2D.h:235
#define Interpolate(a, x, b, y)
Definition geom.c:179
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.
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
constexpr Double_t K()
Boltzmann's constant in : .
Definition TMath.h:247