Initial implementation of Constrained Delaunay Refinement

Include/tesselator.h

@@ -107,11 +107,16 @@ enum TessWindingRule
 //         glEnd();
 //     }
 //
+// TESS_CONSTRAINED_DELAUNAY_TRIANGLES
+//   Similar to TESS_POLYGONS, but we output only triangles and we attempt to provide a valid
+//   Constrained Delaunay triangulation.
+
 enum TessElementType
 {
 	TESS_POLYGONS,
 	TESS_CONNECTED_POLYGONS,
 	TESS_BOUNDARY_CONTOURS,
+	TESS_CONSTRAINED_DELAUNAY_TRIANGLES,
 };
 
 typedef float TESSreal;
@@ -189,7 +194,7 @@ void tessAddContour( TESStesselator *tess, int size, const void* pointer, int st
 //   tess - pointer to tesselator object.
 //   windingRule - winding rules used for tesselation, must be one of TessWindingRule.
 //   elementType - defines the tesselation result element type, must be one of TessElementType.
-//   polySize - defines maximum vertices per polygons if output is polygons.
+//   polySize - defines maximum vertices per polygons if output is polygons. If elementType is TESS_CONSTRAINED_DELAUNAY_TRIANGLES, this parameter is ignored.
 //   vertexSize - defines the number of coordinates in tesselation result vertex, must be 2 or 3.
 //   normal - defines the normal of the input contours, of null the normal is calculated automatically.
 // Returns:

Source/geom.c

@@ -33,6 +33,7 @@
 #include <assert.h>
 #include "mesh.h"
 #include "geom.h"
+#include <math.h>
 
 int tesvertLeq( TESSvertex *u, TESSvertex *v )
 {
@@ -259,3 +260,22 @@ void tesedgeIntersect( TESSvertex *o1, TESSvertex *d1,
 		v->t = Interpolate( z1, o2->t, z2, d2->t );
 	}
 }
+
+/*
+	Calculate the angle between v1-v2 and v1-v0
+ */
+TESSreal calcAngle(TESSvertex *v0, TESSvertex *v1, TESSvertex *v2) {
+	TESSreal a[2] = { v2->s - v1->s, v2->t - v1->t };
+	TESSreal b[2] = { v0->s - v1->s, v0->t - v1->t };
+	return acosf((a[0] * b[0] + a[1] * b[1]) /
+							 (sqrt(a[0] * a[0] + a[1] * a[1]) * sqrt(b[0] * b[0] + b[1] * b[1])));
+}
+
+/*
+	Returns 1 is edge is locally delaunay
+ */
+int tesedgeIsLocallyDelaunay( TESShalfEdge *e )
+{
+	return (calcAngle(e->Lnext->Org, e->Lnext->Lnext->Org, e->Org) +
+					calcAngle(e->Sym->Lnext->Org, e->Sym->Lnext->Lnext->Org, e->Sym->Org)) < (M_PI + 0.01);
+}

Source/geom.h

@@ -59,6 +59,7 @@
 
 #define EdgeGoesLeft(e) VertLeq( (e)->Dst, (e)->Org )
 #define EdgeGoesRight(e) VertLeq( (e)->Org, (e)->Dst )
+#define EdgeIsInternal(e) e->Rface && e->Rface->inside
 
 #define ABS(x) ((x) < 0 ? -(x) : (x))
 #define VertL1dist(u,v) (ABS(u->s - v->s) + ABS(u->t - v->t))
@@ -72,5 +73,6 @@ TESSreal	testransEval( TESSvertex *u, TESSvertex *v, TESSvertex *w );
 TESSreal	testransSign( TESSvertex *u, TESSvertex *v, TESSvertex *w );
 int tesvertCCW( TESSvertex *u, TESSvertex *v, TESSvertex *w );
 void tesedgeIntersect( TESSvertex *o1, TESSvertex *d1, TESSvertex *o2, TESSvertex *d2, TESSvertex *v );
+int tesedgeIsLocallyDelaunay( TESShalfEdge *e );
 
 #endif

Source/mesh.c

@@ -80,6 +80,7 @@ static TESShalfEdge *MakeEdge( TESSmesh* mesh, TESShalfEdge *eNext )
 	e->Lface = NULL;
 	e->winding = 0;
 	e->activeRegion = NULL;
+	e->mark = 0;
 
 	eSym->Sym = e;
 	eSym->Onext = eSym;
@@ -88,6 +89,7 @@ static TESShalfEdge *MakeEdge( TESSmesh* mesh, TESShalfEdge *eNext )
 	eSym->Lface = NULL;
 	eSym->winding = 0;
 	eSym->activeRegion = NULL;
+	eSym->mark = 0;
 
 	return e;
 }
@@ -748,6 +750,85 @@ int tessMeshMergeConvexFaces( TESSmesh *mesh, int maxVertsPerFace )
 	return 1;
 }
 
+#include <stdio.h>
+
+void tessMeshFlipEdge( TESSmesh *mesh, TESShalfEdge *edge )
+{
+	assert(EdgeIsInternal(edge));
+
+	TESShalfEdge *a0 = edge;
+	TESShalfEdge *a1 = a0->Lnext;
+	TESShalfEdge *a2 = a1->Lnext;
+	assert(a2->Lnext == a0);
+	TESShalfEdge *b0 = edge->Sym;
+	TESShalfEdge *b1 = b0->Lnext;
+	TESShalfEdge *b2 = b1->Lnext;
+	assert(b2->Lnext == b0);
+
+	TESSvertex *aOrg = a0->Org;
+	TESSvertex *aOpp = a2->Org;
+	TESSvertex *bOrg = b0->Org;
+	TESSvertex *bOpp = b2->Org;
+
+	TESSface *fa = a0->Lface;
+	TESSface *fb = b0->Lface;
+
+	a0->Org = bOpp;
+	a0->Onext = b1->Sym;
+	b0->Org = aOpp;
+	b0->Onext = a1->Sym;
+	a2->Onext = b0;
+	b2->Onext = a0;
+    b1->Onext = a2->Sym;
+    a1->Onext = b2->Sym;
+
+	a0->Lnext = a2;
+	a2->Lnext = b1;
+	b1->Lnext = a0;
+
+	b0->Lnext = b2;
+	b2->Lnext = a1;
+	a1->Lnext = b0;
+
+	a1->Lface = fb;
+	b1->Lface = fa;
+
+	fa->anEdge = a0;
+	fb->anEdge = b0;
+
+	if (aOrg->anEdge == a0) aOrg->anEdge = b1;
+	if (bOrg->anEdge == b0) bOrg->anEdge = a1;
+
+	assert( a0->Lnext->Onext->Sym == a0 );
+	assert( a0->Onext->Sym->Lnext == a0 );
+	assert( a0->Org->anEdge->Org == a0->Org );
+
+
+	assert( a1->Lnext->Onext->Sym == a1 );
+	assert( a1->Onext->Sym->Lnext == a1 );
+	assert( a1->Org->anEdge->Org == a1->Org );
+
+	assert( a2->Lnext->Onext->Sym == a2 );
+	assert( a2->Onext->Sym->Lnext == a2 );
+	assert( a2->Org->anEdge->Org == a2->Org );
+
+	assert( b0->Lnext->Onext->Sym == b0 );
+	assert( b0->Onext->Sym->Lnext == b0 );
+	assert( b0->Org->anEdge->Org == b0->Org );
+
+	assert( b1->Lnext->Onext->Sym == b1 );
+	assert( b1->Onext->Sym->Lnext == b1 );
+	assert( b1->Org->anEdge->Org == b1->Org );
+
+	assert( b2->Lnext->Onext->Sym == b2 );
+	assert( b2->Onext->Sym->Lnext == b2 );
+	assert( b2->Org->anEdge->Org == b2->Org );
+
+	assert(aOrg->anEdge->Org == aOrg);
+	assert(bOrg->anEdge->Org == bOrg);
+
+	assert(a0->Oprev->Onext->Org == a0->Org);
+}
 
 #ifdef DELETE_BY_ZAPPING
 

Source/mesh.h

@@ -143,6 +143,7 @@ struct TESShalfEdge {
 	ActiveRegion *activeRegion;  /* a region with this upper edge (sweep.c) */
 	int winding;    /* change in winding number when crossing
 						  from the right face to the left face */
+	int mark; /* Used by the Edge Flip algorithm */
 };
 
 #define Rface   Sym->Lface
@@ -155,7 +156,6 @@ struct TESShalfEdge {
 #define Dnext   Rprev->Sym  /* 3 pointers */
 #define Rnext   Oprev->Sym  /* 3 pointers */
 
-
 struct TESSmesh {
 	TESSvertex vHead;      /* dummy header for vertex list */
 	TESSface fHead;      /* dummy header for face list */
@@ -258,6 +258,8 @@ int tessMeshMergeConvexFaces( TESSmesh *mesh, int maxVertsPerFace );
 void tessMeshDeleteMesh( TESSalloc* alloc, TESSmesh *mesh );
 void tessMeshZapFace( TESSmesh *mesh, TESSface *fZap );
 
+void tessMeshFlipEdge( TESSmesh *mesh, TESShalfEdge *edge );
+
 #ifdef NDEBUG
 #define tessMeshCheckMesh( mesh )
 #else

Source/tess.c

@@ -365,7 +365,6 @@ int tessMeshTessellateMonoRegion( TESSmesh *mesh, TESSface *face )
 	return 1;
 }
 
-
 /* tessMeshTessellateInterior( mesh ) tessellates each region of
 * the mesh which is marked "inside" the polygon.  Each such region
 * must be monotone.
@@ -382,6 +381,100 @@ int tessMeshTessellateInterior( TESSmesh *mesh )
 			if ( !tessMeshTessellateMonoRegion( mesh, f ) ) return 0;
 		}
 	}
+	return 1;
+}
+
+
+struct EdgeStackNode {
+	TESShalfEdge *edge;
+	struct EdgeStackNode *next;
+};
+
+struct EdgeStack {
+	struct EdgeStackNode *top;
+};
+
+void stackInit(struct EdgeStack *stack)
+{
+    stack->top = NULL;
+}
+
+int stackEmpty(struct EdgeStack *stack)
+{
+	return stack->top == NULL;
+}
+
+void stackPush(struct EdgeStack *stack, TESShalfEdge *e)
+{
+	struct EdgeStackNode *node = malloc(sizeof(struct EdgeStackNode));
+	node->edge = e;
+	node->next = stack->top;
+	stack->top = node;
+}
+
+TESShalfEdge *stackPop(struct EdgeStack *stack)
+{
+    TESShalfEdge *e = NULL;
+	struct EdgeStackNode *node = stack->top;
+	if (node) {
+		stack->top = node->next;
+		e = node->edge;
+            free(node);
+	}
+	return e;
+}
+
+/*
+	Starting with a valid triangulation, uses the Edge Flip algorithm to
+	refine the triangulation into a Constrained Delaunay Triangulation.
+*/
+int tessMeshRefineDelaunay( TESSmesh *mesh )
+{
+	/* At this point, we have a valid, but not optimal, triangulation.
+		 We refine the triangulation using the Edge Flip algorithm */
+
+/*
+	 1) Find all internal edges
+	 2) Mark all dual edges
+	 3) insert all dual edges into a queue
+*/
+	TESSface *f;
+	struct EdgeStack stack;
+    stackInit(&stack);
+	TESShalfEdge *e;
+	TESShalfEdge *edges[4];
+	for( f = mesh->fHead.next; f != &mesh->fHead; f = f->next ) {
+		if ( f->inside) {
+			e = f->anEdge;
+			do {
+				e->mark = EdgeIsInternal(e); /* Mark internal edges */
+				if (e->mark && !e->Sym->mark) stackPush(&stack, e); /* Insert into queue */
+				e = e->Lnext;
+			} while (e != f->anEdge);
+		}
+	}
+
+	// Pop stack until we find a reversed edge
+	// Flip the reversed edge, and insert any of the four opposite edges
+	// which are internal and not already in the stack (!marked)
+	while (!stackEmpty(&stack)) {
+		e = stackPop(&stack);
+		e->mark = e->Sym->mark = 0;
+		if (!tesedgeIsLocallyDelaunay(e)) {
+			tessMeshFlipEdge(mesh, e);
+			// for each opposite edge
+			edges[0] = e->Lnext;
+			edges[1] = e->Lprev;
+			edges[2] = e->Sym->Lnext;
+			edges[3] = e->Sym->Lprev;
+			for (int i=0;i<3;i++) {
+				if (!edges[i]->mark && EdgeIsInternal(edges[i])) {
+					edges[i]->mark = edges[i]->Sym->mark = 1;
+					stackPush(&stack, edges[i]);
+				}
+			}
+		}
+	}
 
 	return 1;
 }
@@ -926,6 +1019,11 @@ int tessTesselate( TESStesselator *tess, int windingRule, int elementType,
 		rc = tessMeshSetWindingNumber( mesh, 1, TRUE );
 	} else {
 		rc = tessMeshTessellateInterior( mesh ); 
+		if (elementType == TESS_CONSTRAINED_DELAUNAY_TRIANGLES) {
+			rc = tessMeshRefineDelaunay( mesh );
+			elementType = TESS_POLYGONS;
+			polySize = 3;
+		}
 	}
 	if (rc == 0) longjmp(tess->env,1);  /* could've used a label */
 

alg_outline.md

@@ -223,3 +223,11 @@ triangles into fans and strips.  We do this using a greedy approach.
 The triangulation itself is not optimized to reduce the number of
 primitives; we just try to get a reasonable decomposition of the
 computed triangulation.
+
+Optionally, it's possible to output a Constrained Delaunay Triangulation.
+This is done by doing a delaunay refinement with the normal triangulation as
+a basis. The Edge Flip algorithm is used, which is guaranteed to terminate in O(n^2).
+
+Note: We don't use robust predicates to check if edges are locally
+delaunay, but currently us a naive epsilon of 0.01 radians to ensure
+termination.