ROOT   Reference Guide
Searching...
No Matches
Delaunay2D.cxx
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// Implementation file for class Delaunay2D
13
14#include "Math/Delaunay2D.h"
15#include "Rtypes.h"
16
18
19// use the CDT library if we do not use CGAL
20#ifndef HAS_CGAL
21#include "CDT/CDT.h"
22#endif
23
24#include <algorithm>
25#include <cstdlib>
26
27#include <iostream>
28#include <limits>
29
30
31namespace ROOT {
32
33 namespace Math {
34
35
36/// class constructor from array of data points
37Delaunay2D::Delaunay2D(int n, const double * x, const double * y, const double * z,
38 double xmin, double xmax, double ymin, double ymax)
39{
40 // Delaunay2D normal constructor
41
42 fX = x;
43 fY = y;
44 fZ = z;
45 fZout = 0.;
46 fNpoints = n;
47 fOffsetX = 0;
48 fOffsetY = 0;
49 fScaleFactorX = 0;
50 fScaleFactorY = 0;
51 fNdt = 0;
52 fXNmax = 0;
53 fXNmin = 0;
54 fYNmin = 0;
55 fYNmax = 0;
56
58
59}
60
61/// set the input points
62void Delaunay2D::SetInputPoints(int n, const double * x, const double * y, const double * z,
63 double xmin, double xmax, double ymin, double ymax) {
64
65
66 fInit = kFALSE;
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 if ( (xmax-xmin) != 0 ) {
82 fScaleFactorX = 1./(xmax-xmin);
83 fXNmax = Linear_transform(xmax, fOffsetX, fScaleFactorX); //xTransformer(xmax);
84 fXNmin = Linear_transform(xmin, fOffsetX, fScaleFactorX); //xTransformer(xmin);
85 } else {
86 // we can't find triangle in this case
87 fInit = true;
88 }
89
90 if (ymax-ymin != 0) {
91 fScaleFactorY = 1./(ymax-ymin);
92 fYNmax = Linear_transform(ymax, fOffsetY, fScaleFactorY); //yTransformer(ymax);
93 fYNmin = Linear_transform(ymin, fOffsetY, fScaleFactorY); //yTransformer(ymin);
94 } else {
95 fInit = true;
96 }
97
98#ifndef HAS_CGAL
99 fXCellStep = 0.;
100 fYCellStep = 0.;
101#endif
102}
103
104//______________________________________________________________________________
105double Delaunay2D::Interpolate(double x, double y)
106{
107 // Return the interpolated z value corresponding to the given (x,y) point
108
109 // Initialise the Delaunay algorithm if needed.
110 // CreateTrianglesDataStructure computes fXoffset, fYoffset,
111 // fXScaleFactor and fYScaleFactor;
112 // needed in this function.
114
115 // handle case no triangles are found return default value (i.e. 0)
116 // to do: if points are aligned in one direction we could return in
117 // some case the 1d interpolated value
118 if (fNdt == 0) {
119 return fZout;
120 }
121
122 // Find the z value corresponding to the point (x,y).
123 double xx, yy;
124 xx = Linear_transform(x, fOffsetX, fScaleFactorX); //xx = xTransformer(x);
125 yy = Linear_transform(y, fOffsetY, fScaleFactorY); //yy = yTransformer(y);
126 double zz = DoInterpolateNormalized(xx, yy);
127
128 // the case of points on a regular grid (i.e. points on triangle edges) it is now handles in
129 // DoInterpolateNormalized
130
131 return zz;
132}
133
134//______________________________________________________________________________
136{
137
138 if (fInit)
139 return;
140 else
141 fInit = kTRUE;
142
143 // Function used internally only. It creates the data structures needed to
144 // compute the Delaunay triangles.
145
146 // Offset fX and fY so they average zero, and scale so the average
147 // of the X and Y ranges is one. The normalized version of fX and fY used
148 // in Interpolate.
149
150 DoNormalizePoints(); // call backend specific point normalization
151
152 DoFindTriangles(); // call backend specific triangle finding
153
154 fNdt = fTriangles.size();
155}
156
157// backend specific implementations
158
159#ifndef HAS_CGAL
160
161// Triangle implementation (default case)
162
163/// Triangle implementation for points normalization
165 for (Int_t n = 0; n < fNpoints; n++) {
168 }
169
170 //also initialize fXCellStep and FYCellStep
173}
174
175/// Triangle implementation for finding all the triangles
177
178 int i;
179 std::vector<CDT::V2d<double>> points(fNpoints);
180 for (i = 0; i < fNpoints; ++i) points[i] = CDT::V2d<double>::make(fXN[i], fYN[i]);
181 CDT::RemoveDuplicates(points);
182
183 CDT::Triangulation<double> cdt;
184 cdt.insertVertices(points);
185 cdt.eraseSuperTriangle();
186
187 auto AllTriangles = cdt.triangles;
188 auto AllVertices = cdt.vertices;
189 int NumberOfTriangles = cdt.triangles.size();
190
192
193 for(i = 0; i < NumberOfTriangles; i++){
194 Triangle tri;
195 const auto& t = AllTriangles[i];
196
197 const auto& v0 = AllVertices[t.vertices[0]];
198 tri.x[0] = v0.x;
199 tri.y[0] = v0.y;
200 tri.idx[0] = t.vertices[0];
201
202 const auto& v1 = AllVertices[t.vertices[1]];
203 tri.x[1] = v1.x;
204 tri.y[1] = v1.y;
205 tri.idx[1] = t.vertices[1];
206
207 const auto& v2 = AllVertices[t.vertices[2]];
208 tri.x[2] = v2.x;
209 tri.y[2] = v2.y;
210 tri.idx[2] = t.vertices[2];
211
212 // see comment in header for CGAL fallback section
213 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]) );
214
215 fTriangles[i] = tri;
216
217 auto bx = std::minmax({tri.x[0], tri.x[1], tri.x[2]});
218 auto by = std::minmax({tri.y[0], tri.y[1], tri.y[2]});
219
220 unsigned int cellXmin = CellX(bx.first);
221 unsigned int cellXmax = CellX(bx.second);
222
223 unsigned int cellYmin = CellY(by.first);
224 unsigned int cellYmax = CellY(by.second);
225
226 for(unsigned int j = cellXmin; j <= cellXmax; j++) {
227 for(unsigned int k = cellYmin; k <= cellYmax; k++) {
228 fCells[Cell(j,k)].insert(i);
229 }
230 }
231 }
232}
233
234/// Triangle implementation for interpolation
235/// Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
236/// calculate a z-value for it by linearly interpolating the z-values that
237/// make up that triangle.
238/// Relay that all the triangles have been found before
239/// see comment in class description (in Delaunay2D.h) for implementation details:
240/// finding barycentric coordinates and computing the interpolation
241double Delaunay2D::DoInterpolateNormalized(double xx, double yy)
242{
243
244 // compute barycentric coordinates of a point P(xx,yy,zz)
245 auto bayCoords = [&](const unsigned int t) -> std::tuple<double, double, double> {
246 double la = ((fTriangles[t].y[1] - fTriangles[t].y[2]) * (xx - fTriangles[t].x[2]) +
247 (fTriangles[t].x[2] - fTriangles[t].x[1]) * (yy - fTriangles[t].y[2])) *
248 fTriangles[t].invDenom;
249 double lb = ((fTriangles[t].y[2] - fTriangles[t].y[0]) * (xx - fTriangles[t].x[2]) +
250 (fTriangles[t].x[0] - fTriangles[t].x[2]) * (yy - fTriangles[t].y[2])) *
251 fTriangles[t].invDenom;
252
253 return std::make_tuple(la, lb, (1 - la - lb));
254 };
255
256 // function to test if a point with barycentric coordinates (a,b,c) is inside the triangle
257 // If the point is outside one or more of the coordinate are negative.
258 // If the point is on a triangle edge, one of the coordinate (the one not part of the edge) is zero.
259 // Due to numerical error, it can happen that if the point is at the edge the result is a small negative value.
260 // Use then a tolerance (of - eps) to still consider the point within the triangle
261 auto inTriangle = [](const std::tuple<double, double, double> &coords) -> bool {
262 constexpr double eps = -4 * std::numeric_limits<double>::epsilon();
263 return std::get<0>(coords) > eps && std::get<1>(coords) > eps && std::get<2>(coords) > eps;
264 };
265
266 int cX = CellX(xx);
267 int cY = CellY(yy);
268
269 if (cX < 0 || cX > fNCells || cY < 0 || cY > fNCells)
270 return fZout; // TODO some more fancy interpolation here
271
272 for (unsigned int t : fCells[Cell(cX, cY)]) {
273
274 auto coords = bayCoords(t);
275
276 // std::cout << "result of bayCoords " << std::get<0>(coords) <<
277 // " " << std::get<1>(coords) << " " << std::get<2>(coords) << std::endl;
278
279 if (inTriangle(coords)) {
280
281 // we found the triangle -> interpolate using the barycentric interpolation
282
283 return std::get<0>(coords) * fZ[fTriangles[t].idx[0]] + std::get<1>(coords) * fZ[fTriangles[t].idx[1]] +
284 std::get<2>(coords) * fZ[fTriangles[t].idx[2]];
285 }
286 }
287
288 // no triangle found return standard value
289 return fZout;
290}
291
292#else //HAS_CGAL: case of using GCAL
293
294/// CGAL implementation of normalize points
295void Delaunay2D::DonormalizePoints() {
296 for (Int_t n = 0; n < fNpoints; n++) {
297 //Point p(xTransformer(fX[n]), yTransformer(fY[n]));
298 Point p(linear_transform(fX[n], fOffsetX, fScaleFactorX),
299 linear_transform(fY[n], fOffsetY, fScaleFactorY));
300
301 fNormalizedPoints.insert(std::make_pair(p, n));
302 }
303}
304
305/// CGAL implementation for finding triangles
307 fCGALdelaunay.insert(fNormalizedPoints.begin(), fNormalizedPoints.end());
308
309 std::transform(fCGALdelaunay.finite_faces_begin(),
310 fCGALdelaunay.finite_faces_end(), std::back_inserter(fTriangles),
311 [] (const Delaunay::Face face) -> Triangle {
312
313 Triangle tri;
314
315 auto transform = [&] (const unsigned int i) {
316 tri.x[i] = face.vertex(i)->point().x();
317 tri.y[i] = face.vertex(i)->point().y();
318 tri.idx[i] = face.vertex(i)->info();
319 };
320
321 transform(0);
322 transform(1);
323 transform(2);
324
325 return tri;
326
327 });
328}
329
330/// CGAL implementation for interpolation
331double Delaunay2D::DoInterpolateNormalized(double xx, double yy)
332{
333 // Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
334 // calculate a z-value for it by linearly interpolating the z-values that
335 // make up that triangle.
336
337 // initialise the Delaunay algorithm if needed
338 FindAllTriangles();
339
340 //coordinate computation
341 Point p(xx, yy);
342
343 std::vector<std::pair<Point, Coord_type> > coords;
344 auto nn = CGAL::natural_neighbor_coordinates_2(fCGALdelaunay, p,
345 std::back_inserter(coords));
346
347 //std::cout << std::this_thread::get_id() << ": Found " << coords.size() << " points" << std::endl;
348
349 if(!nn.third) // neighbour finding was NOT successful, return standard value
350 return fZout;
351
352 Coord_type res = CGAL::linear_interpolation(coords.begin(), coords.end(),
353 nn.second, Value_access(fNormalizedPoints, fZ));
354
355 //std::cout << std::this_thread::get_id() << ": Result " << res << std::endl;
356
357 return res;
358}
359#endif // HAS_GCAL
360
361} // namespace Math
362} // namespace ROOT
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
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 points
float xmin
float ymin
float xmax
float ymax
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
double Interpolate(double x, double y)
Return the Interpolated z value corresponding to the given (x,y) point.
const double * fZ
! Pointer to Z array
Definition Delaunay2D.h:199
double fXNmax
! Maximum value of fXN
Definition Delaunay2D.h:202
double fZout
! Height for points lying outside the convex hull
Definition Delaunay2D.h:212
double fOffsetY
! Normalization offset Y
Definition Delaunay2D.h:207
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
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.
static const int fNCells
! number of cells to divide the normalized space
Definition Delaunay2D.h:226
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
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
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...