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