Logo ROOT   6.10/09
Reference Guide
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
21 TGraphDelaunay generates a Delaunay triangulation of a TGraph2D. This
22 triangulation code derives from an implementation done by Luke Jones
23 (Royal Holloway, University of London) in April 2002 in the PAW context.
24 
25 This software cannot be guaranteed to work under all circumstances. They
26 were originally written to work with a few hundred points in an XY space
27 with similar X and Y ranges.
28 
29 Definition of Delaunay triangulation (After B. Delaunay):
30 For a set S of points in the Euclidean plane, the unique triangulation DT(S)
31 of S such that no point in S is inside the circumcircle of any triangle in
32 DT(S). DT(S) is the dual of the Voronoi diagram of S. If n is the number of
33 points in S, the Voronoi diagram of S is the partitioning of the plane
34 containing S points into n convex polygons such that each polygon contains
35 exactly one point and every point in a given polygon is closer to its
36 central point than to any other. A Voronoi diagram is sometimes also known
37 as 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
42 practical 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();
93  fNpoints = fGraph2D->GetN();
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);
191  fXNmax = (xmax+fXoffset)*fXScaleFactor;
192  fXNmin = (xmin+fXoffset)*fXScaleFactor;
193  fYNmax = (ymax+fYoffset)*fYScaleFactor;
194  fYNmin = (ymin+fYoffset)*fYScaleFactor;
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.
205  fTriedSize = 2*fNpoints;
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
243 L1:
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 
426  if (!fHullPoints) fHullPoints = new Int_t[fNpoints];
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;
545 L10:
546  deTinhull = kTRUE;
547 L999:
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
567 L1:
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?
722 L1:
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;
814 L2:
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  }
858 L50:
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;
914 L90:
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 }
const int nx
Definition: kalman.C:16
Double_t * fZ
!Pointer to fGraph2D->fZ
float xmin
Definition: THbookFile.cxx:93
void SetMaxIter(Int_t n=100000)
Defines the number of triangles tested for a Delaunay triangle (number of iterations) before abandoni...
static double p3(double t, double a, double b, double c, double d)
Double_t * fYN
!fGraph2D vectors normalized of size fNpoints
Double_t * fX
!Pointer to fGraph2D->fX
return c1
Definition: legend1.C:41
void swap(TDirectoryEntry &e1, TDirectoryEntry &e2) noexcept
constexpr Double_t TwoPi()
Definition: TMath.h:44
float ymin
Definition: THbookFile.cxx:93
Double_t fYNmin
!Minimum value of fYN
Double_t fZout
!Histogram bin height for points lying outside the convex hull
Double_t GetXmax() const
Returns the X maximum.
Definition: TGraph2D.cxx:1197
Int_t fNhull
!Number of points in the hull
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
Double_t fYNmax
!Maximum value of fYN
Int_t fTriedSize
!Real size of the fxTried arrays
Double_t * fDist
!Array used to order mass points by distance
Double_t GetYmax() const
Returns the Y maximum.
Definition: TGraph2D.cxx:1219
Double_t * fXN
!fGraph2D vectors normalized of size fNpoints
Double_t fXoffset
!
TLatex * t1
Definition: textangle.C:20
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
Double_t fYoffset
!Parameters used to normalize user data
Bool_t InHull(Int_t E, Int_t X) const
Is point e inside the hull defined by all points apart from x ?
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.
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
Double_t ComputeZ(Double_t x, Double_t y)
Return the z value corresponding to the (x,y) point in fGraph2D.
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
Double_t fXNmax
!Maximum value of fXN
const int ny
Definition: kalman.C:17
Double_t GetXmin() const
Returns the X minimum.
Definition: TGraph2D.cxx:1208
void SetMarginBinsContent(Double_t z=0.)
Sets the histogram bin height for points lying outside the convex hull ie: the bins in the margin...
static double p2(double t, double a, double b, double c)
TGraphDelaunay()
TGraphDelaunay default constructor.
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMath.h:1151
Double_t ATan2(Double_t, Double_t)
Definition: TMath.h:581
Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y)
Definition: TMath.h:1205
constexpr Double_t Pi()
Definition: TMath.h:40
Int_t * fHullPoints
!Hull points of size fNhull
Int_t * fNTried
!Delaunay triangles storage of size fNdt
float ymax
Definition: THbookFile.cxx:93
Double_t fYScaleFactor
!
SVector< double, 2 > v
Definition: Dict.h:5
Int_t fNpoints
!Number of data points in fGraph2D
Int_t fNdt
!Number of Delaunay triangles found
TMarker * m
Definition: textangle.C:8
Double_t fXScaleFactor
!
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:873
Int_t * fOrder
!Array used to order mass points by distance
TLine * l
Definition: textangle.C:4
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 fXNmin
!Minimum value of fXN
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
void CreateTrianglesDataStructure()
Function used internally only.
constexpr Double_t E()
Definition: TMath.h:74
const Bool_t kFALSE
Definition: RtypesCore.h:92
Double_t * GetY() const
Definition: TGraph2D.h:119
Double_t * fY
!Pointer to fGraph2D->fY
Bool_t fInit
!True if CreateTrianglesDataStructure() and FindHull() have been performed
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 ?
return c2
Definition: legend2.C:14
static const double x1[5]
#define ClassImp(name)
Definition: Rtypes.h:336
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
Double_t * GetZ() const
Definition: TGraph2D.h:120
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
TGraph2D * fGraph2D
!2D graph containing the user data
Bool_t fAllTri
!True if FindAllTriangles() has been performed on fGraph2D
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Int_t * fPTried
!
double f2(const double *x)
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
Double_t GetYmin() const
Returns the Y minimum.
Definition: TGraph2D.cxx:1230
Int_t GetN() const
Definition: TGraph2D.h:117
TF1 * f1
Definition: legend1.C:11
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Int_t * fMTried
!
const int nn
Double_t Sqrt(Double_t x)
Definition: TMath.h:591
void FindAllTriangles()
Attempt to find all the Delaunay triangles of the point set.
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:40
virtual ~TGraphDelaunay()
TGraphDelaunay destructor.
const Bool_t kTRUE
Definition: RtypesCore.h:91
Int_t fMaxIter
!Maximum number of iterations to find Delaunay triangles
const Int_t n
Definition: legend1.C:16
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:859
TGraphDelaunay generates a Delaunay triangulation of a TGraph2D.
Double_t * GetX() const
Definition: TGraph2D.h:118
void FindHull()
Finds those points which make up the convex hull of the set.
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...
static const double x3[11]