Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGraphDelaunay.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id: TGraphDelaunay.cxx,v 1.00
2// Author: Olivier Couet, Luke Jones (Royal Holloway, University of London)
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
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#include "TMath.h"
13#include "TGraph2D.h"
14#include "TGraphDelaunay.h"
15
17
18
19/** \class TGraphDelaunay
20 \ingroup Hist
21TGraphDelaunay generates a Delaunay triangulation of a TGraph2D. This
22triangulation code derives from an implementation done by Luke Jones
23(Royal Holloway, University of London) in April 2002 in the PAW context.
24
25This software cannot be guaranteed to work under all circumstances. They
26were originally written to work with a few hundred points in an XY space
27with similar X and Y ranges.
28
29Definition of Delaunay triangulation (After B. Delaunay):
30For a set S of points in the Euclidean plane, the unique triangulation DT(S)
31of S such that no point in S is inside the circumcircle of any triangle in
32DT(S). DT(S) is the dual of the Voronoi diagram of S. If n is the number of
33points in S, the Voronoi diagram of S is the partitioning of the plane
34containing S points into n convex polygons such that each polygon contains
35exactly one point and every point in a given polygon is closer to its
36central point than to any other. A Voronoi diagram is sometimes also known
37as a Dirichlet tessellation.
38
39\image html tgraph2d_delaunay.png
40
41[This applet](http://www.cs.cornell.edu/Info/People/chew/Delaunay.html) gives a nice
42practical view of Delaunay triangulation and Voronoi diagram.
43*/
44
45
46////////////////////////////////////////////////////////////////////////////////
47/// TGraphDelaunay default constructor
48
50 : TNamed("TGraphDelaunay","TGraphDelaunay")
51{
52 fGraph2D = 0;
53 fX = 0;
54 fY = 0;
55 fZ = 0;
56 fNpoints = 0;
57 fTriedSize = 0;
58 fZout = 0.;
59 fNdt = 0;
60 fNhull = 0;
61 fHullPoints = 0;
62 fXN = 0;
63 fYN = 0;
64 fOrder = 0;
65 fDist = 0;
66 fPTried = 0;
67 fNTried = 0;
68 fMTried = 0;
69 fInit = kFALSE;
70 fXNmin = 0.;
71 fXNmax = 0.;
72 fYNmin = 0.;
73 fYNmax = 0.;
74 fXoffset = 0.;
75 fYoffset = 0.;
76 fXScaleFactor = 0.;
77 fYScaleFactor = 0.;
78
79 SetMaxIter();
80}
81
82
83////////////////////////////////////////////////////////////////////////////////
84/// TGraphDelaunay normal constructor
85
87 : TNamed("TGraphDelaunay","TGraphDelaunay")
88{
89 fGraph2D = g;
90 fX = fGraph2D->GetX();
91 fY = fGraph2D->GetY();
92 fZ = fGraph2D->GetZ();
94 fTriedSize = 0;
95 fZout = 0.;
96 fNdt = 0;
97 fNhull = 0;
98 fHullPoints = 0;
99 fXN = 0;
100 fYN = 0;
101 fOrder = 0;
102 fDist = 0;
103 fPTried = 0;
104 fNTried = 0;
105 fMTried = 0;
106 fInit = kFALSE;
107 fXNmin = 0.;
108 fXNmax = 0.;
109 fYNmin = 0.;
110 fYNmax = 0.;
111 fXoffset = 0.;
112 fYoffset = 0.;
113 fXScaleFactor = 0.;
114 fYScaleFactor = 0.;
115
116 SetMaxIter();
117}
118
119
120////////////////////////////////////////////////////////////////////////////////
121/// TGraphDelaunay destructor.
122
124{
125 if (fPTried) delete [] fPTried;
126 if (fNTried) delete [] fNTried;
127 if (fMTried) delete [] fMTried;
128 if (fHullPoints) delete [] fHullPoints;
129 if (fOrder) delete [] fOrder;
130 if (fDist) delete [] fDist;
131 if (fXN) delete [] fXN;
132 if (fYN) delete [] fYN;
133
134 fPTried = 0;
135 fNTried = 0;
136 fMTried = 0;
137 fHullPoints = 0;
138 fOrder = 0;
139 fDist = 0;
140 fXN = 0;
141 fYN = 0;
142}
143
144
145////////////////////////////////////////////////////////////////////////////////
146/// Return the z value corresponding to the (x,y) point in fGraph2D
147
149{
150 // Initialise the Delaunay algorithm if needed.
151 // CreateTrianglesDataStructure computes fXoffset, fYoffset,
152 // fXScaleFactor and fYScaleFactor;
153 // needed in this function.
154 if (!fInit) {
156 FindHull();
157 fInit = kTRUE;
158 }
159
160 // Find the z value corresponding to the point (x,y).
161 Double_t xx, yy;
162 xx = (x+fXoffset)*fXScaleFactor;
163 yy = (y+fYoffset)*fYScaleFactor;
164 Double_t zz = Interpolate(xx, yy);
165
166 // Wrong zeros may appear when points sit on a regular grid.
167 // The following line try to avoid this problem.
168 if (zz==0) zz = Interpolate(xx+0.0001, yy);
169
170 return zz;
171}
172
173
174////////////////////////////////////////////////////////////////////////////////
175/// Function used internally only. It creates the data structures needed to
176/// compute the Delaunay triangles.
177
179{
180 // Offset fX and fY so they average zero, and scale so the average
181 // of the X and Y ranges is one. The normalized version of fX and fY used
182 // in Interpolate.
187 fXoffset = -(xmax+xmin)/2.;
188 fYoffset = -(ymax+ymin)/2.;
189 fXScaleFactor = 1./(xmax-xmin);
190 fYScaleFactor = 1./(ymax-ymin);
195 fXN = new Double_t[fNpoints+1];
196 fYN = new Double_t[fNpoints+1];
197 for (Int_t n=0; n<fNpoints; n++) {
198 fXN[n+1] = (fX[n]+fXoffset)*fXScaleFactor;
199 fYN[n+1] = (fY[n]+fYoffset)*fYScaleFactor;
200 }
201
202 // If needed, creates the arrays to hold the Delaunay triangles.
203 // A maximum number of 2*fNpoints is guessed. If more triangles will be
204 // find, FillIt will automatically enlarge these arrays.
206 fPTried = new Int_t[fTriedSize];
207 fNTried = new Int_t[fTriedSize];
208 fMTried = new Int_t[fTriedSize];
209}
210
211
212////////////////////////////////////////////////////////////////////////////////
213/// Is point e inside the triangle t1-t2-t3 ?
214
216{
217 Double_t x[4],y[4],xp, yp;
218 x[0] = fXN[t1];
219 x[1] = fXN[t2];
220 x[2] = fXN[t3];
221 x[3] = x[0];
222 y[0] = fYN[t1];
223 y[1] = fYN[t2];
224 y[2] = fYN[t3];
225 y[3] = y[0];
226 xp = fXN[e];
227 yp = fYN[e];
228 return TMath::IsInside(xp, yp, 4, x, y);
229}
230
231
232////////////////////////////////////////////////////////////////////////////////
233/// Files the triangle defined by the 3 vertices p, n and m into the
234/// fxTried arrays. If these arrays are to small they are automatically
235/// expanded.
236
238{
239 Bool_t swap;
240 Int_t tmp, ps = p, ns = n, ms = m;
241
242 // order the vertices before storing them
243L1:
244 swap = kFALSE;
245 if (ns > ps) { tmp = ps; ps = ns; ns = tmp; swap = kTRUE;}
246 if (ms > ns) { tmp = ns; ns = ms; ms = tmp; swap = kTRUE;}
247 if (swap) goto L1;
248
249 // expand the triangles storage if needed
250 if (fNdt>=fTriedSize) {
251 Int_t newN = 2*fTriedSize;
252 Int_t *savep = new Int_t [newN];
253 Int_t *saven = new Int_t [newN];
254 Int_t *savem = new Int_t [newN];
255 memcpy(savep,fPTried,fTriedSize*sizeof(Int_t));
256 memset(&savep[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
257 delete [] fPTried;
258 memcpy(saven,fNTried,fTriedSize*sizeof(Int_t));
259 memset(&saven[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
260 delete [] fNTried;
261 memcpy(savem,fMTried,fTriedSize*sizeof(Int_t));
262 memset(&savem[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
263 delete [] fMTried;
264 fPTried = savep;
265 fNTried = saven;
266 fMTried = savem;
267 fTriedSize = newN;
268 }
269
270 // store a new Delaunay triangle
271 fNdt++;
272 fPTried[fNdt-1] = ps;
273 fNTried[fNdt-1] = ns;
274 fMTried[fNdt-1] = ms;
275}
276
277
278////////////////////////////////////////////////////////////////////////////////
279/// Attempt to find all the Delaunay triangles of the point set. It is not
280/// guaranteed that it will fully succeed, and no check is made that it has
281/// fully succeeded (such a check would be possible by referencing the points
282/// that make up the convex hull). The method is to check if each triangle
283/// shares all three of its sides with other triangles. If not, a point is
284/// generated just outside the triangle on the side(s) not shared, and a new
285/// triangle is found for that point. If this method is not working properly
286/// (many triangles are not being found) it's probably because the new points
287/// are too far beyond or too close to the non-shared sides. Fiddling with
288/// the size of the `alittlebit' parameter may help.
289
291{
292 if (fAllTri) return; else fAllTri = kTRUE;
293
294 Double_t xcntr,ycntr,xm,ym,xx,yy;
295 Double_t sx,sy,nx,ny,mx,my,mdotn,nn,a;
296 Int_t t1,t2,pa,na,ma,pb,nb,mb,p1=0,p2=0,m,n,p3=0;
297 Bool_t s[3];
298 Double_t alittlebit = 0.0001;
299
300 // start with a point that is guaranteed to be inside the hull (the
301 // centre of the hull). The starting point is shifted "a little bit"
302 // otherwise, in case of triangles aligned on a regular grid, we may
303 // found none of them.
304 xcntr = 0;
305 ycntr = 0;
306 for (n=1; n<=fNhull; n++) {
307 xcntr = xcntr+fXN[fHullPoints[n-1]];
308 ycntr = ycntr+fYN[fHullPoints[n-1]];
309 }
310 xcntr = xcntr/fNhull+alittlebit;
311 ycntr = ycntr/fNhull+alittlebit;
312 // and calculate it's triangle
313 Interpolate(xcntr,ycntr);
314
315 // loop over all Delaunay triangles (including those constantly being
316 // produced within the loop) and check to see if their 3 sides also
317 // correspond to the sides of other Delaunay triangles, i.e. that they
318 // have all their neighbours.
319 t1 = 1;
320 while (t1 <= fNdt) {
321 // get the three points that make up this triangle
322 pa = fPTried[t1-1];
323 na = fNTried[t1-1];
324 ma = fMTried[t1-1];
325
326 // produce three integers which will represent the three sides
327 s[0] = kFALSE;
328 s[1] = kFALSE;
329 s[2] = kFALSE;
330 // loop over all other Delaunay triangles
331 for (t2=1; t2<=fNdt; t2++) {
332 if (t2 != t1) {
333 // get the points that make up this triangle
334 pb = fPTried[t2-1];
335 nb = fNTried[t2-1];
336 mb = fMTried[t2-1];
337 // do triangles t1 and t2 share a side?
338 if ((pa==pb && na==nb) || (pa==pb && na==mb) || (pa==nb && na==mb)) {
339 // they share side 1
340 s[0] = kTRUE;
341 } else if ((pa==pb && ma==nb) || (pa==pb && ma==mb) || (pa==nb && ma==mb)) {
342 // they share side 2
343 s[1] = kTRUE;
344 } else if ((na==pb && ma==nb) || (na==pb && ma==mb) || (na==nb && ma==mb)) {
345 // they share side 3
346 s[2] = kTRUE;
347 }
348 }
349 // if t1 shares all its sides with other Delaunay triangles then
350 // forget about it
351 if (s[0] && s[1] && s[2]) continue;
352 }
353 // Looks like t1 is missing a neighbour on at least one side.
354 // For each side, take a point a little bit beyond it and calculate
355 // the Delaunay triangle for that point, this should be the triangle
356 // which shares the side.
357 for (m=1; m<=3; m++) {
358 if (!s[m-1]) {
359 // get the two points that make up this side
360 if (m == 1) {
361 p1 = pa;
362 p2 = na;
363 p3 = ma;
364 } else if (m == 2) {
365 p1 = pa;
366 p2 = ma;
367 p3 = na;
368 } else if (m == 3) {
369 p1 = na;
370 p2 = ma;
371 p3 = pa;
372 }
373 // get the coordinates of the centre of this side
374 xm = (fXN[p1]+fXN[p2])/2.;
375 ym = (fYN[p1]+fYN[p2])/2.;
376 // we want to add a little to these coordinates to get a point just
377 // outside the triangle; (sx,sy) will be the vector that represents
378 // the side
379 sx = fXN[p1]-fXN[p2];
380 sy = fYN[p1]-fYN[p2];
381 // (nx,ny) will be the normal to the side, but don't know if it's
382 // pointing in or out yet
383 nx = sy;
384 ny = -sx;
385 nn = TMath::Sqrt(nx*nx+ny*ny);
386 nx = nx/nn;
387 ny = ny/nn;
388 mx = fXN[p3]-xm;
389 my = fYN[p3]-ym;
390 mdotn = mx*nx+my*ny;
391 if (mdotn > 0) {
392 // (nx,ny) is pointing in, we want it pointing out
393 nx = -nx;
394 ny = -ny;
395 }
396 // increase/decrease xm and ym a little to produce a point
397 // just outside the triangle (ensuring that the amount added will
398 // be large enough such that it won't be lost in rounding errors)
399 a = TMath::Abs(TMath::Max(alittlebit*xm,alittlebit*ym));
400 xx = xm+nx*a;
401 yy = ym+ny*a;
402 // try and find a new Delaunay triangle for this point
403 Interpolate(xx,yy);
404
405 // this side of t1 should now, hopefully, if it's not part of the
406 // hull, be shared with a new Delaunay triangle just calculated by Interpolate
407 }
408 }
409 t1++;
410 }
411}
412
413
414////////////////////////////////////////////////////////////////////////////////
415/// Finds those points which make up the convex hull of the set. If the xy
416/// plane were a sheet of wood, and the points were nails hammered into it
417/// at the respective coordinates, then if an elastic band were stretched
418/// over all the nails it would form the shape of the convex hull. Those
419/// nails in contact with it are the points that make up the hull.
420
422{
423 Int_t n,nhull_tmp;
424 Bool_t in;
425
427
428 nhull_tmp = 0;
429 for(n=1; n<=fNpoints; n++) {
430 // if the point is not inside the hull of the set of all points
431 // bar it, then it is part of the hull of the set of all points
432 // including it
433 in = InHull(n,n);
434 if (!in) {
435 // cannot increment fNhull directly - InHull needs to know that
436 // the hull has not yet been completely found
437 nhull_tmp++;
438 fHullPoints[nhull_tmp-1] = n;
439 }
440 }
441 fNhull = nhull_tmp;
442}
443
444
445////////////////////////////////////////////////////////////////////////////////
446/// Is point e inside the hull defined by all points apart from x ?
447
449{
450 Int_t n1,n2,n,m,ntry;
451 Double_t lastdphi,dd1,dd2,dx1,dx2,dx3,dy1,dy2,dy3;
452 Double_t u,v,vNv1,vNv2,phi1,phi2,dphi,xx,yy;
453
454 Bool_t deTinhull = kFALSE;
455
456 xx = fXN[e];
457 yy = fYN[e];
458
459 if (fNhull > 0) {
460 // The hull has been found - no need to use any points other than
461 // those that make up the hull
462 ntry = fNhull;
463 } else {
464 // The hull has not yet been found, will have to try every point
465 ntry = fNpoints;
466 }
467
468 // n1 and n2 will represent the two points most separated by angle
469 // from point e. Initially the angle between them will be <180 degs.
470 // But subsequent points will increase the n1-e-n2 angle. If it
471 // increases above 180 degrees then point e must be surrounded by
472 // points - it is not part of the hull.
473 n1 = 1;
474 n2 = 2;
475 if (n1 == x) {
476 n1 = n2;
477 n2++;
478 } else if (n2 == x) {
479 n2++;
480 }
481
482 // Get the angle n1-e-n2 and set it to lastdphi
483 dx1 = xx-fXN[n1];
484 dy1 = yy-fYN[n1];
485 dx2 = xx-fXN[n2];
486 dy2 = yy-fYN[n2];
487 phi1 = TMath::ATan2(dy1,dx1);
488 phi2 = TMath::ATan2(dy2,dx2);
489 dphi = (phi1-phi2)-((Int_t)((phi1-phi2)/TMath::TwoPi())*TMath::TwoPi());
490 if (dphi < 0) dphi = dphi+TMath::TwoPi();
491 lastdphi = dphi;
492 for (n=1; n<=ntry; n++) {
493 if (fNhull > 0) {
494 // Try hull point n
495 m = fHullPoints[n-1];
496 } else {
497 m = n;
498 }
499 if ((m!=n1) && (m!=n2) && (m!=x)) {
500 // Can the vector e->m be represented as a sum with positive
501 // coefficients of vectors e->n1 and e->n2?
502 dx1 = xx-fXN[n1];
503 dy1 = yy-fYN[n1];
504 dx2 = xx-fXN[n2];
505 dy2 = yy-fYN[n2];
506 dx3 = xx-fXN[m];
507 dy3 = yy-fYN[m];
508
509 dd1 = (dx2*dy1-dx1*dy2);
510 dd2 = (dx1*dy2-dx2*dy1);
511
512 if (dd1*dd2!=0) {
513 u = (dx2*dy3-dx3*dy2)/dd1;
514 v = (dx1*dy3-dx3*dy1)/dd2;
515 if ((u<0) || (v<0)) {
516 // No, it cannot - point m does not lie in-between n1 and n2 as
517 // viewed from e. Replace either n1 or n2 to increase the
518 // n1-e-n2 angle. The one to replace is the one which makes the
519 // smallest angle with e->m
520 vNv1 = (dx1*dx3+dy1*dy3)/TMath::Sqrt(dx1*dx1+dy1*dy1);
521 vNv2 = (dx2*dx3+dy2*dy3)/TMath::Sqrt(dx2*dx2+dy2*dy2);
522 if (vNv1 > vNv2) {
523 n1 = m;
524 phi1 = TMath::ATan2(dy3,dx3);
525 phi2 = TMath::ATan2(dy2,dx2);
526 } else {
527 n2 = m;
528 phi1 = TMath::ATan2(dy1,dx1);
529 phi2 = TMath::ATan2(dy3,dx3);
530 }
531 dphi = (phi1-phi2)-((Int_t)((phi1-phi2)/TMath::TwoPi())*TMath::TwoPi());
532 if (dphi < 0) dphi = dphi+TMath::TwoPi();
533 if (((dphi-TMath::Pi())*(lastdphi-TMath::Pi())) < 0) {
534 // The addition of point m means the angle n1-e-n2 has risen
535 // above 180 degs, the point is in the hull.
536 goto L10;
537 }
538 lastdphi = dphi;
539 }
540 }
541 }
542 }
543 // Point e is not surrounded by points - it is not in the hull.
544 goto L999;
545L10:
546 deTinhull = kTRUE;
547L999:
548 return deTinhull;
549}
550
551
552////////////////////////////////////////////////////////////////////////////////
553/// Finds the z-value at point e given that it lies
554/// on the plane defined by t1,t2,t3
555
557{
558 Int_t tmp;
559 Bool_t swap;
560 Double_t x1,x2,x3,y1,y2,y3,f1,f2,f3,u,v,w;
561
562 Int_t t1 = TI1;
563 Int_t t2 = TI2;
564 Int_t t3 = TI3;
565
566 // order the vertices
567L1:
568 swap = kFALSE;
569 if (t2 > t1) { tmp = t1; t1 = t2; t2 = tmp; swap = kTRUE;}
570 if (t3 > t2) { tmp = t2; t2 = t3; t3 = tmp; swap = kTRUE;}
571 if (swap) goto L1;
572
573 x1 = fXN[t1];
574 x2 = fXN[t2];
575 x3 = fXN[t3];
576 y1 = fYN[t1];
577 y2 = fYN[t2];
578 y3 = fYN[t3];
579 f1 = fZ[t1-1];
580 f2 = fZ[t2-1];
581 f3 = fZ[t3-1];
582 u = (f1*(y2-y3)+f2*(y3-y1)+f3*(y1-y2))/(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
583 v = (f1*(x2-x3)+f2*(x3-x1)+f3*(x1-x2))/(y1*(x2-x3)+y2*(x3-x1)+y3*(x1-x2));
584 w = f1-u*x1-v*y1;
585
586 return u*fXN[e]+v*fYN[e]+w;
587}
588
589
590////////////////////////////////////////////////////////////////////////////////
591/// Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
592/// calculate a z-value for it by linearly interpolating the z-values that
593/// make up that triangle.
594
596{
597 Double_t thevalue;
598
599 Int_t it, ntris_tried, p, n, m;
600 Int_t i,j,k,l,z,f,d,o1,o2,a,b,t1,t2,t3;
601 Int_t ndegen=0,degen=0,fdegen=0,o1degen=0,o2degen=0;
602 Double_t vxN,vyN;
603 Double_t d1,d2,d3,c1,c2,dko1,dko2,dfo1;
604 Double_t dfo2,sin_sum,cfo1k,co2o1k,co2o1f;
605
606 Bool_t shouldbein;
607 Double_t dx1,dx2,dx3,dy1,dy2,dy3,u,v,dxz[3],dyz[3];
608
609 // initialise the Delaunay algorithm if needed
610 if (!fInit) {
612 FindHull();
613 fInit = kTRUE;
614 }
615
616 // create vectors needed for sorting
617 if (!fOrder) {
618 fOrder = new Int_t[fNpoints];
619 fDist = new Double_t[fNpoints];
620 }
621
622 // the input point will be point zero.
623 fXN[0] = xx;
624 fYN[0] = yy;
625
626 // set the output value to the default value for now
627 thevalue = fZout;
628
629 // some counting
630 ntris_tried = 0;
631
632 // no point in proceeding if xx or yy are silly
633 if ((xx>fXNmax) || (xx<fXNmin) || (yy>fYNmax) || (yy<fYNmin)) return thevalue;
634
635 // check existing Delaunay triangles for a good one
636 for (it=1; it<=fNdt; it++) {
637 p = fPTried[it-1];
638 n = fNTried[it-1];
639 m = fMTried[it-1];
640 // p, n and m form a previously found Delaunay triangle, does it
641 // enclose the point?
642 if (Enclose(p,n,m,0)) {
643 // yes, we have the triangle
644 thevalue = InterpolateOnPlane(p,n,m,0);
645 return thevalue;
646 }
647 }
648
649 // is this point inside the convex hull?
650 shouldbein = InHull(0,-1);
651 if (!shouldbein) return thevalue;
652
653 // it must be in a Delaunay triangle - find it...
654
655 // order mass points by distance in mass plane from desired point
656 for (it=1; it<=fNpoints; it++) {
657 vxN = fXN[it];
658 vyN = fYN[it];
659 fDist[it-1] = TMath::Sqrt((xx-vxN)*(xx-vxN)+(yy-vyN)*(yy-vyN));
660 }
661
662 // sort array 'fDist' to find closest points
664 for (it=0; it<fNpoints; it++) fOrder[it]++;
665
666 // loop over triplets of close points to try to find a triangle that
667 // encloses the point.
668 for (k=3; k<=fNpoints; k++) {
669 m = fOrder[k-1];
670 for (j=2; j<=k-1; j++) {
671 n = fOrder[j-1];
672 for (i=1; i<=j-1; i++) {
673 p = fOrder[i-1];
674 if (ntris_tried > fMaxIter) {
675 // perhaps this point isn't in the hull after all
676/// Warning("Interpolate",
677/// "Abandoning the effort to find a Delaunay triangle (and thus interpolated z-value) for point %g %g"
678/// ,xx,yy);
679 return thevalue;
680 }
681 ntris_tried++;
682 // check the points aren't colinear
683 d1 = TMath::Sqrt((fXN[p]-fXN[n])*(fXN[p]-fXN[n])+(fYN[p]-fYN[n])*(fYN[p]-fYN[n]));
684 d2 = TMath::Sqrt((fXN[p]-fXN[m])*(fXN[p]-fXN[m])+(fYN[p]-fYN[m])*(fYN[p]-fYN[m]));
685 d3 = TMath::Sqrt((fXN[n]-fXN[m])*(fXN[n]-fXN[m])+(fYN[n]-fYN[m])*(fYN[n]-fYN[m]));
686 if ((d1+d2<=d3) || (d1+d3<=d2) || (d2+d3<=d1)) goto L90;
687
688 // does the triangle enclose the point?
689 if (!Enclose(p,n,m,0)) goto L90;
690
691 // is it a Delaunay triangle? (ie. are there any other points
692 // inside the circle that is defined by its vertices?)
693
694 // test the triangle for Delaunay'ness
695
696 // loop over all other points testing each to see if it's
697 // inside the triangle's circle
698 ndegen = 0;
699 for ( z=1; z<=fNpoints; z++) {
700 if ((z==p) || (z==n) || (z==m)) goto L50;
701 // An easy first check is to see if point z is inside the triangle
702 // (if it's in the triangle it's also in the circle)
703
704 // point z cannot be inside the triangle if it's further from (xx,yy)
705 // than the furthest pointing making up the triangle - test this
706 for (l=1; l<=fNpoints; l++) {
707 if (fOrder[l-1] == z) {
708 if ((l<i) || (l<j) || (l<k)) {
709 // point z is nearer to (xx,yy) than m, n or p - it could be in the
710 // triangle so call enclose to find out
711
712 // if it is inside the triangle this can't be a Delaunay triangle
713 if (Enclose(p,n,m,z)) goto L90;
714 } else {
715 // there's no way it could be in the triangle so there's no point
716 // calling enclose
717 goto L1;
718 }
719 }
720 }
721 // is point z colinear with any pair of the triangle points?
722L1:
723 if (((fXN[p]-fXN[z])*(fYN[p]-fYN[n])) == ((fYN[p]-fYN[z])*(fXN[p]-fXN[n]))) {
724 // z is colinear with p and n
725 a = p;
726 b = n;
727 } else if (((fXN[p]-fXN[z])*(fYN[p]-fYN[m])) == ((fYN[p]-fYN[z])*(fXN[p]-fXN[m]))) {
728 // z is colinear with p and m
729 a = p;
730 b = m;
731 } else if (((fXN[n]-fXN[z])*(fYN[n]-fYN[m])) == ((fYN[n]-fYN[z])*(fXN[n]-fXN[m]))) {
732 // z is colinear with n and m
733 a = n;
734 b = m;
735 } else {
736 a = 0;
737 b = 0;
738 }
739 if (a != 0) {
740 // point z is colinear with 2 of the triangle points, if it lies
741 // between them it's in the circle otherwise it's outside
742 if (fXN[a] != fXN[b]) {
743 if (((fXN[z]-fXN[a])*(fXN[z]-fXN[b])) < 0) {
744 goto L90;
745 } else if (((fXN[z]-fXN[a])*(fXN[z]-fXN[b])) == 0) {
746 // At least two points are sitting on top of each other, we will
747 // treat these as one and not consider this a 'multiple points lying
748 // on a common circle' situation. It is a sign something could be
749 // wrong though, especially if the two coincident points have
750 // different fZ's. If they don't then this is harmless.
751 Warning("Interpolate", "Two of these three points are coincident %d %d %d",a,b,z);
752 }
753 } else {
754 if (((fYN[z]-fYN[a])*(fYN[z]-fYN[b])) < 0) {
755 goto L90;
756 } else if (((fYN[z]-fYN[a])*(fYN[z]-fYN[b])) == 0) {
757 // At least two points are sitting on top of each other - see above.
758 Warning("Interpolate", "Two of these three points are coincident %d %d %d",a,b,z);
759 }
760 }
761 // point is outside the circle, move to next point
762 goto L50;
763 }
764
765 // if point z were to look at the triangle, which point would it see
766 // lying between the other two? (we're going to form a quadrilateral
767 // from the points, and then demand certain properties of that
768 // quadrilateral)
769 dxz[0] = fXN[p]-fXN[z];
770 dyz[0] = fYN[p]-fYN[z];
771 dxz[1] = fXN[n]-fXN[z];
772 dyz[1] = fYN[n]-fYN[z];
773 dxz[2] = fXN[m]-fXN[z];
774 dyz[2] = fYN[m]-fYN[z];
775 for(l=1; l<=3; l++) {
776 dx1 = dxz[l-1];
777 dx2 = dxz[l%3];
778 dx3 = dxz[(l+1)%3];
779 dy1 = dyz[l-1];
780 dy2 = dyz[l%3];
781 dy3 = dyz[(l+1)%3];
782
783 // u et v are used only to know their sign. The previous
784 // code computed them with a division which was long and
785 // might be a division by 0. It is now replaced by a
786 // multiplication.
787 u = (dy3*dx2-dx3*dy2)*(dy1*dx2-dx1*dy2);
788 v = (dy3*dx1-dx3*dy1)*(dy2*dx1-dx2*dy1);
789
790 if ((u>=0) && (v>=0)) {
791 // vector (dx3,dy3) is expressible as a sum of the other two vectors
792 // with positive coefficients -> i.e. it lies between the other two vectors
793 if (l == 1) {
794 f = m;
795 o1 = p;
796 o2 = n;
797 } else if (l == 2) {
798 f = p;
799 o1 = n;
800 o2 = m;
801 } else {
802 f = n;
803 o1 = m;
804 o2 = p;
805 }
806 goto L2;
807 }
808 }
809/// Error("Interpolate", "Should not get to here");
810 // may as well soldier on
811 f = m;
812 o1 = p;
813 o2 = n;
814L2:
815 // this is not a valid quadrilateral if the diagonals don't cross,
816 // check that points f and z lie on opposite side of the line o1-o2,
817 // this is true if the angle f-o1-z is greater than o2-o1-z and o2-o1-f
818 cfo1k = ((fXN[f]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[f]-fYN[o1])*(fYN[z]-fYN[o1]))/
819 TMath::Sqrt(((fXN[f]-fXN[o1])*(fXN[f]-fXN[o1])+(fYN[f]-fYN[o1])*(fYN[f]-fYN[o1]))*
820 ((fXN[z]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[z]-fYN[o1])*(fYN[z]-fYN[o1])));
821 co2o1k = ((fXN[o2]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[z]-fYN[o1]))/
822 TMath::Sqrt(((fXN[o2]-fXN[o1])*(fXN[o2]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[o2]-fYN[o1]))*
823 ((fXN[z]-fXN[o1])*(fXN[z]-fXN[o1]) + (fYN[z]-fYN[o1])*(fYN[z]-fYN[o1])));
824 co2o1f = ((fXN[o2]-fXN[o1])*(fXN[f]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[f]-fYN[o1]))/
825 TMath::Sqrt(((fXN[o2]-fXN[o1])*(fXN[o2]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[o2]-fYN[o1]))*
826 ((fXN[f]-fXN[o1])*(fXN[f]-fXN[o1]) + (fYN[f]-fYN[o1])*(fYN[f]-fYN[o1]) ));
827 if ((cfo1k>co2o1k) || (cfo1k>co2o1f)) {
828 // not a valid quadrilateral - point z is definitely outside the circle
829 goto L50;
830 }
831 // calculate the 2 internal angles of the quadrangle formed by joining
832 // points z and f to points o1 and o2, at z and f. If they sum to less
833 // than 180 degrees then z lies outside the circle
834 dko1 = TMath::Sqrt((fXN[z]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[z]-fYN[o1])*(fYN[z]-fYN[o1]));
835 dko2 = TMath::Sqrt((fXN[z]-fXN[o2])*(fXN[z]-fXN[o2])+(fYN[z]-fYN[o2])*(fYN[z]-fYN[o2]));
836 dfo1 = TMath::Sqrt((fXN[f]-fXN[o1])*(fXN[f]-fXN[o1])+(fYN[f]-fYN[o1])*(fYN[f]-fYN[o1]));
837 dfo2 = TMath::Sqrt((fXN[f]-fXN[o2])*(fXN[f]-fXN[o2])+(fYN[f]-fYN[o2])*(fYN[f]-fYN[o2]));
838 c1 = ((fXN[z]-fXN[o1])*(fXN[z]-fXN[o2])+(fYN[z]-fYN[o1])*(fYN[z]-fYN[o2]))/dko1/dko2;
839 c2 = ((fXN[f]-fXN[o1])*(fXN[f]-fXN[o2])+(fYN[f]-fYN[o1])*(fYN[f]-fYN[o2]))/dfo1/dfo2;
840 sin_sum = c1*TMath::Sqrt(1-c2*c2)+c2*TMath::Sqrt(1-c1*c1);
841
842 // sin_sum doesn't always come out as zero when it should do.
843 if (sin_sum < -1.E-6) {
844 // z is inside the circle, this is not a Delaunay triangle
845 goto L90;
846 } else if (TMath::Abs(sin_sum) <= 1.E-6) {
847 // point z lies on the circumference of the circle (within rounding errors)
848 // defined by the triangle, so there is potential for degeneracy in the
849 // triangle set (Delaunay triangulation does not give a unique way to split
850 // a polygon whose points lie on a circle into constituent triangles). Make
851 // a note of the additional point number.
852 ndegen++;
853 degen = z;
854 fdegen = f;
855 o1degen = o1;
856 o2degen = o2;
857 }
858L50:
859 continue;
860 }
861 // This is a good triangle
862 if (ndegen > 0) {
863 // but is degenerate with at least one other,
864 // haven't figured out what to do if more than 4 points are involved
865/// if (ndegen > 1) {
866/// Error("Interpolate",
867/// "More than 4 points lying on a circle. No decision making process formulated for triangulating this region in a non-arbitrary way %d %d %d %d",
868/// p,n,m,degen);
869/// return thevalue;
870/// }
871
872 // we have a quadrilateral which can be split down either diagonal
873 // (d<->f or o1<->o2) to form valid Delaunay triangles. Choose diagonal
874 // with highest average z-value. Whichever we choose we will have
875 // verified two triangles as good and two as bad, only note the good ones
876 d = degen;
877 f = fdegen;
878 o1 = o1degen;
879 o2 = o2degen;
880 if ((fZ[o1-1]+fZ[o2-1]) > (fZ[d-1]+fZ[f-1])) {
881 // best diagonalisation of quadrilateral is current one, we have
882 // the triangle
883 t1 = p;
884 t2 = n;
885 t3 = m;
886 // file the good triangles
887 FileIt(p, n, m);
888 FileIt(d, o1, o2);
889 } else {
890 // use other diagonal to split quadrilateral, use triangle formed by
891 // point f, the degnerate point d and whichever of o1 and o2 create
892 // an enclosing triangle
893 t1 = f;
894 t2 = d;
895 if (Enclose(f,d,o1,0)) {
896 t3 = o1;
897 } else {
898 t3 = o2;
899 }
900 // file the good triangles
901 FileIt(f, d, o1);
902 FileIt(f, d, o2);
903 }
904 } else {
905 // this is a Delaunay triangle, file it
906 FileIt(p, n, m);
907 t1 = p;
908 t2 = n;
909 t3 = m;
910 }
911 // do the interpolation
912 thevalue = InterpolateOnPlane(t1,t2,t3,0);
913 return thevalue;
914L90:
915 continue;
916 }
917 }
918 }
919 if (shouldbein) {
920 Error("Interpolate",
921 "Point outside hull when expected inside: this point could be dodgy %g %g %d",
922 xx, yy, ntris_tried);
923 }
924 return thevalue;
925}
926
927
928////////////////////////////////////////////////////////////////////////////////
929/// Defines the number of triangles tested for a Delaunay triangle
930/// (number of iterations) before abandoning the search
931
933{
934 fAllTri = kFALSE;
935 fMaxIter = n;
936}
937
938
939////////////////////////////////////////////////////////////////////////////////
940/// Sets the histogram bin height for points lying outside the convex hull ie:
941/// the bins in the margin.
942
944{
945 fZout = z;
946}
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
static const double x2[5]
static const double x1[5]
static const double x3[11]
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:91
#define ClassImp(name)
Definition Rtypes.h:364
float xmin
float ymin
float xmax
float ymax
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition TGraph2D.h:41
Double_t * GetY() const
Definition TGraph2D.h:122
Double_t GetYmin() const
Returns the Y minimum.
Double_t * GetX() const
Definition TGraph2D.h:121
Double_t GetXmin() const
Returns the X minimum.
Double_t GetXmax() const
Returns the X maximum.
Int_t GetN() const
Definition TGraph2D.h:120
Double_t GetYmax() const
Returns the Y maximum.
Double_t * GetZ() const
Definition TGraph2D.h:123
TGraphDelaunay generates a Delaunay triangulation of a TGraph2D.
Int_t fNpoints
!Number of data points in fGraph2D
Double_t ComputeZ(Double_t x, Double_t y)
Return the z value corresponding to the (x,y) point in fGraph2D.
Int_t * fNTried
!Delaunay triangles storage of size fNdt
Double_t * fY
!Pointer to fGraph2D->fY
Double_t fXNmax
!Maximum value of fXN
Bool_t InHull(Int_t E, Int_t X) const
Is point e inside the hull defined by all points apart from x ?
Double_t fXNmin
!Minimum value of fXN
Double_t InterpolateOnPlane(Int_t TI1, Int_t TI2, Int_t TI3, Int_t E) const
Finds the z-value at point e given that it lies on the plane defined by t1,t2,t3.
Double_t * fYN
!fGraph2D vectors normalized of size fNpoints
Double_t fYoffset
!Parameters used to normalize user data
void FindHull()
Finds those points which make up the convex hull of the set.
void SetMarginBinsContent(Double_t z=0.)
Sets the histogram bin height for points lying outside the convex hull ie: the bins in the margin.
Double_t fXoffset
!
Double_t * fDist
!Array used to order mass points by distance
void FileIt(Int_t P, Int_t N, Int_t M)
Files the triangle defined by the 3 vertices p, n and m into the fxTried arrays.
Double_t fYNmin
!Minimum value of fYN
Int_t fTriedSize
!Real size of the fxTried arrays
Double_t * fX
!Pointer to fGraph2D->fX
TGraph2D * fGraph2D
!2D graph containing the user data
void FindAllTriangles()
Attempt to find all the Delaunay triangles of the point set.
Bool_t fAllTri
!True if FindAllTriangles() has been performed on fGraph2D
void SetMaxIter(Int_t n=100000)
Defines the number of triangles tested for a Delaunay triangle (number of iterations) before abandoni...
Int_t fMaxIter
!Maximum number of iterations to find Delaunay triangles
Bool_t fInit
!True if CreateTrianglesDataStructure() and FindHull() have been performed
Int_t * fOrder
!Array used to order mass points by distance
Bool_t Enclose(Int_t T1, Int_t T2, Int_t T3, Int_t Ex) const
Is point e inside the triangle t1-t2-t3 ?
Double_t fXScaleFactor
!
virtual ~TGraphDelaunay()
TGraphDelaunay destructor.
Int_t fNdt
!Number of Delaunay triangles found
Double_t * fZ
!Pointer to fGraph2D->fZ
Double_t fYNmax
!Maximum value of fYN
Double_t fZout
!Histogram bin height for points lying outside the convex hull
Double_t fYScaleFactor
!
void CreateTrianglesDataStructure()
Function used internally only.
Int_t fNhull
!Number of points in the hull
Double_t Interpolate(Double_t x, Double_t y)
Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and calculate a z-value for it by...
TGraphDelaunay()
TGraphDelaunay default constructor.
Int_t * fHullPoints
!Hull points of size fNhull
Double_t * fXN
!fGraph2D vectors normalized of size fNpoints
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:879
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:893
#define Interpolate(a, x, b, y)
Definition geom.c:179
return c1
Definition legend1.C:41
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TF1 * f1
Definition legend1.C:11
return c2
Definition legend2.C:14
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:212
Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y)
Function which returns kTRUE if point xp,yp lies inside the polygon defined by the np points in array...
Definition TMath.h:1209
Double_t ATan2(Double_t y, Double_t x)
Definition TMath.h:679
Double_t Sqrt(Double_t x)
Definition TMath.h:691
constexpr Double_t Pi()
Definition TMath.h:37
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition TMathBase.h:362
Short_t Abs(Short_t d)
Definition TMathBase.h:120
constexpr Double_t TwoPi()
Definition TMath.h:44
auto * m
Definition textangle.C:8
auto * l
Definition textangle.C:4
auto * t1
Definition textangle.C:20