diff options
Diffstat (limited to 'Build/source/utils/asymptote/Delaunay.cc')
-rw-r--r-- | Build/source/utils/asymptote/Delaunay.cc | 212 |
1 files changed, 212 insertions, 0 deletions
diff --git a/Build/source/utils/asymptote/Delaunay.cc b/Build/source/utils/asymptote/Delaunay.cc new file mode 100644 index 00000000000..3e9e5507f21 --- /dev/null +++ b/Build/source/utils/asymptote/Delaunay.cc @@ -0,0 +1,212 @@ +// Robust version of Gilles Dumoulin's C++ port of Paul Bourke's +// triangulation code available from +// http://astronomy.swin.edu.au/~pbourke/papers/triangulate +// Used with permission of Paul Bourke. +// Segmentation fault and numerical robustness improvements by John C. Bowman + +#include <cassert> +#include "Delaunay.h" +#include "predicates.h" + +inline double max(double a, double b) +{ + return (a > b) ? a : b; +} + +int XYZCompare(const void *v1, const void *v2) +{ + double x1=((XYZ*)v1)->p[0]; + double x2=((XYZ*)v2)->p[0]; + if(x1 < x2) + return(-1); + else if(x1 > x2) + return(1); + else + return(0); +} + +inline double hypot2(double *x) +{ + return x[0]*x[0]+x[1]*x[1]; +} + +/////////////////////////////////////////////////////////////////////////////// +// Triangulate(): +// Triangulation subroutine +// Takes as input NV vertices in array pxyz +// Returned is a list of ntri triangular faces in the array v +// These triangles are arranged in a consistent clockwise order. +// The triangle array v should be allocated to 4 * nv +// The vertex array pxyz must be big enough to hold 3 additional points. +// By default, the array pxyz is automatically presorted and postsorted. +/////////////////////////////////////////////////////////////////////////////// + +Int Triangulate(Int nv, XYZ pxyz[], ITRIANGLE v[], Int &ntri, + bool presort, bool postsort) +{ + Int emax = 200; + + if(presort) qsort(pxyz,nv,sizeof(XYZ),XYZCompare); + else postsort=false; + +/* Allocate memory for the completeness list, flag for each triangle */ + Int trimax = 4 * nv; + Int *complete = new Int[trimax]; +/* Allocate memory for the edge list */ + IEDGE *edges = new IEDGE[emax]; +/* + Find the maximum and minimum vertex bounds. + This is to allow calculation of the bounding triangle +*/ + double xmin = pxyz[0].p[0]; + double ymin = pxyz[0].p[1]; + double xmax = xmin; + double ymax = ymin; + for(Int i = 1; i < nv; i++) { + XYZ *pxyzi=pxyz+i; + double x=pxyzi->p[0]; + double y=pxyzi->p[1]; + if (x < xmin) xmin = x; + if (x > xmax) xmax = x; + if (y < ymin) ymin = y; + if (y > ymax) ymax = y; + } + double dx = xmax - xmin; + double dy = ymax - ymin; +/* + Set up the supertriangle. + This is a triangle which encompasses all the sample points. + The supertriangle coordinates are added to the end of the + vertex list. The supertriangle is the first triangle in + the triangle list. +*/ + static const double margin=0.01; + double xmargin=margin*dx; + double ymargin=margin*dy; + pxyz[nv+0].p[0] = xmin-xmargin; + pxyz[nv+0].p[1] = ymin-ymargin; + pxyz[nv+1].p[0] = xmin-xmargin; + pxyz[nv+1].p[1] = ymax+ymargin+dx; + pxyz[nv+2].p[0] = xmax+xmargin+dy; + pxyz[nv+2].p[1] = ymin-ymargin; + v->p1 = nv; + v->p2 = nv+1; + v->p3 = nv+2; + complete[0] = false; + ntri = 1; +/* + Include each point one at a time into the existing mesh +*/ + for(Int i = 0; i < nv; i++) { + Int nedge = 0; + double *d=pxyz[i].p; +/* + Set up the edge buffer. + If the point d lies inside the circumcircle then the + three edges of that triangle are added to the edge buffer + and that triangle is removed. +*/ + for(Int j = 0; j < ntri; j++) { + if(complete[j]) + continue; + ITRIANGLE *vj=v+j; + + double *a=pxyz[vj->p1].p; + double *b=pxyz[vj->p2].p; + double *c=pxyz[vj->p3].p; + + if(incircle(a,b,c,d) <= 0.0) { // Point d is inside or on circumcircle +/* Check that we haven't exceeded the edge list size */ + if(nedge + 3 >= emax) { + emax += 100; + IEDGE *p_EdgeTemp = new IEDGE[emax]; + for (Int i = 0; i < nedge; i++) { + p_EdgeTemp[i] = edges[i]; + } + delete[] edges; + edges = p_EdgeTemp; + } + ITRIANGLE *vj=v+j; + Int p1=vj->p1; + Int p2=vj->p2; + Int p3=vj->p3; + edges[nedge].p1 = p1; + edges[nedge].p2 = p2; + edges[++nedge].p1 = p2; + edges[nedge].p2 = p3; + edges[++nedge].p1 = p3; + edges[nedge].p2 = p1; + ++nedge; + ntri--; + v[j] = v[ntri]; + complete[j] = complete[ntri]; + j--; + } else { + double A=hypot2(a); + double B=hypot2(b); + double C=hypot2(c); + double a0=orient2d(a,b,c); + // Is d[0] > xc+r for circumcircle abc of radius r about (xc,yc)? + if(d[0]*a0 < 0.5*orient2d(A,a[1],B,b[1],C,c[1])) + complete[j]= + incircle(a[0]*a0,a[1]*a0,b[0]*a0,b[1]*a0,c[0]*a0,c[1]*a0, + d[0]*a0,0.5*orient2d(a[0],A,b[0],B,c[0],C)) > 0.0; + } + } +/* + Tag multiple edges + Note: if all triangles are specified anticlockwise then all + interior edges are opposite pointing in direction. +*/ + for(Int j = 0; j < nedge - 1; j++) { + for(Int k = j + 1; k < nedge; k++) { + if((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) { + edges[j].p1 = -1; + edges[j].p2 = -1; + edges[k].p1 = -1; + edges[k].p2 = -1; + } + } + } +/* + Form new triangles for the current point + Skipping over any tagged edges. + All edges are arranged in clockwise order. +*/ + for(Int j = 0; j < nedge; j++) { + if(edges[j].p1 < 0 || edges[j].p2 < 0) + continue; + v[ntri].p1 = edges[j].p1; + v[ntri].p2 = edges[j].p2; + v[ntri].p3 = i; + complete[ntri] = false; + ntri++; + assert(ntri < trimax); + } + } +/* + Remove triangles with supertriangle vertices + These are triangles which have a vertex number greater than nv +*/ + for(Int i = 0; i < ntri; i++) { + ITRIANGLE *vi=v+i; + if(vi->p1 >= nv || vi->p2 >= nv || vi->p3 >= nv) { + ntri--; + *vi = v[ntri]; + i--; + } + } + delete[] edges; + delete[] complete; + + if(postsort) { + for(Int i = 0; i < ntri; i++) { + ITRIANGLE *vi=v+i; + vi->p1=pxyz[vi->p1].i; + vi->p2=pxyz[vi->p2].i; + vi->p3=pxyz[vi->p3].i; + } + } + + return 0; +} |