Initial implementation of Constrained Delaunay Refinement
diff --git a/Include/tesselator.h b/Include/tesselator.h
index c27541e..1cc4166 100755
--- a/Include/tesselator.h
+++ b/Include/tesselator.h
@@ -107,11 +107,16 @@
// 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 @@
// 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:
diff --git a/Source/geom.c b/Source/geom.c
index a49ea17..4c49d58 100755
--- a/Source/geom.c
+++ b/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 @@
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);
+}
diff --git a/Source/geom.h b/Source/geom.h
index cf83897..c29a932 100755
--- a/Source/geom.h
+++ b/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 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
diff --git a/Source/mesh.c b/Source/mesh.c
index 06a8ced..1b0ff9e 100755
--- a/Source/mesh.c
+++ b/Source/mesh.c
@@ -80,6 +80,7 @@
e->Lface = NULL;
e->winding = 0;
e->activeRegion = NULL;
+ e->mark = 0;
eSym->Sym = e;
eSym->Onext = eSym;
@@ -88,6 +89,7 @@
eSym->Lface = NULL;
eSym->winding = 0;
eSym->activeRegion = NULL;
+ eSym->mark = 0;
return e;
}
@@ -748,6 +750,85 @@
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
diff --git a/Source/mesh.h b/Source/mesh.h
index c492916..479c66f 100755
--- a/Source/mesh.h
+++ b/Source/mesh.h
@@ -143,6 +143,7 @@
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 @@
#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 @@
void tessMeshDeleteMesh( TESSalloc* alloc, TESSmesh *mesh );
void tessMeshZapFace( TESSmesh *mesh, TESSface *fZap );
+void tessMeshFlipEdge( TESSmesh *mesh, TESShalfEdge *edge );
+
#ifdef NDEBUG
#define tessMeshCheckMesh( mesh )
#else
diff --git a/Source/tess.c b/Source/tess.c
index 013b84c..e1d3b63 100755
--- a/Source/tess.c
+++ b/Source/tess.c
@@ -365,7 +365,6 @@
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 @@
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 @@
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 */
diff --git a/alg_outline.md b/alg_outline.md
index 449ca2a..9a4a2bd 100644
--- a/alg_outline.md
+++ b/alg_outline.md
@@ -223,3 +223,11 @@
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.
