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