/* dt.c - Don Yang (uguu.org) Delaunay triangulation. 07/14/01 */ #include #include #include #include #include"dt.h" #include"grid.h" #include"triangle.h" /* Local functions */ static int AddPoint(Triangle **triangles, Vertex *point, int update); static int Degenerate(Triangle **triangles, Vertex *point, Triangle *s, int update); static Triangle *FindEncloseTriangle(Triangle *root, Vertex *point); static Triangle *FindOppositeTriangle(Triangle *current, Vertex *point, Vertex **opposite); static int Repartition(Triangle **triangles, Triangle **current, Vertex *point, int update); /************************************************************* Triangulate */ int Triangulate(Triangle **triangles, Vertex **points, Grid *grid) { Vertex *slist, *v[4], *p; int i; assert(triangles); assert(points); assert(grid); /* Create supertriangles 3 2 0 1 */ slist = NULL; if( AddVertex(&slist, &v[0], grid->minx - 16, grid->miny - 16) || AddVertex(&slist, &v[1], grid->maxx + 16, grid->miny - 16) || AddVertex(&slist, &v[2], grid->maxx + 16, grid->maxy + 16) ) return 1; if( AddTriangle(triangles, NULL, v) ) return 1; v[3] = v[1]; if( AddVertex(&slist, &v[1], grid->minx - 16, grid->maxy + 16) ) return 1; if( AddTriangle(triangles, NULL, v) ) return 1; /* Add points */ i = 1; for(p = *points; p; p = p->next) { if( AddPoint(triangles, p, i++) ) return 1; } /* Remove triangles sharing same edge as the supertriangle i.e. if one of the four supertriangle vertices is used */ for(i = 0; i < 4; i++) { while( v[i]->tlist ) DeleteTriangle(triangles, v[i]->tlist->t); } /* Delete supertriangle vertices */ DeleteVertexList(&slist); /* Success */ return 0; } /* Triangulate() */ /**************************************************************** AddPoint */ static int AddPoint(Triangle **triangles, Vertex *point, int update) { Triangle *s, *div[3]; Vertex *sv[3], *v[3]; int i, j; assert(triangles); /* Find enclosing triangle */ s = FindEncloseTriangle(*triangles, point); if( s == NULL ) { fprintf(stderr, "No enclosing triangle for (%g, %g)\n", point->x, point->y); return 0; } /* Break triangle */ sv[0] = s->v[0]; sv[1] = s->v[1]; sv[2] = s->v[2]; for(i = 0; i < 3; i++) { for(j = 0; j < 3; j++) v[j] = i == j ? point : sv[j]; if( AddTriangle(triangles, &div[i], v) ) return 1; if( div[i] == NULL ) { /* Degenerate found */ for(j = 0; j < i; j++) DeleteTriangle(triangles, div[j]); return Degenerate(triangles, point, s, update); } } /* Triangle inserted properly, remove original */ DeleteTriangle(triangles, s); /* Repartition triangles */ for(i = 0; i < 3; i++) { if( Repartition(triangles, &div[i], point, update) ) return 1; } return 0; } /* AddPoint() */ /************************************************************** Degenerate */ static int Degenerate(Triangle **triangles, Vertex *point, Triangle *s, int update) { Triangle *div[4], *otriangle; Vertex *sv[2], *v[3], *overtex; double dx[3], dy[3]; int i, j, x; /* Silently drop points that are too close to any of the vertices */ for(i = 0; i < 3; i++) { dx[i] = s->v[i]->x - point->x; dy[i] = s->v[i]->y - point->y; if( dx[i] * dx[i] + dy[i] * dy[i] < MIN_VERTEX_DISTANCE_2 ) return 0; } /* Find colinear edge */ if( fabs(dx[0] * dy[1] - dx[1] * dy[0]) < COLINEAR_EPSILON ) x = 2; else if( fabs(dx[0] * dy[2] - dx[2] * dy[0]) < COLINEAR_EPSILON ) x = 1; else x = 0; /* Break triangle on colinear edge */ if( (otriangle = FindOppositeTriangle(s, s->v[x], &overtex)) == NULL ) return 0; assert(overtex); for(i = j = 0; i < 3; i++) { if( i != x ) sv[j++] = s->v[i]; } for(i = 0; i < 2; i++) { v[0] = s->v[x]; v[1] = sv[i]; v[2] = point; if( AddTriangle(triangles, &div[i], v) ) return 1; v[0] = overtex; if( AddTriangle(triangles, &div[i + 2], v) ) return 1; } if( !div[0] || !div[1] || !div[2] || !div[3] ) { for(i = 0; i < 4; i++) { if( div[i] ) DeleteTriangle(triangles, div[i]); } fprintf(stderr, "Can not insert point (%g, %g)\n" "Contained in: (%g, %g), (%g, %g), (%g, %g), (%g, %g)\n", point->x, point->y, s->v[0]->x, s->v[0]->y, s->v[1]->x, s->v[1]->y, s->v[2]->x, s->v[2]->y, overtex->x, overtex->y); return 0; } /* Delete original triangles */ DeleteTriangle(triangles, s); DeleteTriangle(triangles, otriangle); /* Repartition triangles */ for(i = 0; i < 4; i++) { if( Repartition(triangles, &div[i], point, update) ) return 1; } /* Success */ return 0; } /* Degenerate() */ /***************************************************** FindEncloseTriangle */ static Triangle *FindEncloseTriangle(Triangle *root, Vertex *point) { Triangle *ret; while( root != NULL ) { /* Check if current triangle is the enclosing one */ if( Enclosed(root->v, point) ) return root; /* Find enclosing triangle in one of the subtrees. Potentially it could search through all the triangles, but the algorithm will try the triangles with closer X first. */ if( point->x > root->key ) { if( (ret = FindEncloseTriangle(root->right, point)) != NULL ) return ret; /* Instead of recurse through left subtree, apply tail recursion */ root = root->left; } else { if( (ret = FindEncloseTriangle(root->left, point)) != NULL ) return ret; root = root->right; } } return NULL; } /* FindEncloseTriangle() */ /**************************************************** FindOppositeTriangle */ static Triangle *FindOppositeTriangle(Triangle *current, Vertex *point, Vertex **opposite) { Vertex *shared[2]; Triangle *t; TList *i, *j; assert(opposite); *opposite = NULL; /* Find shared vertices */ if( point == current->v[0] ) { shared[0] = current->v[1]; shared[1] = current->v[2]; } else if( point == current->v[1] ) { shared[0] = current->v[0]; shared[1] = current->v[2]; } else { shared[0] = current->v[0]; shared[1] = current->v[1]; } /* Find triangle using these vertices that is not current triangle. Guaranteed to be at most one, because there can only be at most one triangle on each side of the shared edge (triangles can't overlap), and one of them is current. */ for(i = shared[0]->tlist; i; i = i->next) { if( i->t == current ) continue; for(j = shared[1]->tlist; j; j = j->next) { if( i->t == j->t ) { t = (Triangle*)(j->t); /* Find opposite vertex */ if( t->v[0] != shared[0] && t->v[0] != shared[1] ) *opposite = t->v[0]; if( t->v[1] != shared[0] && t->v[1] != shared[1] ) *opposite = t->v[1]; if( t->v[2] != shared[0] && t->v[2] != shared[1] ) *opposite = t->v[2]; return t; } } } return NULL; } /* FindOppositeTriangle() */ /************************************************************* Repartition */ static int Repartition(Triangle **triangles, Triangle **current, Vertex *point, int update) { Triangle *ntriangle[2], *otriangle; Vertex *svertex[2], *nvertex[3], *overtex; double dx, dy; int i, j; assert(triangles); assert(current); if( *current == NULL ) return 0; /* Do nothing if current triangle has already been updated in this generation. */ if( (*current)->id >= update ) return 0; (*current)->id = update; /* Find vertex opposite of current triangle */ if( (otriangle = FindOppositeTriangle(*current, point, &overtex)) == NULL ) return 0; assert(overtex); /* Repartition triangle if opposite point is contained in current triangle's circumcircle. */ dx = overtex->x - (*current)->cx; dy = overtex->y - (*current)->cy; if( dx * dx + dy * dy >= (*current)->cr2 ) return 0; /* Find shared vertices */ for(i = j = 0; i < 3; i++) { if( (*current)->v[i] != point ) svertex[j++] = (*current)->v[i]; } /* Repartition */ nvertex[0] = point; nvertex[1] = overtex; for(i = 0; i < 2; i++) { nvertex[2] = svertex[i]; if( AddTriangle(triangles, &ntriangle[i], nvertex) ) return 1; if( ntriangle[i] == NULL ) { if( i == 1 ) DeleteTriangle(triangles, ntriangle[0]); return 0; } } DeleteTriangle(triangles, *current); DeleteTriangle(triangles, otriangle); *current = NULL; /* Repartition neighboring triangles */ if( Repartition(triangles, &ntriangle[0], point, update) || Repartition(triangles, &ntriangle[1], point, update) || Repartition(triangles, &ntriangle[0], overtex, update) || Repartition(triangles, &ntriangle[1], overtex, update) ) return 1; /* Success */ return 0; } /* Repartition() */