2 Efficient Triangulation Algorithm Suitable for Terrain Modelling
4 An Algorithm for Interpolating Irregularly-Spaced Data
5 with Applications in Terrain Modelling
8 Presented at Pan Pacific Computer Conference, Beijing, China.
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.
18 http://paulbourke.net/papers/triangulate/
19 http://paulbourke.net/papers/triangulate/triangulate.c
33 #define EPSILON 0.000001
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
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)
46 double m1,m2,mx1,mx2,my1,my2;
48 double fabsy1y2 = fabs(y1-y2);
49 double fabsy2y3 = fabs(y2-y3);
51 /* Check for coincident points */
52 if (fabsy1y2 < EPSILON && fabsy2y3 < EPSILON)
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;
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;
78 *yc = m2 * (*xc - mx2) + my2;
84 *rsqr = dx*dx + dy*dy;
88 drsqr = dx*dx + dy*dy;
91 return((drsqr <= *rsqr) ? TRUE : FALSE);
92 Proposed by Chuck Morris */
93 return((drsqr - *rsqr) <= EPSILON ? TRUE : FALSE);
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);
108 delaunay (int nv,XYZ *pxyz,ITRIANGLE *v,int *ntri)
110 int *complete = NULL;
113 int trimax,emax = 200;
118 double xp,yp,x1,y1,x2,y2,x3,y3,xc=0,yc=0,r=0;
119 double xmin,xmax,ymin,ymax,xmid,ymid;
122 /* Allocate memory for the completeness list, flag for each triangle */
124 if ((complete = malloc(trimax*sizeof(int))) == NULL) {
129 /* Allocate memory for the edge list */
130 if ((edges = malloc(emax*(long)sizeof(IEDGE))) == NULL) {
136 Find the maximum and minimum vertex bounds.
137 This is to allow calculation of the bounding triangle
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;
151 dmax = (dx > dy) ? dx : dy;
152 xmid = (xmax + xmin) / 2.0;
153 ymid = (ymax + ymin) / 2.0;
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
162 pxyz[nv+0].x = xmid - 20 * dmax;
163 pxyz[nv+0].y = ymid - dmax;
166 pxyz[nv+1].y = ymid + 20 * dmax;
168 pxyz[nv+2].x = xmid + 20 * dmax;
169 pxyz[nv+2].y = ymid - dmax;
178 Include each point one at a time into the existing mesh
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.
192 for (j=0;j<(*ntri);j++) {
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)
205 /* Check that we haven't exceeded the edge list size */
206 if (nedge+3 >= emax) {
208 if ((edges = realloc(edges,emax*(long)sizeof(IEDGE))) == NULL) {
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;
221 complete[j] = complete[(*ntri)-1];
229 Note: if all triangles are specified anticlockwise then all
230 interior edges are opposite pointing in direction.
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)) {
240 /* Shouldn't need the following, see note above */
241 if ((edges[j].p1 == edges[k].p1) && (edges[j].p2 == edges[k].p2)) {
251 Form new triangles for the current point
252 Skipping over any tagged edges.
253 All edges are arranged in clockwise order.
255 for (j=0;j<nedge;j++) {
256 if (edges[j].p1 < 0 || edges[j].p2 < 0)
258 if ((*ntri) >= trimax) {
262 v[*ntri].p1 = edges[j].p1;
263 v[*ntri].p2 = edges[j].p2;
265 complete[*ntri] = FALSE;
271 Remove triangles with supertriangle vertices
272 These are triangles which have a vertex number greater than nv
274 for (i=0;i<(*ntri);i++) {
275 if (v[i].p1 >= nv || v[i].p2 >= nv || v[i].p3 >= nv) {
290 delaunay_xyzcompare (const void *v1, const void *v2)
297 else if (p1->x > p2->x)