Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
28namespace ROOT {
29
30 namespace Math {
31
32
33/// class constructor from array of data points
34Delaunay2D::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
55
56}
57
58/// set the input points
59void 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
64 fInit = Initialization::UNINITIALIZED;
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//______________________________________________________________________________
104double 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
174void 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
210double 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
217 FindAllTriangles();
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
244void Delaunay2D::DoNormalizePoints() {
245 for (Int_t n = 0; n < fNpoints; n++) {
246 fXN.push_back(Linear_transform(fX[n], fOffsetX, fScaleFactorX));
247 fYN.push_back(Linear_transform(fY[n], fOffsetY, fScaleFactorY));
248 }
249
250 //also initialize fXCellStep and FYCellStep
251 fXCellStep = fNCells / (fXNmax - fXNmin);
252 fYCellStep = fNCells / (fYNmax - fYNmin);
253}
254
255/// Triangle implementation for finding all the triangles
256void Delaunay2D::DoFindTriangles() {
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
316 in.numberofpoints = fNpoints;
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.
379double 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#endif //HAS_CGAL
456
457} // namespace Math
458} // namespace ROOT
459
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:101
const Bool_t kTRUE
Definition RtypesCore.h:100
float xmin
float ymin
float xmax
float ymax
#define free
Definition civetweb.c:1539
#define malloc
Definition civetweb.c:1536
double fXNmin
Pointer to Z array.
Definition Delaunay2D.h:157
double DoInterpolateNormalized(double x, double y)
internal method to compute the interpolation
double Interpolate(double x, double y)
Return the Interpolated z value corresponding to the (x,y) point.
const double * fZ
Pointer to Y array.
Definition Delaunay2D.h:155
double fXNmax
Minimum value of fXN.
Definition Delaunay2D.h:158
double fZout
Normalization factor Y.
Definition Delaunay2D.h:171
double fOffsetY
Normalization offset X.
Definition Delaunay2D.h:166
Triangles fTriangles
True if FindAllTriangels() has been performed.
Definition Delaunay2D.h:183
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
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
void FindAllTriangles()
Find all triangles.
double fXCellStep
number of cells to divide the normalized space
Definition Delaunay2D.h:270
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
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
const double * fX
Number of data points.
Definition Delaunay2D.h:153
void DoNormalizePoints()
internal function to normalize the points
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
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.
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
int numberofpoints
Definition triangle.h:292
int numberoftriangles
Definition triangle.h:299
int * trianglelist
Definition triangle.h:295
REAL * pointlist
Definition triangle.h:289
int numberofcorners
Definition triangle.h:300
void triangulate(char *triswitches, struct triangulateio *in, struct triangulateio *out, struct triangulateio *vorout)
Definition triangle.c:15653
#define REAL
Definition triangle.h:277