Merge branch 'pr/7'
diff --git a/Example/example.c b/Example/example.c
index acbc388..2d4e23b 100644
--- a/Example/example.c
+++ b/Example/example.c
@@ -56,6 +56,7 @@
 
 
 int run = 1;
+int cdt = 0;
 
 static void key(GLFWwindow* window, int key, int scancode, int action, int mods)
 {
@@ -65,6 +66,8 @@
 		glfwSetWindowShouldClose(window, GL_TRUE);
 	if (key == GLFW_KEY_SPACE && action == GLFW_PRESS)
 		run = !run;
+	if (key == GLFW_KEY_C && action == GLFW_PRESS)
+		cdt = !cdt;
 }
 
 int main(int argc, char *argv[])
@@ -79,14 +82,15 @@
 	float t = 0.0f, pt = 0.0f;
 	TESSalloc ma;
 	TESStesselator* tess = 0;
-	const int nvp = 6;
+	const int nvp = 3;
 	unsigned char* vflags = 0;
-	int nvflags = 0;
 #ifdef USE_POOL
 	struct MemPool pool;
 	unsigned char mem[1024*1024];
+	int nvflags = 0;
 #else
 	int allocated = 0;
+	double t0 = 0, t1 = 0;
 #endif
 	TESS_NOTUSED(argc);
 	TESS_NOTUSED(argv);
@@ -104,7 +108,7 @@
 	if (!fg) return -1;
 
 	printf("go...\n");
-	
+
 	// Flip y
 	for (it = bg; it != NULL; it = it->next)
 		for (i = 0; i < it->npts; ++i)
@@ -138,7 +142,7 @@
 			it->pts[i*2+1] -= cy;
 		}
 	}
-			
+
 	// Find BG bounds.
 	bounds[0] = bounds[2] = bg->pts[0];
 	bounds[1] = bounds[3] = bg->pts[1];
@@ -154,7 +158,7 @@
 			if (y > bounds[3]) bounds[3] = y;
 		}
 	}
-		
+
 #ifdef USE_POOL
 
 	pool.size = 0;
@@ -168,6 +172,8 @@
 
 #else
 
+	t0 = glfwGetTime();
+
 	memset(&ma, 0, sizeof(ma));
 	ma.memalloc = stdAlloc;
 	ma.memfree = stdFree;
@@ -178,6 +184,8 @@
 	if (!tess)
 		return -1;
 
+	tessSetOption(tess, TESS_CONSTRAINED_DELAUNAY_TRIANGULATION, 1);
+
 	// Offset the foreground shape to center of the bg.
 	offx = (bounds[2]+bounds[0])/2;
 	offy = (bounds[3]+bounds[1])/2;
@@ -189,7 +197,7 @@
 			it->pts[i*2+1] += offy;
 		}
 	}
-	
+
 	// Add contours.
 	for (it = bg; it != NULL; it = it->next)
 		tessAddContour(tess, 2, it->pts, sizeof(float)*2, it->npts);
@@ -197,10 +205,14 @@
 		tessAddContour(tess, 2, it->pts, sizeof(float)*2, it->npts);
 	if (!tessTesselate(tess, TESS_WINDING_POSITIVE, TESS_POLYGONS, nvp, 2, 0))
 		return -1;
+
+	t1 = glfwGetTime();
+
+	printf("Time: %.3f ms\n", (t1 - t0) * 1000.0f);
 	printf("Memory used: %.1f kB\n", allocated/1024.0f);
-	
+
 #endif
-	
+
 	mode = glfwGetVideoMode(glfwGetPrimaryMonitor());
 	width = mode->width - 40;
 	height = mode->height - 80;
@@ -221,17 +233,27 @@
 	view[2] = cx + w*1.2f;
 	view[1] = cy - w*1.2f*(float)height/(float)width;
 	view[3] = cy + w*1.2f*(float)height/(float)width;
-		
+
 	glfwSetTime(0);
 
 	while (!glfwWindowShouldClose(window))
 	{
-		float ct = (float)glfwGetTime();
+		int winWidth, winHeight;
+		int fbWidth, fbHeight;
+		float pxr, ct;
+
+		glfwGetWindowSize(window, &winWidth, &winHeight);
+		glfwGetFramebufferSize(window, &fbWidth, &fbHeight);
+
+		// Calculate pixel ration for hi-dpi devices.
+		pxr = (float)fbWidth / (float)winWidth;
+
+		ct = (float)glfwGetTime();
 		if (run) t += ct - pt;
 		pt = ct;
-		
+
 		// Update and render
-		glViewport(0, 0, width, height);
+		glViewport(0, 0, fbWidth, fbHeight);
 		glClearColor(0.3f, 0.3f, 0.32f, 1.0f);
 		glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT);
 		glEnable(GL_BLEND);
@@ -251,9 +273,11 @@
 		tess = tessNewTess(&ma);
 		if (tess)
 		{
+			tessSetOption(tess, TESS_CONSTRAINED_DELAUNAY_TRIANGULATION, cdt);
+
 			offx = (view[2]+view[0])/2 + sinf(t) * (view[2]-view[0])/2;
 			offy = (view[3]+view[1])/2 + cosf(t*3.13f) * (view[3]-view[1])/6;
-			
+
 			for (it = fg; it != NULL; it = it->next)
 			{
 				for (i = 0; i < it->npts; ++i)
@@ -293,7 +317,7 @@
 					nvflags = nverts;
 					vflags = (unsigned char*)malloc(sizeof(unsigned char)*nvflags);
 				}
-				
+
 				if (vflags)
 				{
 					// Vertex indices describe the order the indices were added and can be used
@@ -303,7 +327,7 @@
 					for (i = 0; i < nverts; ++i)
 						vflags[i] = vinds[i] == TESS_UNDEF ? 1 : 0;
 				}
-				
+
 				for (i = 0; i < nelems; ++i)
 				{
 					int b = elems[i*2];
@@ -314,9 +338,9 @@
 					tess = 0;
 			}
 			else
-				tess = 0;				
+				tess = 0;
 		}
-#endif		
+#endif
 
 		// Draw tesselated pieces.
 		if (tess)
@@ -326,7 +350,7 @@
 			const int* elems = tessGetElements(tess);
 			const int nverts = tessGetVertexCount(tess);
 			const int nelems = tessGetElementCount(tess);
-			
+
 			// Draw polygons.
 			glColor4ub(255,255,255,128);
 			for (i = 0; i < nelems; ++i)
@@ -337,7 +361,10 @@
 					glVertex2f(verts[p[j]*2], verts[p[j]*2+1]);
 				glEnd();
 			}
-			
+
+			glLineWidth(1.0f * pxr);
+			glPointSize(3.0f * pxr);
+
 			glColor4ub(0,0,0,16);
 			for (i = 0; i < nelems; ++i)
 			{
@@ -347,9 +374,8 @@
 					glVertex2f(verts[p[j]*2], verts[p[j]*2+1]);
 				glEnd();
 			}
-			
+
 			glColor4ub(0,0,0,128);
-			glPointSize(3.0f);
 			glBegin(GL_POINTS);
 			for (i = 0; i < nverts; ++i)
 			{
@@ -360,21 +386,22 @@
 				glVertex2f(verts[i*2], verts[i*2+1]);
 			}
 			glEnd();
+
 			glPointSize(1.0f);
 		}
-		
+
 		glEnable(GL_DEPTH_TEST);
 		glfwSwapBuffers(window);
 		glfwPollEvents();
 	}
-	
+
 	if (tess) tessDeleteTess(tess);
-	
+
 	if (vflags)
 		free(vflags);
-	
-	svgDelete(bg);	
-	svgDelete(fg);	
+
+	svgDelete(bg);
+	svgDelete(fg);
 
 	glfwTerminate();
 	return 0;
diff --git a/Include/tesselator.h b/Include/tesselator.h
index c27541e..38e05f1 100755
--- a/Include/tesselator.h
+++ b/Include/tesselator.h
@@ -1,5 +1,5 @@
 /*
-** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008) 
+** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
 ** Copyright (C) [dates of first publication] Silicon Graphics, Inc.
 ** All Rights Reserved.
 **
@@ -9,10 +9,10 @@
 ** to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
 ** of the Software, and to permit persons to whom the Software is furnished to do so,
 ** subject to the following conditions:
-** 
+**
 ** The above copyright notice including the dates of first publication and either this
 ** permission notice or a reference to http://oss.sgi.com/projects/FreeB/ shall be
-** included in all copies or substantial portions of the Software. 
+** included in all copies or substantial portions of the Software.
 **
 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
 ** INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
@@ -20,7 +20,7 @@
 ** BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
 ** OR OTHER DEALINGS IN THE SOFTWARE.
-** 
+**
 ** Except as contained in this notice, the name of Silicon Graphics, Inc. shall not
 ** be used in advertising or otherwise to promote the sale, use or other dealings in
 ** this Software without prior written authorization from Silicon Graphics, Inc.
@@ -106,7 +106,7 @@
 //         }
 //         glEnd();
 //     }
-//
+
 enum TessElementType
 {
 	TESS_POLYGONS,
@@ -114,6 +114,16 @@
 	TESS_BOUNDARY_CONTOURS,
 };
 
+
+// TESS_CONSTRAINED_DELAUNAY_TRIANGULATION
+//   If enabled, the initial triagulation is improved with non-robust Constrained Delayney triangulation.
+//   Disable by default.
+
+enum TessOption
+{
+	TESS_CONSTRAINED_DELAUNAY_TRIANGULATION,
+};
+
 typedef float TESSreal;
 typedef int TESSindex;
 typedef struct TESStesselator TESStesselator;
@@ -133,12 +143,12 @@
 // how often to allocate memory from the system versus how much extra space the system
 // should allocate. Reasonable defaults are show in commects below, they will be used if
 // the bucket sizes are zero.
-// 
+//
 // The use may left the memrealloc to be null. In that case, the tesselator will not try to
 // dynamically grow int's internal arrays. The tesselator only needs the reallocation when it
 // has found intersecting segments and needs to add new vertex. This defency can be cured by
 // allocating some extra vertices beforehand. The 'extraVertices' variable allows to specify
-// number of expected extra vertices.  
+// number of expected extra vertices.
 struct TESSalloc
 {
 	void *(*memalloc)( void *userData, unsigned int size );
@@ -184,6 +194,12 @@
 //   count - number of vertices in contour.
 void tessAddContour( TESStesselator *tess, int size, const void* pointer, int stride, int count );
 
+// tessSetOption() - Toggles optional tessellation parameters
+// Parameters:
+//  option - one of TessOption
+//  value - 1 if enabled, 0 if disabled.
+void tessSetOption( TESStesselator *tess, int option, int value );
+
 // tessTesselate() - tesselate contours.
 // Parameters:
 //   tess - pointer to tesselator object.
@@ -207,7 +223,7 @@
 // Every point added using tessAddContour() will get a new index starting at 0.
 // New vertices generated at the intersections of segments are assigned value TESS_UNDEF.
 const TESSindex* tessGetVertexIndices( TESStesselator *tess );
-	
+
 // tessGetElementCount() - Returns number of elements in the the tesselated output.
 int tessGetElementCount( TESStesselator *tess );
 
diff --git a/Source/geom.c b/Source/geom.c
index a49ea17..97a5046 100755
--- a/Source/geom.c
+++ b/Source/geom.c
@@ -1,5 +1,5 @@
 /*

-** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008) 

+** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)

 ** Copyright (C) [dates of first publication] Silicon Graphics, Inc.

 ** All Rights Reserved.

 **

@@ -9,10 +9,10 @@
 ** to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies

 ** of the Software, and to permit persons to whom the Software is furnished to do so,

 ** subject to the following conditions:

-** 

+**

 ** The above copyright notice including the dates of first publication and either this

 ** permission notice or a reference to http://oss.sgi.com/projects/FreeB/ shall be

-** included in all copies or substantial portions of the Software. 

+** included in all copies or substantial portions of the Software.

 **

 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,

 ** INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A

@@ -20,7 +20,7 @@
 ** BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,

 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE

 ** OR OTHER DEALINGS IN THE SOFTWARE.

-** 

+**

 ** Except as contained in this notice, the name of Silicon Graphics, Inc. shall not

 ** be used in advertising or otherwise to promote the sale, use or other dealings in

 ** this Software without prior written authorization from Silicon Graphics, Inc.

@@ -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,30 @@
 		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 num, den, ax, ay, bx, by;

+	ax = v2->s - v1->s;

+	ay = v2->t - v1->t;

+	bx = v0->s - v1->s;

+	by = v0->t - v1->t;

+	num = ax * bx + ay * by;

+	den = sqrt( ax * ax + ay * ay ) * sqrt( bx * bx + by * by );

+	if ( den > 0.0 ) num /= den;

+	if ( num < -1.0 ) num = -1.0;

+	if ( num >  1.0 ) num =  1.0;

+	return acos( num );

+}

+

+/*

+	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..15ce620 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,83 @@
 	return 1;
 }
 
+void tessMeshFlipEdge( TESSmesh *mesh, TESShalfEdge *edge )
+{
+	TESShalfEdge *a0 = edge;
+	TESShalfEdge *a1 = a0->Lnext;
+	TESShalfEdge *a2 = a1->Lnext;
+	TESShalfEdge *b0 = edge->Sym;
+	TESShalfEdge *b1 = b0->Lnext;
+	TESShalfEdge *b2 = b1->Lnext;
+
+	TESSvertex *aOrg = a0->Org;
+	TESSvertex *aOpp = a2->Org;
+	TESSvertex *bOrg = b0->Org;
+	TESSvertex *bOpp = b2->Org;
+
+	TESSface *fa = a0->Lface;
+	TESSface *fb = b0->Lface;
+
+	assert(EdgeIsInternal(edge));
+	assert(a2->Lnext == a0);
+	assert(b2->Lnext == b0);
+
+	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 67c038f..8059bf4 100755
--- a/Source/tess.c
+++ b/Source/tess.c
@@ -1,5 +1,5 @@
 /*
-** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008) 
+** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
 ** Copyright (C) [dates of first publication] Silicon Graphics, Inc.
 ** All Rights Reserved.
 **
@@ -9,10 +9,10 @@
 ** to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
 ** of the Software, and to permit persons to whom the Software is furnished to do so,
 ** subject to the following conditions:
-** 
+**
 ** The above copyright notice including the dates of first publication and either this
 ** permission notice or a reference to http://oss.sgi.com/projects/FreeB/ shall be
-** included in all copies or substantial portions of the Software. 
+** included in all copies or substantial portions of the Software.
 **
 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
 ** INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
@@ -20,7 +20,7 @@
 ** BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
 ** OR OTHER DEALINGS IN THE SOFTWARE.
-** 
+**
 ** Except as contained in this notice, the name of Silicon Graphics, Inc. shall not
 ** be used in advertising or otherwise to promote the sale, use or other dealings in
 ** this Software without prior written authorization from Silicon Graphics, Inc.
@@ -186,7 +186,7 @@
 #define S_UNIT_X	(RandomSweep ? (2*drand48()-1) : 1.0)
 #define S_UNIT_Y	(RandomSweep ? (2*drand48()-1) : 0.0)
 #else
-#if defined(SLANTED_SWEEP) 
+#if defined(SLANTED_SWEEP)
 /* The "feature merging" is not intended to be complete.  There are
 * special cases where edges are nearly parallel to the sweep line
 * which are not implemented.  The algorithm should still behave
@@ -295,7 +295,7 @@
 * (what else would it do??)  The region must consist of a single
 * loop of half-edges (see mesh.h) oriented CCW.  "Monotone" in this
 * case means that any vertical line intersects the interior of the
-* region in a single interval.  
+* region in a single interval.
 *
 * Tessellation consists of adding interior edges (actually pairs of
 * half-edges), to split the region into non-overlapping triangles.
@@ -374,7 +374,6 @@
 	return 1;
 }
 
-
 /* tessMeshTessellateInterior( mesh ) tessellates each region of
 * the mesh which is marked "inside" the polygon.  Each such region
 * must be monotone.
@@ -391,11 +390,125 @@
 			if ( !tessMeshTessellateMonoRegion( mesh, f ) ) return 0;
 		}
 	}
-
 	return 1;
 }
 
 
+typedef struct EdgeStackNode EdgeStackNode;
+typedef struct EdgeStack EdgeStack;
+
+struct EdgeStackNode {
+	TESShalfEdge *edge;
+	EdgeStackNode *next;
+};
+
+struct EdgeStack {
+	EdgeStackNode *top;
+	struct BucketAlloc *nodeBucket;
+};
+
+int stackInit( EdgeStack *stack, TESSalloc *alloc )
+{
+	stack->top = NULL;
+	stack->nodeBucket = createBucketAlloc( alloc, "CDT nodes", sizeof(EdgeStackNode), 512 );
+	return stack->nodeBucket != NULL;
+}
+
+void stackDelete( EdgeStack *stack )
+{
+    deleteBucketAlloc( stack->nodeBucket );
+}
+
+int stackEmpty( EdgeStack *stack )
+{
+	return stack->top == NULL;
+}
+
+void stackPush( EdgeStack *stack, TESShalfEdge *e )
+{
+	EdgeStackNode *node = (EdgeStackNode *)bucketAlloc( stack->nodeBucket );
+	if ( ! node ) return;
+	node->edge = e;
+	node->next = stack->top;
+	stack->top = node;
+}
+
+TESShalfEdge *stackPop( EdgeStack *stack )
+{
+	TESShalfEdge *e = NULL;
+	EdgeStackNode *node = stack->top;
+	if (node) {
+		stack->top = node->next;
+		e = node->edge;
+		bucketFree( stack->nodeBucket, node );
+	}
+	return e;
+}
+
+
+//	Starting with a valid triangulation, uses the Edge Flip algorithm to
+//	refine the triangulation into a Constrained Delaunay Triangulation.
+void tessMeshRefineDelaunay( TESSmesh *mesh, TESSalloc *alloc )
+{
+	// 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;
+	EdgeStack stack;
+	TESShalfEdge *e;
+	int maxFaces = 0, maxIter = 0, iter = 0;
+
+	stackInit(&stack, alloc);
+
+	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);
+			maxFaces++;
+		}
+	}
+
+	// The algorithm should converge on O(n^2), since the predicate is not robust,
+	// we'll save guard against infinite loop.
+	maxIter = maxFaces * maxFaces;
+
+	// 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) && iter < maxIter) {
+		e = stackPop(&stack);
+		e->mark = e->Sym->mark = 0;
+		if (!tesedgeIsLocallyDelaunay(e)) {
+			TESShalfEdge *edges[4];
+			int i;
+			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 (i = 0; i < 4; i++) {
+				if (!edges[i]->mark && EdgeIsInternal(edges[i])) {
+					edges[i]->mark = edges[i]->Sym->mark = 1;
+					stackPush(&stack, edges[i]);
+				}
+			}
+		}
+		iter++;
+	}
+
+	stackDelete(&stack);
+}
+
+
 /* tessMeshDiscardExterior( mesh ) zaps (ie. sets to NULL) all faces
 * which are not marked "inside" the polygon.  Since further mesh operations
 * on NULL faces are not allowed, the main purpose is to clean up the
@@ -485,7 +598,7 @@
 
 	if (alloc == NULL)
 		alloc = &defaulAlloc;
-	
+
 	/* Only initialize fields which can be changed by the api.  Other fields
 	* are initialized where they are used.
 	*/
@@ -506,7 +619,7 @@
 		tess->alloc.dictNodeBucketSize = 512;
 	if (tess->alloc.regionBucketSize == 0)
 		tess->alloc.regionBucketSize = 256;
-	
+
 	tess->normal[0] = 0;
 	tess->normal[1] = 0;
 	tess->normal[2] = 0;
@@ -530,7 +643,7 @@
 
 	tess->outOfMemory = 0;
 	tess->vertexIndexCounter = 0;
-	
+
 	tess->vertices = 0;
 	tess->vertexIndices = 0;
 	tess->vertexCount = 0;
@@ -542,9 +655,9 @@
 
 void tessDeleteTess( TESStesselator *tess )
 {
-	
+
 	struct TESSalloc alloc = tess->alloc;
-	
+
 	deleteBucketAlloc( tess->regionPool );
 
 	if( tess->mesh != NULL ) {
@@ -623,7 +736,7 @@
 			edge = edge->Lnext;
 		}
 		while (edge != f->anEdge);
-		
+
 		assert( faceVerts <= polySize );
 
 		f->n = maxFaceCount;
@@ -640,7 +753,7 @@
 		tess->outOfMemory = 1;
 		return;
 	}
-	
+
 	tess->vertexCount = maxVertexCount;
 	tess->vertices = (TESSreal*)tess->alloc.memalloc( tess->alloc.userData,
 													 sizeof(TESSreal) * tess->vertexCount * vertexSize );
@@ -657,7 +770,7 @@
 		tess->outOfMemory = 1;
 		return;
 	}
-	
+
 	// Output vertices.
 	for ( v = mesh->vHead.next; v != &mesh->vHead; v = v->next )
 	{
@@ -679,7 +792,7 @@
 	for ( f = mesh->fHead.next; f != &mesh->fHead; f = f->next )
 	{
 		if ( !f->inside ) continue;
-		
+
 		// Store polygon
 		edge = f->anEdge;
 		faceVerts = 0;
@@ -748,7 +861,7 @@
 		tess->outOfMemory = 1;
 		return;
 	}
-	
+
 	tess->vertices = (TESSreal*)tess->alloc.memalloc( tess->alloc.userData,
 													  sizeof(TESSreal) * tess->vertexCount * vertexSize );
 	if (!tess->vertices)
@@ -764,7 +877,7 @@
 		tess->outOfMemory = 1;
 		return;
 	}
-	
+
 	verts = tess->vertices;
 	elements = tess->elements;
 	vertInds = tess->vertexIndices;
@@ -865,6 +978,17 @@
 	}
 }
 
+void tessSetOption( TESStesselator *tess, int option, int value )
+{
+	switch(option)
+	{
+	case TESS_CONSTRAINED_DELAUNAY_TRIANGULATION:
+		tess->processCDT = value > 0 ? 1 : 0;
+		break;
+	}
+}
+
+
 int tessTesselate( TESStesselator *tess, int windingRule, int elementType,
 				  int polySize, int vertexSize, const TESSreal* normal )
 {
@@ -885,7 +1009,7 @@
 	}
 
 	tess->vertexIndexCounter = 0;
-	
+
 	if (normal)
 	{
 		tess->normal[0] = normal[0];
@@ -900,7 +1024,7 @@
 	if (vertexSize > 3)
 		vertexSize = 3;
 
-	if (setjmp(tess->env) != 0) { 
+	if (setjmp(tess->env) != 0) {
 		/* come back here if out of memory */
 		return 0;
 	}
@@ -934,7 +1058,9 @@
 	if (elementType == TESS_BOUNDARY_CONTOURS) {
 		rc = tessMeshSetWindingNumber( mesh, 1, TRUE );
 	} else {
-		rc = tessMeshTessellateInterior( mesh ); 
+		rc = tessMeshTessellateInterior( mesh );
+		if (rc != 0 && tess->processCDT != 0)
+			tessMeshRefineDelaunay( mesh, &tess->alloc );
 	}
 	if (rc == 0) longjmp(tess->env,1);  /* could've used a label */
 
diff --git a/Source/tess.h b/Source/tess.h
index c8c747b..c1f981f 100755
--- a/Source/tess.h
+++ b/Source/tess.h
@@ -1,5 +1,5 @@
 /*
-** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008) 
+** SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
 ** Copyright (C) [dates of first publication] Silicon Graphics, Inc.
 ** All Rights Reserved.
 **
@@ -9,10 +9,10 @@
 ** to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
 ** of the Software, and to permit persons to whom the Software is furnished to do so,
 ** subject to the following conditions:
-** 
+**
 ** The above copyright notice including the dates of first publication and either this
 ** permission notice or a reference to http://oss.sgi.com/projects/FreeB/ shall be
-** included in all copies or substantial portions of the Software. 
+** included in all copies or substantial portions of the Software.
 **
 ** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
 ** INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
@@ -20,7 +20,7 @@
 ** BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
 ** OR OTHER DEALINGS IN THE SOFTWARE.
-** 
+**
 ** Except as contained in this notice, the name of Silicon Graphics, Inc. shall not
 ** be used in advertising or otherwise to promote the sale, use or other dealings in
 ** this Software without prior written authorization from Silicon Graphics, Inc.
@@ -61,6 +61,8 @@
 	TESSreal bmin[2];
 	TESSreal bmax[2];
 
+	int processCDT;	/* option to run Constrained Delayney pass. */
+
 	/*** state needed for the line sweep ***/
 	int	windingRule;	/* rule for determining polygon interior */
 
@@ -71,7 +73,7 @@
 	struct BucketAlloc* regionPool;
 
 	TESSindex vertexIndexCounter;
-	
+
 	TESSreal *vertices;
 	TESSindex *vertexIndices;
 	int vertexCount;
@@ -79,7 +81,7 @@
 	int elementCount;
 
 	TESSalloc alloc;
-	
+
 	jmp_buf env;			/* place to jump to when memAllocs fail */
 };
 
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.