ROOT   6.08/07 Reference Guide
Delaunay2D.cxx
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 // Implementation file for class Delaunay2D
13
14 #include "Math/Delaunay2D.h"
15 #include "Rtypes.h"
16
18
19 // in case we do not use CGAL
20 #ifndef HAS_CGAL
21 // use the triangle library
22 #include "triangle.h"
23 #endif
24
25 #include <algorithm>
26 #include <stdlib.h>
27
28 namespace ROOT {
29
30  namespace Math {
31
32
33 /// class constructor from array of data points
34 Delaunay2D::Delaunay2D(int n, const double * x, const double * y, const double * z,
35  double xmin, double xmax, double ymin, double ymax)
36 {
37  // Delaunay2D normal constructor
38
39  fX = x;
40  fY = y;
41  fZ = z;
42  fNpoints = n;
43  fOffsetX = 0;
44  fOffsetY = 0;
45  fScaleFactorX = 0;
46  fScaleFactorY = 0;
47  fNdt = 0;
48  fXNmax = 0;
49  fXNmin = 0;
50  fYNmin = 0;
51  fYNmax = 0;
52
53  SetInputPoints(n,x,y,z,xmin,xmax,ymin,ymax);
54
55 }
56
57 /// set the input points
58 void Delaunay2D::SetInputPoints(int n, const double * x, const double * y, const double * z,
59  double xmin, double xmax, double ymin, double ymax) {
60
61
64 #else
65  fInit = kFALSE;
66 #endif
67
68  if (n == 0 || !x || !y || !z ) return;
69
70  if (xmin >= xmax) {
71  xmin = *std::min_element(x, x+n);
72  xmax = *std::max_element(x, x+n);
73
74  ymin = *std::min_element(y, y+n);
75  ymax = *std::max_element(y, y+n);
76  }
77
78  fOffsetX = -(xmax+xmin)/2.;
79  fOffsetY = -(ymax+ymin)/2.;
80
81  fScaleFactorX = 1./(xmax-xmin);
82  fScaleFactorY = 1./(ymax-ymin);
83
84  //xTransformer = std::bind(linear_transform, std::placeholders::_1, Xoffset, XScaleFactor);
85  //yTransformer = std::bind(linear_transform, std::placeholders::_1, Yoffset, YScaleFactor);
86
87  fXNmax = Linear_transform(xmax, fOffsetX, fScaleFactorX); //xTransformer(xmax);
88  fXNmin = Linear_transform(xmin, fOffsetX, fScaleFactorX); //xTransformer(xmin);
89
90  fYNmax = Linear_transform(ymax, fOffsetY, fScaleFactorY); //yTransformer(ymax);
91  fYNmin = Linear_transform(ymin, fOffsetY, fScaleFactorY); //yTransformer(ymin);
92
93  //printf("Normalized space extends from (%f,%f) to (%f,%f)\n", fXNmin, fYNmin, fXNmax, fYNmax);
94
95
96 #ifndef HAS_CGAL
97  fXCellStep = 0.;
98  fYCellStep = 0.;
99 #endif
100 }
101
102 //______________________________________________________________________________
103 double Delaunay2D::Interpolate(double x, double y)
104 {
105  // Return the z value corresponding to the (x,y) point in fGraph2D
106
107  // Initialise the Delaunay algorithm if needed.
108  // CreateTrianglesDataStructure computes fXoffset, fYoffset,
109  // fXScaleFactor and fYScaleFactor;
110  // needed in this function.
112
113  // Find the z value corresponding to the point (x,y).
114  double xx, yy;
115  xx = Linear_transform(x, fOffsetX, fScaleFactorX); //xx = xTransformer(x);
116  yy = Linear_transform(y, fOffsetY, fScaleFactorY); //yy = yTransformer(y);
117  double zz = DoInterpolateNormalized(xx, yy);
118
119  // Wrong zeros may appear when points sit on a regular grid.
120  // The following line try to avoid this problem.
121  if (zz==0) zz = DoInterpolateNormalized(xx+0.0001, yy);
122
123  return zz;
124 }
125
126 //______________________________________________________________________________
128 {
129
131  //treat the common case first
133  return;
134
135  Initialization cState = Initialization::UNINITIALIZED;
136  if(fInit.compare_exchange_strong(cState, Initialization::INITIALIZING,
137  std::memory_order::memory_order_release, std::memory_order::memory_order_relaxed))
138  {
139  // the value of fInit was indeed UNINIT, we replaced it atomically with initializing
140  // performing the initialzing now
141 #else
142  if (fInit) return; else fInit = kTRUE;
143 #endif
144
145  // Function used internally only. It creates the data structures needed to
146  // compute the Delaunay triangles.
147
148  // Offset fX and fY so they average zero, and scale so the average
149  // of the X and Y ranges is one. The normalized version of fX and fY used
150  // in Interpolate.
151
152  DoNormalizePoints(); // call backend specific point normalization
153
154  DoFindTriangles(); // call backend specific triangle finding
155
156  fNdt = fTriangles.size();
157
159  fInit = Initialization::INITIALIZED;
160  } else while(cState != Initialization::INITIALIZED) {
161  //the value of fInit was NOT UNINIT, so we have to spin until we reach INITIALEZED
163  }
164 #endif
165
166 }
167
168
169 // backend specific implementations
170
171 #ifdef HAS_CGAL
172 /// CGAL implementation of normalize points
173 void Delaunay2D::DonormalizePoints() {
174  for (Int_t n = 0; n < fNpoints; n++) {
175  //Point p(xTransformer(fX[n]), yTransformer(fY[n]));
176  Point p(linear_transform(fX[n], fOffsetX, fScaleFactorX),
177  linear_transform(fY[n], fOffsetY, fScaleFactorY));
178
179  fNormalizedPoints.insert(std::make_pair(p, n));
180  }
181 }
182
183 /// CGAL implementation for finding triangles
185  fCGALdelaunay.insert(fNormalizedPoints.begin(), fNormalizedPoints.end());
186
187  std::transform(fCGALdelaunay.finite_faces_begin(),
188  fCGALdelaunay.finite_faces_end(), std::back_inserter(fTriangles),
189  [] (const Delaunay::Face face) -> Triangle {
190
191  Triangle tri;
192
193  auto transform = [&] (const unsigned int i) {
194  tri.x[i] = face.vertex(i)->point().x();
195  tri.y[i] = face.vertex(i)->point().y();
196  tri.idx[i] = face.vertex(i)->info();
197  };
198
199  transform(0);
200  transform(1);
201  transform(2);
202
203  return tri;
204
205  });
206 }
207
208 /// CGAL implementation for interpolation
209 double Delaunay2D::DoInterpolateNormalized(double xx, double yy)
210 {
211  // Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
212  // calculate a z-value for it by linearly interpolating the z-values that
213  // make up that triangle.
214
215  // initialise the Delaunay algorithm if needed
217
218  //coordinate computation
219  Point p(xx, yy);
220
221  std::vector<std::pair<Point, Coord_type> > coords;
222  auto nn = CGAL::natural_neighbor_coordinates_2(fCGALdelaunay, p,
223  std::back_inserter(coords));
224
225  //std::cout << std::this_thread::get_id() << ": Found " << coords.size() << " points" << std::endl;
226
227  if(!nn.third) //neighbor finding was NOT successfull, return standard value
228  return fZout;
229
230  //printf("found neighbors %u\n", coords.size());
231
232  Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(),
233  nn.second, Value_access(fNormalizedPoints, fZ));
234
235  //std::cout << std::this_thread::get_id() << ": Result " << res << std::endl;
236
237  return res;
238 }
239
240 #else // HAS_CGAL
241
242 /// Triangle implementation for normalizing the points
244  for (Int_t n = 0; n < fNpoints; n++) {
246  fYN.push_back(Linear_transform(fY[n], fOffsetY, fScaleFactorY));
247  }
248
249  //also initialize fXCellStep and FYCellStep
252 }
253
254 /// Triangle implementation for finding all the triangles
256
257  auto initStruct = [] (triangulateio & s) {
258  s.pointlist = nullptr; /* In / out */
259  s.pointattributelist = nullptr; /* In / out */
260  s.pointmarkerlist = nullptr; /* In / out */
261  s.numberofpoints = 0; /* In / out */
262  s.numberofpointattributes = 0; /* In / out */
263
264  s.trianglelist = nullptr; /* In / out */
265  s.triangleattributelist = nullptr; /* In / out */
266  s.trianglearealist = nullptr; /* In only */
267  s.neighborlist = nullptr; /* Out only */
268  s.numberoftriangles = 0; /* In / out */
269  s.numberofcorners = 0; /* In / out */
270  s.numberoftriangleattributes = 0; /* In / out */
271
272  s.segmentlist = nullptr; /* In / out */
273  s.segmentmarkerlist = nullptr; /* In / out */
274  s.numberofsegments = 0; /* In / out */
275
276  s.holelist = nullptr; /* In / pointer to array copied out */
277  s.numberofholes = 0; /* In / copied out */
278
279  s.regionlist = nullptr; /* In / pointer to array copied out */
280  s.numberofregions = 0; /* In / copied out */
281
282  s.edgelist = nullptr; /* Out only */
283  s.edgemarkerlist = nullptr; /* Not used with Voronoi diagram; out only */
284  s.normlist = nullptr; /* Used only with Voronoi diagram; out only */
285  s.numberofedges = 0; /* Out only */
286  };
287
288  auto freeStruct = [] (triangulateio & s) {
289  if(s.pointlist != nullptr) free(s.pointlist); /* In / out */
290  if(s.pointattributelist != nullptr) free(s.pointattributelist); /* In / out */
291  if(s.pointmarkerlist != nullptr) free(s.pointmarkerlist); /* In / out */
292
293  if(s.trianglelist != nullptr) free(s.trianglelist); /* In / out */
294  if(s.triangleattributelist != nullptr) free(s.triangleattributelist); /* In / out */
295  if(s.trianglearealist != nullptr) free(s.trianglearealist); /* In only */
296  if(s.neighborlist != nullptr) free(s.neighborlist); /* Out only */
297
298  if(s.segmentlist != nullptr) free(s.segmentlist); /* In / out */
299  if(s.segmentmarkerlist != nullptr) free(s.segmentmarkerlist); /* In / out */
300
301  if(s.holelist != nullptr) free(s.holelist); /* In / pointer to array copied out */
302
303  if(s.regionlist != nullptr) free(s.regionlist); /* In / pointer to array copied out */
304
305  if(s.edgelist != nullptr) free(s.edgelist); /* Out only */
306  if(s.edgemarkerlist != nullptr) free(s.edgemarkerlist); /* Not used with Voronoi diagram; out only */
307  if(s.normlist != nullptr) free(s.normlist); /* Used only with Voronoi diagram; out only */
308  };
309
310  struct triangulateio in, out;
311  initStruct(in); initStruct(out);
312
313  /* Define input points. */
314
316  in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
317
318  for (Int_t i = 0; i < fNpoints; ++i) {
319  in.pointlist[2 * i] = fXN[i];
320  in.pointlist[2 * i + 1] = fYN[i];
321  }
322
323  triangulate((char *) "zQN", &in, &out, nullptr);
324
325  fTriangles.resize(out.numberoftriangles);
326  for(int t = 0; t < out.numberoftriangles; ++t){
327  Triangle tri;
328
329  auto transform = [&] (const unsigned int v) {
330  //each triangle as numberofcorners vertices ( = 3)
331  tri.idx[v] = out.trianglelist[t*out.numberofcorners + v];
332
333  //printf("triangle %u vertex %u: point %u/%i\n", t, v, tri.idx[v], out.numberofpoints);
334
335  //pointlist is [x0 y0 x1 y1 ...]
336  tri.x[v] = in.pointlist[tri.idx[v] * 2 + 0];
337
338  //printf("\t x: %f\n", tri.x[v]);
339
340  tri.y[v] = in.pointlist[tri.idx[v] * 2 + 1];
341
342  //printf("\t y: %f\n", tri.y[v]);
343  };
344
345  transform(0);
346  transform(1);
347  transform(2);
348
349  //see comment in header for CGAL fallback section
350  tri.invDenom = 1 / ( (tri.y[1] - tri.y[2])*(tri.x[0] - tri.x[2]) + (tri.x[2] - tri.x[1])*(tri.y[0] - tri.y[2]) );
351
352  fTriangles[t] = tri;
353
354  auto bx = std::minmax({tri.x[0], tri.x[1], tri.x[2]});
355  auto by = std::minmax({tri.y[0], tri.y[1], tri.y[2]});
356
357  unsigned int cellXmin = CellX(bx.first);
358  unsigned int cellXmax = CellX(bx.second);
359
360  unsigned int cellYmin = CellY(by.first);
361  unsigned int cellYmax = CellY(by.second);
362
363  for(unsigned int i = cellXmin; i <= cellXmax; ++i) {
364  for(unsigned int j = cellYmin; j <= cellYmax; ++j) {
365  //printf("(%u,%u) = %u\n", i, j, Cell(i,j));
366  fCells[Cell(i,j)].insert(t);
367  }
368  }
369  }
370
371  freeStruct(in); freeStruct(out);
372 }
373
374 /// Triangle implementation for interpolation
375 /// Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
376 /// calculate a z-value for it by linearly interpolating the z-values that
377 /// make up that triangle.
378 double Delaunay2D::DoInterpolateNormalized(double xx, double yy)
379 {
380
381  // relay that ll the triangles have been found
382  /// FindAllTriangles();
383
384  //see comment in header for CGAL fallback section
385  auto bayCoords = [&] (const unsigned int t) -> std::tuple<double, double, double> {
386  double la = ( (fTriangles[t].y[1] - fTriangles[t].y[2])*(xx - fTriangles[t].x[2])
387  + (fTriangles[t].x[2] - fTriangles[t].x[1])*(yy - fTriangles[t].y[2]) ) * fTriangles[t].invDenom;
388  double lb = ( (fTriangles[t].y[2] - fTriangles[t].y[0])*(xx - fTriangles[t].x[2])
389  + (fTriangles[t].x[0] - fTriangles[t].x[2])*(yy - fTriangles[t].y[2]) ) * fTriangles[t].invDenom;
390
391  return std::make_tuple(la, lb, (1 - la - lb));
392  };
393
394  auto inTriangle = [] (const std::tuple<double, double, double> & coords) -> bool {
395  return std::get<0>(coords) >= 0 && std::get<1>(coords) >= 0 && std::get<2>(coords) >= 0;
396  };
397
398  int cX = CellX(xx);
399  int cY = CellY(yy);
400
401  if(cX < 0 || cX > fNCells || cY < 0 || cY > fNCells)
402  return fZout; //TODO some more fancy interpolation here
403
404  for(unsigned int t : fCells[Cell(cX, cY)]){
405  auto coords = bayCoords(t);
406
407  if(inTriangle(coords)){
408  //we found the triangle -> interpolate using the barycentric interpolation
409  return std::get<0>(coords) * fZ[fTriangles[t].idx[0]]
410  + std::get<1>(coords) * fZ[fTriangles[t].idx[1]]
411  + std::get<2>(coords) * fZ[fTriangles[t].idx[2]];
412
413  }
414  }
415
416  //debugging
417
418  /*
419  for(unsigned int t = 0; t < fNdt; ++t) {
420  auto coords = bayCoords(t);
421
422  if(inTriangle(coords)){
423
424  //brute force found a triangle -> grid not
425  printf("Found triangle %u for (%f,%f) -> (%u,%u)\n", t, xx,yy, cX, cY);
426  printf("Triangles in grid cell: ");
427  for(unsigned int x : fCells[Cell(cX, cY)])
428  printf("%u ", x);
429  printf("\n");
430
431  printf("Triangle %u is in cells: ", t);
432  for(unsigned int i = 0; i <= fNCells; ++i)
433  for(unsigned int j = 0; j <= fNCells; ++j)
434  if(fCells[Cell(i,j)].count(t))
435  printf("(%u,%u) ", i, j);
436  printf("\n");
437  for(unsigned int i = 0; i < 3; ++i)
438  printf("\tpoint %u (%u): (%f,%f) -> (%u,%u)\n", i, fTriangles[t].idx[i], fTriangles[t].x[i], fTriangles[t].y[i], CellX(fTriangles[t].x[i]), CellY(fTriangles[t].y[i]));
439
440  //we found the triangle -> interpolate using the barycentric interpolation
441  return std::get<0>(coords) * fZ[fTriangles[t].idx[0]]
442  + std::get<1>(coords) * fZ[fTriangles[t].idx[1]]
443  + std::get<2>(coords) * fZ[fTriangles[t].idx[2]];
444
445  }
446  }
447
448  printf("Could not find a triangle for point (%f,%f)\n", xx, yy);
449  */
450
451  //no triangle found return standard value
452  return fZout;
453 }
454
455 } // namespace Math
456 } // namespace ROOT
457
458
459 #endif //HAS_CGAL
unsigned int Cell(UInt_t x, UInt_t y) const
grid cells with containing triangles
Definition: Delaunay2D.h:276
double fZout
Normalization factor Y.
Definition: Delaunay2D.h:173
void DoNormalizePoints()
internal function to normalize the points
Definition: Delaunay2D.cxx:243
float xmin
Definition: THbookFile.cxx:93
This namespace contains pre-defined functions to be used in conjuction with TExecutor::Map and TExecu...
Definition: StringConv.hxx:21
int * trianglelist
Definition: triangle.h:295
float ymin
Definition: THbookFile.cxx:93
Bool_t fInit
Height for points lying outside the convex hull.
Definition: Delaunay2D.h:181
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: Rtypes.h:92
void FindAllTriangles()
Find all triangles.
Definition: Delaunay2D.cxx:127
int CellY(double y) const
Definition: Delaunay2D.h:284
#define malloc
Definition: civetweb.c:818
double fScaleFactorX
Normalization offset Y.
Definition: Delaunay2D.h:170
#define REAL
Definition: triangle.h:277
double DoInterpolateNormalized(double x, double y)
internal method to compute the interpolation
Definition: Delaunay2D.cxx:378
Double_t x[n]
Definition: legend1.C:17
double fOffsetY
Normalization offset X.
Definition: Delaunay2D.h:168
double fScaleFactorY
Normalization factor X.
Definition: Delaunay2D.h:171
double fOffsetX
Maximum value of fYN.
Definition: Delaunay2D.h:167
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:58
std::set< UInt_t > fCells[(fNCells+1) *(fNCells+1)]
inverse denominator to calculate X cell = fNCells / (fYNmax - fYNmin)
Definition: Delaunay2D.h:274
double fXNmin
Pointer to Z array.
Definition: Delaunay2D.h:159
float ymax
Definition: THbookFile.cxx:93
SVector< double, 2 > v
Definition: Dict.h:5
const double * fY
Pointer to X array (managed externally)
Definition: Delaunay2D.h:156
double fXNmax
Minimum value of fXN.
Definition: Delaunay2D.h:160
int numberofcorners
Definition: triangle.h:300
std::vector< double > fYN
normalized X
Definition: Delaunay2D.h:264
float xmax
Definition: THbookFile.cxx:93
double fYNmax
Minimum value of fYN.
Definition: Delaunay2D.h:162
double fXCellStep
number of cells to divide the normalized space
Definition: Delaunay2D.h:272
int numberofpoints
Definition: triangle.h:292
Int_t fNpoints
Number of Delaunay triangles found.
Definition: Delaunay2D.h:153
double fYCellStep
inverse denominator to calculate X cell = fNCells / (fXNmax - fXNmin)
Definition: Delaunay2D.h:273
#define free
Definition: civetweb.c:821
Double_t y[n]
Definition: legend1.C:17
Namespace for new Math classes and functions.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
static const Int_t UNINITIALIZED
static const int fNCells
normalized Y
Definition: Delaunay2D.h:271
void DoFindTriangles()
internal function to find the triangle use Triangle or CGAL if flag is set
Definition: Delaunay2D.cxx:255
std::vector< double > fXN
Triangles of Triangulation.
Definition: Delaunay2D.h:263
void triangulate(char *triswitches, struct triangulateio *in, struct triangulateio *out, struct triangulateio *vorout)
Definition: triangle.c:15651
const double * fX
Number of data points.
Definition: Delaunay2D.h:155
const int nn
double fYNmin
Maximum value of fXN.
Definition: Delaunay2D.h:161
const Bool_t kTRUE
Definition: Rtypes.h:91
int CellX(double x) const
Definition: Delaunay2D.h:280
REAL * pointlist
Definition: triangle.h:289
double Interpolate(double x, double y)
Return the Interpolated z value corresponding to the (x,y) point.
Definition: Delaunay2D.cxx:103
const Int_t n
Definition: legend1.C:16
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:127
const double * fZ
Pointer to Y array.
Definition: Delaunay2D.h:157
Triangles fTriangles
True if FindAllTriangels() has been performed.
Definition: Delaunay2D.h:185
int numberoftriangles
Definition: triangle.h:299