Delaunay.cc

// 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;
}