// @(#)root/hist:$Name: $:$Id: TGraphDelaunay.cxx,v 1.00
// Author: Olivier Couet, Luke Jones (Royal Holloway, University of London)
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
#include "TROOT.h"
#include "TMath.h"
#include "TGraphDelaunay.h"
ClassImp(TGraphDelaunay)
//______________________________________________________________________________
//
// TGraphDelaunay generates a Delaunay triangulation of a TGraph2D. This
// triangulation code derives from an implementation done by Luke Jones
// (Royal Holloway, University of London) in April 2002 in the PAW context.
//
// This software cannot be guaranteed to work under all circumstances. They
// were originally written to work with a few hundred points in an XY space
// with similar X and Y ranges.
//
// Definition of Delaunay triangulation (After B. Delaunay):
// For a set S of points in the Euclidean plane, the unique triangulation DT(S)
// of S such that no point in S is inside the circumcircle of any triangle in
// DT(S). DT(S) is the dual of the Voronoi diagram of S. If n is the number of
// points in S, the Voronoi diagram of S is the partitioning of the plane
// containing S points into n convex polygons such that each polygon contains
// exactly one point and every point in a given polygon is closer to its
// central point than to any other. A Voronoi diagram is sometimes also known
// as a Dirichlet tessellation.
//
/*
This applet
gives a nice practical view of Delaunay triangulation and Voronoi diagram.
*/
//
//______________________________________________________________________________
TGraphDelaunay::TGraphDelaunay()
: TNamed("TGraphDelaunay","TGraphDelaunay")
{
// TGraphDelaunay default constructor
fGraph2D = 0;
fX = 0;
fY = 0;
fZ = 0;
fNpoints = 0;
fTriedSize = 0;
fZout = 0.;
fNdt = 0;
fNhull = 0;
fHullPoints = 0;
fXN = 0;
fYN = 0;
fOrder = 0;
fDist = 0;
fPTried = 0;
fNTried = 0;
fMTried = 0;
fInit = kFALSE;
SetMaxIter();
}
//______________________________________________________________________________
TGraphDelaunay::TGraphDelaunay(TGraph2D *g)
: TNamed("TGraphDelaunay","TGraphDelaunay")
{
// TGraphDelaunay default constructor
fGraph2D = g;
fX = fGraph2D->GetX();
fY = fGraph2D->GetY();
fZ = fGraph2D->GetZ();
fNpoints = fGraph2D->GetN();
fTriedSize = 0;
fZout = 0.;
fNdt = 0;
fNhull = 0;
fHullPoints = 0;
fXN = 0;
fYN = 0;
fOrder = 0;
fDist = 0;
fPTried = 0;
fNTried = 0;
fMTried = 0;
fInit = kFALSE;
SetMaxIter();
}
//______________________________________________________________________________
TGraphDelaunay::~TGraphDelaunay()
{
// TGraphDelaunay destructor.
if (fPTried) delete [] fPTried;
if (fNTried) delete [] fNTried;
if (fMTried) delete [] fMTried;
if (fHullPoints) delete [] fHullPoints;
if (fOrder) delete [] fOrder;
if (fDist) delete [] fDist;
if (fXN) delete [] fXN;
if (fYN) delete [] fYN;
fPTried = 0;
fNTried = 0;
fMTried = 0;
fHullPoints = 0;
fOrder = 0;
fDist = 0;
fXN = 0;
fYN = 0;
}
//______________________________________________________________________________
Double_t TGraphDelaunay::ComputeZ(Double_t x, Double_t y)
{
// Return the z value corresponding to the (x,y) point in fGraph2D
Double_t xx, yy;
xx = (x+fXoffset)*fScaleFactor;
yy = (y+fYoffset)*fScaleFactor;
return Interpolate(xx, yy);
}
//______________________________________________________________________________
void TGraphDelaunay::CreateTrianglesDataStructure()
{
// Fonction used internally only. It creates the data structures needed to
// compute the Delaunay triangles.
// Offset fX and fY so they average zero, and scale so the average
// of the X and Y ranges is one. The normalized version of fX and fY used
// in Interpolate.
Double_t xmax = fGraph2D->GetXmax();
Double_t ymax = fGraph2D->GetYmax();
Double_t xmin = fGraph2D->GetXmin();
Double_t ymin = fGraph2D->GetYmin();
fXoffset = -(xmax+xmin)/2.;
fYoffset = -(ymax+ymin)/2.;
fScaleFactor = 2./((xmax-xmin)+(ymax-ymin));
fXNmax = (xmax+fXoffset)*fScaleFactor;
fXNmin = (xmin+fXoffset)*fScaleFactor;
fYNmax = (ymax+fYoffset)*fScaleFactor;
fYNmin = (ymin+fYoffset)*fScaleFactor;
fXN = new Double_t[fNpoints+1];
fYN = new Double_t[fNpoints+1];
for (Int_t n=0; n<fNpoints; n++) {
fXN[n+1] = (fX[n]+fXoffset)*fScaleFactor;
fYN[n+1] = (fY[n]+fYoffset)*fScaleFactor;
}
// If needed, creates the arrays to hold the Delaunay triangles.
// A maximum number of 2*fNpoints is guessed. If more triangles will be
// find, FillIt will automatically enlarge these arrays.
fTriedSize = 2*fNpoints;
fPTried = new Int_t[fTriedSize];
fNTried = new Int_t[fTriedSize];
fMTried = new Int_t[fTriedSize];
}
//______________________________________________________________________________
Bool_t TGraphDelaunay::Enclose(Int_t t1, Int_t t2, Int_t t3, Int_t e) const
{
// Is point e inside the triangle t1-t2-t3 ?
Int_t a = 0, b = 0;
Double_t dx1,dx2,dx3,dy1,dy2,dy3,u,v;
// First ask if point e is colinear with any pair of the triangle points
if (((fXN[t1]-fXN[e])*(fYN[t1]-fYN[t2])) == ((fYN[t1]-fYN[e])*(fXN[t1]-fXN[t2]))) {
// e is colinear with t1 and t2
a = t1;
b = t2;
} else if (((fXN[t1]-fXN[e])*(fYN[t1]-fYN[t3])) == ((fYN[t1]-fYN[e])*(fXN[t1]-fXN[t3]))) {
// e is colinear with t1 and t3
a = t1;
b = t3;
} else if (((fXN[t2]-fXN[e])*(fYN[t2]-fYN[t3])) == ((fYN[t2]-fYN[e])*(fXN[t2]-fXN[t3]))) {
// e is colinear with t2 and t3
a = t2;
b = t3;
}
if (a != 0) {
// point e is colinear with 2 of the triangle points, if it lies
// between them it's in the circle otherwise it's outside
if (fXN[a] != fXN[b]) {
if (((fXN[e]-fXN[a])*(fXN[e]-fXN[b])) <= 0) return kTRUE;
} else {
if (((fYN[e]-fYN[a])*(fYN[e]-fYN[b])) <= 0) return kTRUE;
}
// point is outside the triangle
return kFALSE;
}
// e is not colinear with any pair of triangle points, if it is inside
// the triangle then the vector from e to one of the corners must be
// expressible as a sum with positive coefficients of the vectors from
// the two other corners to e. Say vector3=u*vector1+v*vector2
// vector1==t1->e
dx1 = fXN[e]-fXN[t1];
dy1 = fYN[e]-fYN[t1];
// vector2==t2->e
dx2 = fXN[e]-fXN[t2];
dy2 = fYN[e]-fYN[t2];
// vector3==e->t3
dx3 = fXN[t3]-fXN[e];
dy3 = fYN[t3]-fYN[e];
u = (dx2*dy3-dx3*dy2)/(dx2*dy1-dx1*dy2);
v = (dx1*dy3-dx3*dy1)/(dx1*dy2-dx2*dy1);
if ((u>=0) && (v>=0)) return kTRUE;
return kFALSE;
}
//______________________________________________________________________________
void TGraphDelaunay::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. If these arrays are to small they are automatically
// expanded.
Bool_t swap;
Int_t tmp, ps = p, ns = n, ms = m;
// order the vertices before storing them
L1:
swap = kFALSE;
if (ns > ps) { tmp = ps; ps = ns; ns = tmp; swap = kTRUE;}
if (ms > ns) { tmp = ns; ns = ms; ms = tmp; swap = kTRUE;}
if (swap) goto L1;
// expand the triangles storage if needed
if (fNdt> fTriedSize) {
Int_t newN = 2*fTriedSize;
Int_t *savep = new Int_t [newN];
Int_t *saven = new Int_t [newN];
Int_t *savem = new Int_t [newN];
memcpy(savep,fPTried,fTriedSize*sizeof(Int_t));
memset(&savep[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
delete [] fPTried;
memcpy(saven,fNTried,fTriedSize*sizeof(Int_t));
memset(&saven[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
delete [] fNTried;
memcpy(savem,fMTried,fTriedSize*sizeof(Int_t));
memset(&savem[fTriedSize],0,(newN-fTriedSize)*sizeof(Int_t));
delete [] fMTried;
fPTried = savep;
fNTried = saven;
fMTried = savem;
fTriedSize = newN;
}
// store a new Delaunay triangle
fNdt++;
fPTried[fNdt-1] = ps;
fNTried[fNdt-1] = ns;
fMTried[fNdt-1] = ms;
}
//______________________________________________________________________________
void TGraphDelaunay::FindAllTriangles()
{
// Attempt to find all the Delaunay triangles of the point set. It is not
// guaranteed that it will fully succeed, and no check is made that it has
// fully succeeded (such a check would be possible by referencing the points
// that make up the convex hull). The method is to check if each triangle
// shares all three of its sides with other triangles. If not, a point is
// generated just outside the triangle on the side(s) not shared, and a new
// triangle is found for that point. If this method is not working properly
// (many triangles are not being found) it's probably because the new points
// are too far beyond or too close to the non-shared sides. Fiddling with
// the size of the `alittlebit' parameter may help.
if (fAllTri) return; else fAllTri = kTRUE;
Double_t xcntr,ycntr,xm,ym,xx,yy;
Double_t sx,sy,nx,ny,mx,my,mdotn,nn,a;
Int_t t1,t2,pa,na,ma,pb,nb,mb,p1=0,p2=0,m,n,p3=0;
Bool_t s[3];
Double_t alittlebit = 0.0001;
// start with a point that is guaranteed to be inside the hull (the
// centre of the hull). The starting point is shifted "a little bit"
// otherwise, in case of triangles aligned on a regular grid, we may
// found none of them.
xcntr = 0;
ycntr = 0;
for (n=1; n<=fNhull; n++) {
xcntr = xcntr+fXN[fHullPoints[n-1]];
ycntr = ycntr+fYN[fHullPoints[n-1]];
}
xcntr = xcntr/fNhull+alittlebit;
ycntr = ycntr/fNhull+alittlebit;
// and calculate it's triangle
Interpolate(xcntr,ycntr);
// loop over all Delaunay triangles (including those constantly being
// produced within the loop) and check to see if their 3 sides also
// correspond to the sides of other Delaunay triangles, i.e. that they
// have all their neighbours.
t1 = 1;
while (t1 <= fNdt) {
// get the three points that make up this triangle
pa = fPTried[t1-1];
na = fNTried[t1-1];
ma = fMTried[t1-1];
// produce three integers which will represent the three sides
s[0] = kFALSE;
s[1] = kFALSE;
s[2] = kFALSE;
// loop over all other Delaunay triangles
for (t2=1; t2<=fNdt; t2++) {
if (t2 != t1) {
// get the points that make up this triangle
pb = fPTried[t2-1];
nb = fNTried[t2-1];
mb = fMTried[t2-1];
// do triangles t1 and t2 share a side?
if ((pa==pb && na==nb) || (pa==pb && na==mb) || (pa==nb && na==mb)) {
// they share side 1
s[0] = kTRUE;
} else if ((pa==pb && ma==nb) || (pa==pb && ma==mb) || (pa==nb && ma==mb)) {
// they share side 2
s[1] = kTRUE;
} else if ((na==pb && ma==nb) || (na==pb && ma==mb) || (na==nb && ma==mb)) {
// they share side 3
s[2] = kTRUE;
}
}
// if t1 shares all its sides with other Delaunay triangles then
// forget about it
if (s[0] && s[1] && s[2]) continue;
}
// Looks like t1 is missing a neighbour on at least one side.
// For each side, take a point a little bit beyond it and calculate
// the Delaunay triangle for that point, this should be the triangle
// which shares the side.
for (m=1; m<=3; m++) {
if (!s[m-1]) {
// get the two points that make up this side
if (m == 1) {
p1 = pa;
p2 = na;
p3 = ma;
} else if (m == 2) {
p1 = pa;
p2 = ma;
p3 = na;
} else if (m == 3) {
p1 = na;
p2 = ma;
p3 = pa;
}
// get the coordinates of the centre of this side
xm = (fXN[p1]+fXN[p2])/2.;
ym = (fYN[p1]+fYN[p2])/2.;
// we want to add a little to these coordinates to get a point just
// outside the triangle; (sx,sy) will be the vector that represents
// the side
sx = fXN[p1]-fXN[p2];
sy = fYN[p1]-fYN[p2];
// (nx,ny) will be the normal to the side, but don't know if it's
// pointing in or out yet
nx = sy;
ny = -sx;
nn = TMath::Sqrt(nx*nx+ny*ny);
nx = nx/nn;
ny = ny/nn;
mx = fXN[p3]-xm;
my = fYN[p3]-ym;
mdotn = mx*nx+my*ny;
if (mdotn > 0) {
// (nx,ny) is pointing in, we want it pointing out
nx = -nx;
ny = -ny;
}
// increase/decrease xm and ym a little to produce a point
// just outside the triangle (ensuring that the amount added will
// be large enough such that it won't be lost in rounding errors)
a = TMath::Abs(TMath::Max(alittlebit*xm,alittlebit*ym));
xx = xm+nx*a;
yy = ym+ny*a;
// try and find a new Delaunay triangle for this point
Interpolate(xx,yy);
// this side of t1 should now, hopefully, if it's not part of the
// hull, be shared with a new Delaunay triangle just calculated by Interpolate
}
}
t1++;
}
}
//______________________________________________________________________________
void TGraphDelaunay::FindHull()
{
// Finds those points which make up the convex hull of the set. If the xy
// plane were a sheet of wood, and the points were nails hammered into it
// at the respective coordinates, then if an elastic band were stretched
// over all the nails it would form the shape of the convex hull. Those
// nails in contact with it are the points that make up the hull.
Int_t n,nhull_tmp;
Bool_t in;
if (!fHullPoints) fHullPoints = new Int_t[fNpoints];
nhull_tmp = 0;
for(n=1; n<=fNpoints; n++) {
// if the point is not inside the hull of the set of all points
// bar it, then it is part of the hull of the set of all points
// including it
in = InHull(n,n);
if (!in) {
// cannot increment fNhull directly - InHull needs to know that
// the hull has not yet been completely found
nhull_tmp++;
fHullPoints[nhull_tmp-1] = n;
}
}
fNhull = nhull_tmp;
}
//______________________________________________________________________________
Bool_t TGraphDelaunay::InHull(Int_t e, Int_t x) const
{
// Is point e inside the hull defined by all points apart from x ?
Int_t n1,n2,n,m,ntry;
Double_t lastdphi,dd1,dd2,dx1,dx2,dx3,dy1,dy2,dy3;
Double_t u,v,vNv1,vNv2,phi1,phi2,dphi,xx,yy;
Bool_t deTinhull = kFALSE;
xx = fXN[e];
yy = fYN[e];
if (fNhull > 0) {
// The hull has been found - no need to use any points other than
// those that make up the hull
ntry = fNhull;
} else {
// The hull has not yet been found, will have to try every point
ntry = fNpoints;
}
// n1 and n2 will represent the two points most separated by angle
// from point e. Initially the angle between them will be <180 degs.
// But subsequent points will increase the n1-e-n2 angle. If it
// increases above 180 degrees then point e must be surrounded by
// points - it is not part of the hull.
n1 = 1;
n2 = 2;
if (n1 == x) {
n1 = n2;
n2++;
} else if (n2 == x) {
n2++;
}
// Get the angle n1-e-n2 and set it to lastdphi
dx1 = xx-fXN[n1];
dy1 = yy-fYN[n1];
dx2 = xx-fXN[n2];
dy2 = yy-fYN[n2];
phi1 = TMath::ATan2(dy1,dx1);
phi2 = TMath::ATan2(dy2,dx2);
dphi = (phi1-phi2)-((Int_t)((phi1-phi2)/TMath::TwoPi())*TMath::TwoPi());
if (dphi < 0) dphi = dphi+TMath::TwoPi();
lastdphi = dphi;
for (n=1; n<=ntry; n++) {
if (fNhull > 0) {
// Try hull point n
m = fHullPoints[n-1];
} else {
m = n;
}
if ((m!=n1) && (m!=n2) && (m!=x)) {
// Can the vector e->m be represented as a sum with positive
// coefficients of vectors e->n1 and e->n2?
dx1 = xx-fXN[n1];
dy1 = yy-fYN[n1];
dx2 = xx-fXN[n2];
dy2 = yy-fYN[n2];
dx3 = xx-fXN[m];
dy3 = yy-fYN[m];
dd1 = (dx2*dy1-dx1*dy2);
dd2 = (dx1*dy2-dx2*dy1);
if (dd1*dd2!=0) {
u = (dx2*dy3-dx3*dy2)/dd1;
v = (dx1*dy3-dx3*dy1)/dd2;
if ((u<0) || (v<0)) {
// No, it cannot - point m does not lie inbetween n1 and n2 as
// viewed from e. Replace either n1 or n2 to increase the
// n1-e-n2 angle. The one to replace is the one which makes the
// smallest angle with e->m
vNv1 = (dx1*dx3+dy1*dy3)/TMath::Sqrt(dx1*dx1+dy1*dy1);
vNv2 = (dx2*dx3+dy2*dy3)/TMath::Sqrt(dx2*dx2+dy2*dy2);
if (vNv1 > vNv2) {
n1 = m;
phi1 = TMath::ATan2(dy3,dx3);
phi2 = TMath::ATan2(dy2,dx2);
} else {
n2 = m;
phi1 = TMath::ATan2(dy1,dx1);
phi2 = TMath::ATan2(dy3,dx3);
}
dphi = (phi1-phi2)-((Int_t)((phi1-phi2)/TMath::TwoPi())*TMath::TwoPi());
if (dphi < 0) dphi = dphi+TMath::TwoPi();
if (((dphi-TMath::Pi())*(lastdphi-TMath::Pi())) < 0) {
// The addition of point m means the angle n1-e-n2 has risen
// above 180 degs, the point is in the hull.
goto L10;
}
lastdphi = dphi;
}
}
}
}
// Point e is not surrounded by points - it is not in the hull.
goto L999;
L10:
deTinhull = kTRUE;
L999:
return deTinhull;
}
//______________________________________________________________________________
Double_t TGraphDelaunay::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
Int_t tmp;
Bool_t swap;
Double_t x1,x2,x3,y1,y2,y3,f1,f2,f3,u,v,w;
Int_t t1 = TI1;
Int_t t2 = TI2;
Int_t t3 = TI3;
// order the vertices
L1:
swap = kFALSE;
if (t2 > t1) { tmp = t1; t1 = t2; t2 = tmp; swap = kTRUE;}
if (t3 > t2) { tmp = t2; t2 = t3; t3 = tmp; swap = kTRUE;}
if (swap) goto L1;
x1 = fXN[t1];
x2 = fXN[t2];
x3 = fXN[t3];
y1 = fYN[t1];
y2 = fYN[t2];
y3 = fYN[t3];
f1 = fZ[t1-1];
f2 = fZ[t2-1];
f3 = fZ[t3-1];
u = (f1*(y2-y3)+f2*(y3-y1)+f3*(y1-y2))/(x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2));
v = (f1*(x2-x3)+f2*(x3-x1)+f3*(x1-x2))/(y1*(x2-x3)+y2*(x3-x1)+y3*(x1-x2));
w = f1-u*x1-v*y1;
return u*fXN[e]+v*fYN[e]+w;
}
//______________________________________________________________________________
Double_t TGraphDelaunay::Interpolate(Double_t xx, Double_t yy)
{
// Finds the Delaunay triangle that the point (xi,yi) sits in (if any) and
// calculate a z-value for it by linearly interpolating the z-values that
// make up that triangle.
Double_t thevalue;
Int_t it, ntris_tried, p, n, m;
Int_t i,j,k,l,z,f,d,o1,o2,a,b,t1,t2,t3;
Int_t ndegen=0,degen=0,fdegen=0,o1degen=0,o2degen=0;
Double_t vxN,vyN;
Double_t d1,d2,d3,c1,c2,dko1,dko2,dfo1;
Double_t dfo2,sin_sum,cfo1k,co2o1k,co2o1f;
Bool_t shouldbein;
Double_t dx1,dx2,dx3,dy1,dy2,dy3,u,v,dxz[3],dyz[3];
// initialise the Delaunay algorithm
if (!fInit) {
CreateTrianglesDataStructure();
FindHull();
fInit = kTRUE;
}
// create vectors needed for sorting
if (!fOrder) {
fOrder = new Int_t[fNpoints];
fDist = new Double_t[fNpoints];
}
// the input point will be point zero.
fXN[0] = xx;
fYN[0] = yy;
// set the output value to the default value for now
thevalue = fZout;
// some counting
ntris_tried = 0;
// no point in proceeding if xx or yy are silly
if ((xx>fXNmax) || (xx<fXNmin) || (yy>fYNmax) || (yy<fYNmin)) return thevalue;
// check existing Delaunay triangles for a good one
for (it=1; it<=fNdt; it++) {
p = fPTried[it-1];
n = fNTried[it-1];
m = fMTried[it-1];
// p, n and m form a previously found Delaunay triangle, does it
// enclose the point?
if (Enclose(p,n,m,0)) {
// yes, we have the triangle
thevalue = InterpolateOnPlane(p,n,m,0);
return thevalue;
}
}
// is this point inside the convex hull?
shouldbein = InHull(0,-1);
if (!shouldbein) return thevalue;
// it must be in a Delaunay triangle - find it...
// order mass points by distance in mass plane from desired point
for (it=1; it<=fNpoints; it++) {
vxN = fXN[it];
vyN = fYN[it];
fDist[it-1] = TMath::Sqrt((xx-vxN)*(xx-vxN)+(yy-vyN)*(yy-vyN));
}
// sort array 'fDist' to find closest points
TMath::Sort(fNpoints, fDist, fOrder, kFALSE);
for (it=0; it<fNpoints; it++) fOrder[it]++;
// loop over triplets of close points to try to find a triangle that
// encloses the point.
for (k=3; k<=fNpoints; k++) {
m = fOrder[k-1];
for (j=2; j<=k-1; j++) {
n = fOrder[j-1];
for (i=1; i<=j-1; i++) {
p = fOrder[i-1];
if (ntris_tried > fMaxIter) {
// perhaps this point isn't in the hull after all
/// Warning("Interpolate",
/// "Abandoning the effort to find a Delaunay triangle (and thus interpolated z-value) for point %g %g"
/// ,xx,yy);
return thevalue;
}
ntris_tried++;
// check the points aren't colinear
d1 = TMath::Sqrt((fXN[p]-fXN[n])*(fXN[p]-fXN[n])+(fYN[p]-fYN[n])*(fYN[p]-fYN[n]));
d2 = TMath::Sqrt((fXN[p]-fXN[m])*(fXN[p]-fXN[m])+(fYN[p]-fYN[m])*(fYN[p]-fYN[m]));
d3 = TMath::Sqrt((fXN[n]-fXN[m])*(fXN[n]-fXN[m])+(fYN[n]-fYN[m])*(fYN[n]-fYN[m]));
if ((d1+d2<=d3) || (d1+d3<=d2) || (d2+d3<=d1)) goto L90;
// does the triangle enclose the point?
if (!Enclose(p,n,m,0)) goto L90;
// is it a Delaunay triangle? (ie. are there any other points
// inside the circle that is defined by its vertices?)
// test the triangle for Delaunay'ness
// loop over all other points testing each to see if it's
// inside the triangle's circle
ndegen = 0;
for ( z=1; z<=fNpoints; z++) {
if ((z==p) || (z==n) || (z==m)) goto L50;
// An easy first check is to see if point z is inside the triangle
// (if it's in the triangle it's also in the circle)
// point z cannot be inside the triangle if it's further from (xx,yy)
// than the furthest pointing making up the triangle - test this
for (l=1; l<=fNpoints; l++) {
if (fOrder[l-1] == z) {
if ((l<i) || (l<j) || (l<k)) {
// point z is nearer to (xx,yy) than m, n or p - it could be in the
// triangle so call enclose to find out
// if it is inside the triangle this can't be a Delaunay triangle
if (Enclose(p,n,m,z)) goto L90;
} else {
// there's no way it could be in the triangle so there's no point
// calling enclose
goto L1;
}
}
}
// is point z colinear with any pair of the triangle points?
L1:
if (((fXN[p]-fXN[z])*(fYN[p]-fYN[n])) == ((fYN[p]-fYN[z])*(fXN[p]-fXN[n]))) {
// z is colinear with p and n
a = p;
b = n;
} else if (((fXN[p]-fXN[z])*(fYN[p]-fYN[m])) == ((fYN[p]-fYN[z])*(fXN[p]-fXN[m]))) {
// z is colinear with p and m
a = p;
b = m;
} else if (((fXN[n]-fXN[z])*(fYN[n]-fYN[m])) == ((fYN[n]-fYN[z])*(fXN[n]-fXN[m]))) {
// z is colinear with n and m
a = n;
b = m;
} else {
a = 0;
b = 0;
}
if (a != 0) {
// point z is colinear with 2 of the triangle points, if it lies
// between them it's in the circle otherwise it's outside
if (fXN[a] != fXN[b]) {
if (((fXN[z]-fXN[a])*(fXN[z]-fXN[b])) < 0) {
goto L90;
} else if (((fXN[z]-fXN[a])*(fXN[z]-fXN[b])) == 0) {
// At least two points are sitting on top of each other, we will
// treat these as one and not consider this a 'multiple points lying
// on a common circle' situation. It is a sign something could be
// wrong though, especially if the two coincident points have
// different fZ's. If they don't then this is harmless.
Warning("Interpolate", "Two of these three points are coincident %d %d %d",a,b,z);
}
} else {
if (((fYN[z]-fYN[a])*(fYN[z]-fYN[b])) < 0) {
goto L90;
} else if (((fYN[z]-fYN[a])*(fYN[z]-fYN[b])) == 0) {
// At least two points are sitting on top of each other - see above.
Warning("Interpolate", "Two of these three points are coincident %d %d %d",a,b,z);
}
}
// point is outside the circle, move to next point
goto L50;
}
// if point z were to look at the triangle, which point would it see
// lying between the other two? (we're going to form a quadrilateral
// from the points, and then demand certain properties of that
// quadrilateral)
dxz[0] = fXN[p]-fXN[z];
dyz[0] = fYN[p]-fYN[z];
dxz[1] = fXN[n]-fXN[z];
dyz[1] = fYN[n]-fYN[z];
dxz[2] = fXN[m]-fXN[z];
dyz[2] = fYN[m]-fYN[z];
for(l=1; l<=3; l++) {
dx1 = dxz[l-1];
dx2 = dxz[l%3];
dx3 = dxz[(l+1)%3];
dy1 = dyz[l-1];
dy2 = dyz[l%3];
dy3 = dyz[(l+1)%3];
u = (dy3*dx2-dx3*dy2)/(dy1*dx2-dx1*dy2);
v = (dy3*dx1-dx3*dy1)/(dy2*dx1-dx2*dy1);
if ((u>=0) && (v>=0)) {
// vector (dx3,dy3) is expressible as a sum of the other two vectors
// with positive coefficents -> i.e. it lies between the other two vectors
if (l == 1) {
f = m;
o1 = p;
o2 = n;
} else if (l == 2) {
f = p;
o1 = n;
o2 = m;
} else {
f = n;
o1 = m;
o2 = p;
}
goto L2;
}
}
/// Error("Interpolate", "Should not get to here");
// may as well soldier on
f = m;
o1 = p;
o2 = n;
L2:
// this is not a valid quadrilateral if the diagonals don't cross,
// check that points f and z lie on opposite side of the line o1-o2,
// this is true if the angle f-o1-z is greater than o2-o1-z and o2-o1-f
cfo1k = ((fXN[f]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[f]-fYN[o1])*(fYN[z]-fYN[o1]))/
TMath::Sqrt(((fXN[f]-fXN[o1])*(fXN[f]-fXN[o1])+(fYN[f]-fYN[o1])*(fYN[f]-fYN[o1]))*
((fXN[z]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[z]-fYN[o1])*(fYN[z]-fYN[o1])));
co2o1k = ((fXN[o2]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[z]-fYN[o1]))/
TMath::Sqrt(((fXN[o2]-fXN[o1])*(fXN[o2]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[o2]-fYN[o1]))*
((fXN[z]-fXN[o1])*(fXN[z]-fXN[o1]) + (fYN[z]-fYN[o1])*(fYN[z]-fYN[o1])));
co2o1f = ((fXN[o2]-fXN[o1])*(fXN[f]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[f]-fYN[o1]))/
TMath::Sqrt(((fXN[o2]-fXN[o1])*(fXN[o2]-fXN[o1])+(fYN[o2]-fYN[o1])*(fYN[o2]-fYN[o1]))*
((fXN[f]-fXN[o1])*(fXN[f]-fXN[o1]) + (fYN[f]-fYN[o1])*(fYN[f]-fYN[o1]) ));
if ((cfo1k>co2o1k) || (cfo1k>co2o1f)) {
// not a valid quadrilateral - point z is definitely outside the circle
goto L50;
}
// calculate the 2 internal angles of the quadrangle formed by joining
// points z and f to points o1 and o2, at z and f. If they sum to less
// than 180 degrees then z lies outside the circle
dko1 = TMath::Sqrt((fXN[z]-fXN[o1])*(fXN[z]-fXN[o1])+(fYN[z]-fYN[o1])*(fYN[z]-fYN[o1]));
dko2 = TMath::Sqrt((fXN[z]-fXN[o2])*(fXN[z]-fXN[o2])+(fYN[z]-fYN[o2])*(fYN[z]-fYN[o2]));
dfo1 = TMath::Sqrt((fXN[f]-fXN[o1])*(fXN[f]-fXN[o1])+(fYN[f]-fYN[o1])*(fYN[f]-fYN[o1]));
dfo2 = TMath::Sqrt((fXN[f]-fXN[o2])*(fXN[f]-fXN[o2])+(fYN[f]-fYN[o2])*(fYN[f]-fYN[o2]));
c1 = ((fXN[z]-fXN[o1])*(fXN[z]-fXN[o2])+(fYN[z]-fYN[o1])*(fYN[z]-fYN[o2]))/dko1/dko2;
c2 = ((fXN[f]-fXN[o1])*(fXN[f]-fXN[o2])+(fYN[f]-fYN[o1])*(fYN[f]-fYN[o2]))/dfo1/dfo2;
sin_sum = c1*TMath::Sqrt(1-c2*c2)+c2*TMath::Sqrt(1-c1*c1);
// sin_sum doesn't always come out as zero when it should do.
if (sin_sum < -1.E-6) {
// z is inside the circle, this is not a Delaunay triangle
goto L90;
} else if (TMath::Abs(sin_sum) <= 1.E-6) {
// point z lies on the circumference of the circle (within rounding errors)
// defined by the triangle, so there is potential for degeneracy in the
// triangle set (Delaunay triangulation does not give a unique way to split
// a polygon whose points lie on a circle into constituent triangles). Make
// a note of the additional point number.
ndegen++;
degen = z;
fdegen = f;
o1degen = o1;
o2degen = o2;
}
L50:
continue;
}
// This is a good triangle
if (ndegen > 0) {
// but is degenerate with at least one other,
// haven't figured out what to do if more than 4 points are involved
/// if (ndegen > 1) {
/// Error("Interpolate",
/// "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",
/// p,n,m,degen);
/// return thevalue;
/// }
// we have a quadrilateral which can be split down either diagonal
// (d<->f or o1<->o2) to form valid Delaunay triangles. Choose diagonal
// with highest average z-value. Whichever we choose we will have
// verified two triangles as good and two as bad, only note the good ones
d = degen;
f = fdegen;
o1 = o1degen;
o2 = o2degen;
if ((fZ[o1-1]+fZ[o2-1]) > (fZ[d-1]+fZ[f-1])) {
// best diagonalisation of quadrilateral is current one, we have
// the triangle
t1 = p;
t2 = n;
t3 = m;
// file the good triangles
FileIt(p, n, m);
FileIt(d, o1, o2);
} else {
// use other diagonal to split quadrilateral, use triangle formed by
// point f, the degnerate point d and whichever of o1 and o2 create
// an enclosing triangle
t1 = f;
t2 = d;
if (Enclose(f,d,o1,0)) {
t3 = o1;
} else {
t3 = o2;
}
// file the good triangles
FileIt(f, d, o1);
FileIt(f, d, o2);
}
} else {
// this is a Delaunay triangle, file it
FileIt(p, n, m);
t1 = p;
t2 = n;
t3 = m;
}
// do the interpolation
thevalue = InterpolateOnPlane(t1,t2,t3,0);
return thevalue;
L90:
continue;
}
}
}
if (shouldbein) {
Error("Interpolate",
"Point outside hull when expected inside: this point could be dodgy %g %g %d",
xx, yy, ntris_tried);
}
return thevalue;
}
//______________________________________________________________________________
void TGraphDelaunay::SetMaxIter(Int_t n)
{
// Defines the number of triangles tested for a Delaunay triangle
// (number of iterations) before abandoning the search
fAllTri = kFALSE;
fMaxIter = n;
}
//______________________________________________________________________________
void TGraphDelaunay::SetMarginBinsContent(Double_t z)
{
// Sets the histogram bin height for points lying outside the convex hull ie:
// the bins in the margin.
fZout = z;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.