git.hungrycats.org Git - xscreensaver/blob - hacks/delaunay.c

   1 /* Triangulate
   2    Efficient Triangulation Algorithm Suitable for Terrain Modelling
   3    or
   4    An Algorithm for Interpolating Irregularly-Spaced Data
   5    with Applications in Terrain Modelling
   6
   7    Written by Paul Bourke
   8    Presented at Pan Pacific Computer Conference, Beijing, China.
   9    January 1989
  10    Abstract
  11
  12    A discussion of a method that has been used with success in terrain
  13    modelling to estimate the height at any point on the land surface
  14    from irregularly distributed samples. The special requirements of
  15    terrain modelling are discussed as well as a detailed description
  16    of the algorithm and an example of its application.
  17
  18    http://paulbourke.net/papers/triangulate/
  19    http://paulbourke.net/papers/triangulate/triangulate.c
  20  */
  21
  22 #include <stdlib.h>
  23 #include <math.h>
  24
  25 #include "delaunay.h"
  26
  27 typedef struct {
  28    int p1,p2;
  29 } IEDGE;
  30
  31 #define TRUE 1
  32 #define FALSE 0
  33 #define EPSILON 0.000001
  34
  35 /*
  36    Return TRUE if a point (xp,yp) is inside the circumcircle made up
  37    of the points (x1,y1), (x2,y2), (x3,y3)
  38    The circumcircle centre is returned in (xc,yc) and the radius r
  39    NOTE: A point on the edge is inside the circumcircle
  40 */
  41 static int
  42 circumcircle (double xp,double yp,
  43               double x1,double y1,double x2,double y2,double x3,double y3,
  44               double *xc,double *yc,double *rsqr)
  45 {
  46   double m1,m2,mx1,mx2,my1,my2;
  47   double dx,dy,drsqr;
  48   double fabsy1y2 = fabs(y1-y2);
  49   double fabsy2y3 = fabs(y2-y3);
  50
  51   /* Check for coincident points */
  52   if (fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
  53     return(FALSE);
  54
  55   if (fabsy1y2 < EPSILON) {
  56     m2 = - (x3-x2) / (y3-y2);
  57     mx2 = (x2 + x3) / 2.0;
  58     my2 = (y2 + y3) / 2.0;
  59     *xc = (x2 + x1) / 2.0;
  60     *yc = m2 * (*xc - mx2) + my2;
  61   } else if (fabsy2y3 < EPSILON) {
  62     m1 = - (x2-x1) / (y2-y1);
  63     mx1 = (x1 + x2) / 2.0;
  64     my1 = (y1 + y2) / 2.0;
  65     *xc = (x3 + x2) / 2.0;
  66     *yc = m1 * (*xc - mx1) + my1;
  67   } else {
  68     m1 = - (x2-x1) / (y2-y1);
  69     m2 = - (x3-x2) / (y3-y2);
  70     mx1 = (x1 + x2) / 2.0;
  71     mx2 = (x2 + x3) / 2.0;
  72     my1 = (y1 + y2) / 2.0;
  73     my2 = (y2 + y3) / 2.0;
  74     *xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
  75     if (fabsy1y2 > fabsy2y3) {
  76       *yc = m1 * (*xc - mx1) + my1;
  77     } else {
  78       *yc = m2 * (*xc - mx2) + my2;
  79     }
  80   }
  81
  82   dx = x2 - *xc;
  83   dy = y2 - *yc;
  84   *rsqr = dx*dx + dy*dy;
  85
  86   dx = xp - *xc;
  87   dy = yp - *yc;
  88   drsqr = dx*dx + dy*dy;
  89
  90   // Original
  91   //return((drsqr <= *rsqr) ? TRUE : FALSE);
  92   // Proposed by Chuck Morris
  93   return((drsqr - *rsqr) <= EPSILON ? TRUE : FALSE);
  94 }
  95
  96
  97 /*
  98    Triangulation subroutine
  99    Takes as input NV vertices in array pxyz
 100    Returned is a list of ntri triangular faces in the array v
 101    These triangles are arranged in a consistent clockwise order.
 102    The triangle array 'v' should be malloced to 3 * nv
 103    The vertex array pxyz must be big enough to hold 3 more points
 104    The vertex array must be sorted in increasing x values say
 105    qsort(p,nv,sizeof(XYZ),XYZCompare);
 106 */
 107 int
 108 delaunay (int nv,XYZ *pxyz,ITRIANGLE *v,int *ntri)
 109 {
 110   int *complete = NULL;
 111   IEDGE *edges = NULL;
 112   int nedge = 0;
 113   int trimax,emax = 200;
 114   int status = 0;
 115
 116   int inside;
 117   int i,j,k;
 118   double xp,yp,x1,y1,x2,y2,x3,y3,xc,yc,r;
 119   double xmin,xmax,ymin,ymax,xmid,ymid;
 120   double dx,dy,dmax;
 121
 122   /* Allocate memory for the completeness list, flag for each triangle */
 123   trimax = 4 * nv;
 124   if ((complete = malloc(trimax*sizeof(int))) == NULL) {
 125     status = 1;
 126     goto skip;
 127   }
 128
 129   /* Allocate memory for the edge list */
 130   if ((edges = malloc(emax*(long)sizeof(IEDGE))) == NULL) {
 131     status = 2;
 132     goto skip;
 133   }
 134
 135   /*
 136     Find the maximum and minimum vertex bounds.
 137     This is to allow calculation of the bounding triangle
 138   */
 139   xmin = pxyz[0].x;
 140   ymin = pxyz[0].y;
 141   xmax = xmin;
 142   ymax = ymin;
 143   for (i=1;i<nv;i++) {
 144     if (pxyz[i].x < xmin) xmin = pxyz[i].x;
 145     if (pxyz[i].x > xmax) xmax = pxyz[i].x;
 146     if (pxyz[i].y < ymin) ymin = pxyz[i].y;
 147     if (pxyz[i].y > ymax) ymax = pxyz[i].y;
 148   }
 149   dx = xmax - xmin;
 150   dy = ymax - ymin;
 151   dmax = (dx > dy) ? dx : dy;
 152   xmid = (xmax + xmin) / 2.0;
 153   ymid = (ymax + ymin) / 2.0;
 154
 155   /*
 156     Set up the supertriangle
 157     This is a triangle which encompasses all the sample points.
 158     The supertriangle coordinates are added to the end of the
 159     vertex list. The supertriangle is the first triangle in
 160     the triangle list.
 161   */
 162   pxyz[nv+0].x = xmid - 20 * dmax;
 163   pxyz[nv+0].y = ymid - dmax;
 164   pxyz[nv+0].z = 0.0;
 165   pxyz[nv+1].x = xmid;
 166   pxyz[nv+1].y = ymid + 20 * dmax;
 167   pxyz[nv+1].z = 0.0;
 168   pxyz[nv+2].x = xmid + 20 * dmax;
 169   pxyz[nv+2].y = ymid - dmax;
 170   pxyz[nv+2].z = 0.0;
 171   v[0].p1 = nv;
 172   v[0].p2 = nv+1;
 173   v[0].p3 = nv+2;
 174   complete[0] = FALSE;
 175   *ntri = 1;
 176
 177   /*
 178     Include each point one at a time into the existing mesh
 179   */
 180   for (i=0;i<nv;i++) {
 181
 182     xp = pxyz[i].x;
 183     yp = pxyz[i].y;
 184     nedge = 0;
 185
 186     /*
 187       Set up the edge buffer.
 188       If the point (xp,yp) lies inside the circumcircle then the
 189       three edges of that triangle are added to the edge buffer
 190       and that triangle is removed.
 191     */
 192     for (j=0;j<(*ntri);j++) {
 193       if (complete[j])
 194         continue;
 195       x1 = pxyz[v[j].p1].x;
 196       y1 = pxyz[v[j].p1].y;
 197       x2 = pxyz[v[j].p2].x;
 198       y2 = pxyz[v[j].p2].y;
 199       x3 = pxyz[v[j].p3].x;
 200       y3 = pxyz[v[j].p3].y;
 201       inside = circumcircle(xp,yp,x1,y1,x2,y2,x3,y3,&xc,&yc,&r);
 202       if (xc < xp && ((xp-xc)*(xp-xc)) > r)
 203         complete[j] = TRUE;
 204       if (inside) {
 205         /* Check that we haven't exceeded the edge list size */
 206         if (nedge+3 >= emax) {
 207           emax += 100;
 208           if ((edges = realloc(edges,emax*(long)sizeof(IEDGE))) == NULL) {
 209             status = 3;
 210             goto skip;
 211           }
 212         }
 213         edges[nedge+0].p1 = v[j].p1;
 214         edges[nedge+0].p2 = v[j].p2;
 215         edges[nedge+1].p1 = v[j].p2;
 216         edges[nedge+1].p2 = v[j].p3;
 217         edges[nedge+2].p1 = v[j].p3;
 218         edges[nedge+2].p2 = v[j].p1;
 219         nedge += 3;
 220         v[j] = v[(*ntri)-1];
 221         complete[j] = complete[(*ntri)-1];
 222         (*ntri)--;
 223         j--;
 224       }
 225     }
 226
 227     /*
 228       Tag multiple edges
 229       Note: if all triangles are specified anticlockwise then all
 230       interior edges are opposite pointing in direction.
 231     */
 232     for (j=0;j<nedge-1;j++) {
 233       for (k=j+1;k<nedge;k++) {
 234         if ((edges[j].p1 == edges[k].p2) && (edges[j].p2 == edges[k].p1)) {
 235           edges[j].p1 = -1;
 236           edges[j].p2 = -1;
 237           edges[k].p1 = -1;
 238           edges[k].p2 = -1;
 239         }
 240         /* Shouldn't need the following, see note above */
 241         if ((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)) {
 242           edges[j].p1 = -1;
 243           edges[j].p2 = -1;
 244           edges[k].p1 = -1;
 245           edges[k].p2 = -1;
 246         }
 247       }
 248     }
 249
 250     /*
 251       Form new triangles for the current point
 252       Skipping over any tagged edges.
 253       All edges are arranged in clockwise order.
 254     */
 255     for (j=0;j<nedge;j++) {
 256       if (edges[j].p1 < 0 || edges[j].p2 < 0)
 257         continue;
 258       if ((*ntri) >= trimax) {
 259         status = 4;
 260         goto skip;
 261       }
 262       v[*ntri].p1 = edges[j].p1;
 263       v[*ntri].p2 = edges[j].p2;
 264       v[*ntri].p3 = i;
 265       complete[*ntri] = FALSE;
 266       (*ntri)++;
 267     }
 268   }
 269
 270   /*
 271     Remove triangles with supertriangle vertices
 272     These are triangles which have a vertex number greater than nv
 273   */
 274   for (i=0;i<(*ntri);i++) {
 275     if (v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv) {
 276       v[i] = v[(*ntri)-1];
 277       (*ntri)--;
 278       i--;
 279     }
 280   }
 281
 282  skip:
 283   free(edges);
 284   free(complete);
 285   return(status);
 286 }
 287
 288
 289 int
 290 delaunay_xyzcompare (const void *v1, const void *v2)
 291 {
 292   const XYZ *p1,*p2;
 293   p1 = v1;
 294   p2 = v2;
 295   if (p1->x < p2->x)
 296     return(-1);
 297   else if (p1->x > p2->x)
 298     return(1);
 299   else
 300     return(0);
 301 }