ROOT   6.08/07 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
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
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->fY.
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
Number of points in the hull.
return c1
Definition: legend1.C:41
void swap(TDirectoryEntry &e1, TDirectoryEntry &e2) noexcept
float ymin
Definition: THbookFile.cxx:93
Double_t fYNmin
Maximum value of fXN.
Double_t GetXmax() const
Returns the X maximum.
Definition: TGraph2D.cxx:1197
Int_t fNhull
Number of data points in fGraph2D.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
Double_t fYNmax
Minimum value of fYN.
Int_t fTriedSize
Maximum number of iterations to find Delaunay triangles.
Double_t * fDist
Histogram bin height for points lying outside the convex hull.
Double_t GetYmax() const
Returns the Y maximum.
Definition: TGraph2D.cxx:1219
Double_t * fXN
Pointer to fGraph2D->fZ.
Double_t fXoffset
Maximum value of fYN.
TLatex * t1
Definition: textangle.C:20
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
Double_t fYoffset
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:33
Double_t fXNmax
Minimum 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:989
Double_t ATan2(Double_t, Double_t)
Definition: TMath.h:454
Bool_t IsInside(T xp, T yp, Int_t np, T *x, T *y)
Definition: TMath.h:1043
Double_t TwoPi()
Definition: TMath.h:45
Int_t * fHullPoints
float ymax
Definition: THbookFile.cxx:93
Double_t fYScaleFactor
SVector< double, 2 > v
Definition: Dict.h:5
Int_t fNpoints
Number of Delaunay triangles found.
TMarker * m
Definition: textangle.C:8
Double_t fXScaleFactor
Parameters used to normalize user data.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
Int_t * fOrder
Hull points of size fNhull.
Double_t E()
Definition: TMath.h:54
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
fGraph2D vectors normalized of size fNpoints
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
void CreateTrianglesDataStructure()
2D graph containing the user data
Double_t * GetY() const
Definition: TGraph2D.h:129
Double_t Pi()
Definition: TMath.h:44
Double_t * fY
Pointer to fGraph2D->fX.
Bool_t fInit
True if FindAllTriangles() has been performed on fGraph2D.
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:279
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:130
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
True if CreateTrianglesDataStructure() and FindHull() have been performed.
Bool_t fAllTri
Array used to order mass points by distance.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Int_t * fPTried
Real size of the fxTried arrays.
double f2(const double *x)
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
Double_t GetYmin() const
Returns the Y minimum.
Definition: TGraph2D.cxx:1230
Int_t GetN() const
Definition: TGraph2D.h:127
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
Delaunay triangles storage of size fNdt.
const int nn
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
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:50
const Bool_t kTRUE
Definition: Rtypes.h:91
virtual ~TGraphDelaunay()
TGraphDelaunay destructor.
Int_t fMaxIter
Array used to order mass points by distance.
const Int_t n
Definition: legend1.C:16
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:911
TGraphDelaunay generates a Delaunay triangulation of a TGraph2D.
Double_t * GetX() const
Definition: TGraph2D.h:128
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]