summaryrefslogtreecommitdiff
path: root/Build/source/utils/asymptote/Delaunay.cc
diff options
context:
space:
mode:
Diffstat (limited to 'Build/source/utils/asymptote/Delaunay.cc')
-rw-r--r--Build/source/utils/asymptote/Delaunay.cc212
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..2090ddfeb3e
--- /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;
+}