diff --git a/Tools/Construct/Clipper/Makefile.am b/Tools/Construct/Clipper/Makefile.am
index abeabbf9f..8d0745b0c 100644
--- a/Tools/Construct/Clipper/Makefile.am
+++ b/Tools/Construct/Clipper/Makefile.am
@@ -6,11 +6,16 @@ noinst_PROGRAMS = testclipper
 
 testclipper_SOURCES = testclipper.cxx
 
-testclipper_LDADD = $(top_builddir)/Tools/Construct/Clipper/libClipper.a \
+testclipper_LDADD = \
+	$(top_builddir)/Tools/Construct/Clipper/libClipper.a \
+	$(top_builddir)/Tools/Construct/Triangulate/libTriangulate.a \
 	$(top_builddir)/Tools/Lib/Polygon/libPolygon.a \
 	$(top_builddir)/Lib/Debug/libDebug.a \
 	$(top_builddir)/Lib/Misc/libMisc.a \
 	$(top_builddir)/Lib/zlib/libz.a \
 	-lgfc -lgpc
 
-INCLUDES += -I$(top_builddir) -I$(top_builddir)/Lib -I$(top_builddir)/Tools/Lib
+INCLUDES += -I$(top_builddir) \
+	-I$(top_builddir)/Lib \
+	-I$(top_builddir)/Tools/Lib \
+	-I$(top_builddir)/Tools/Construct
diff --git a/Tools/Construct/Clipper/clipper.cxx b/Tools/Construct/Clipper/clipper.cxx
index 78fc11d08..f3c69cb8d 100644
--- a/Tools/Construct/Clipper/clipper.cxx
+++ b/Tools/Construct/Clipper/clipper.cxx
@@ -44,8 +44,8 @@ FGClipper::~FGClipper( void ) {
 
 // Initialize Clipper (allocate and/or connect structures)
 bool FGClipper::init() {
-    v_list.num_vertices = 0;
-    v_list.vertex = new gpc_vertex[FG_MAX_VERTICES];;
+    // v_list.num_vertices = 0;
+    // v_list.vertex = new gpc_vertex[FG_MAX_VERTICES];;
 
     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 	polys_in.polys[i].clear();
@@ -72,10 +72,12 @@ bool FGClipper::load_polys(const string& path) {
 	exit(-1);
     }
 
-    gpc_polygon *poly = new gpc_polygon;
-    poly->num_contours = 0;
-    poly->contour = NULL;
+    // gpc_polygon *poly = new gpc_polygon;
+    // poly->num_contours = 0;
+    // poly->contour = NULL;
+    FGPolygon poly;
 
+    Point3D p;
     in >> skipcomment;
     while ( !in.eof() ) {
 	in >> poly_name;
@@ -85,6 +87,8 @@ bool FGClipper::load_polys(const string& path) {
 	in >> contours;
 	cout << "num contours = " << contours << endl;
 
+	poly.erase();
+
 	for ( i = 0; i < contours; ++i ) {
 	    in >> count;
 
@@ -98,19 +102,18 @@ bool FGClipper::load_polys(const string& path) {
 
 	    in >> startx;
 	    in >> starty;
-	    v_list.vertex[0].x = startx;
-	    v_list.vertex[0].y = starty;
+	    p = Point3D(startx, starty, 0.0);
+	    poly.add_node( i, p );
 	    FG_LOG( FG_CLIPPER, FG_BULK, "0 = " 
 		    << startx << ", " << starty );
 
 	    for ( j = 1; j < count - 1; ++j ) {
 		in >> x;
 		in >> y;
-		v_list.vertex[j].x = x;
-		v_list.vertex[j].y = y;
+		p = Point3D( x, y, 0.0 );
+		poly.add_node( i, p );
 		FG_LOG( FG_CLIPPER, FG_BULK, j << " = " << x << ", " << y );
 	    }
-	    v_list.num_vertices = count - 1;
 
 	    in >> lastx;
 	    in >> lasty;
@@ -119,14 +122,13 @@ bool FGClipper::load_polys(const string& path) {
 		 && (fabs(starty - lasty) < FG_EPSILON) ) {
 		// last point same as first, discard
 	    } else {
-		v_list.vertex[count - 1].x = lastx;
-		v_list.vertex[count - 1].y = lasty;
-		++v_list.num_vertices;
+		p = Point3D( lastx, lasty, 0.0 );
+		poly.add_node( i, p );
 		FG_LOG( FG_CLIPPER, FG_BULK, count - 1 << " = " 
 			<< lastx << ", " << lasty );
 	    }
 
-	    gpc_add_contour( poly, &v_list, hole_flag );
+	    // gpc_add_contour( poly, &v_list, hole_flag );
 	}
 
 	in >> skipcomment;
@@ -152,15 +154,15 @@ bool FGClipper::load_polys(const string& path) {
 
 // merge any slivers in the specified polygon with larger
 // neighboring polygons of higher priorigy
-void FGClipper::merge_slivers(gpc_polygon *poly) {
+void FGClipper::merge_slivers(FGPolygon& poly) {
     cout << "Begin merge slivers" << endl;
     // traverse each contour of the polygon and attempt to identify
     // likely slivers
-    for ( int i = 0; i < poly->num_contours; i++ ) {
+    for ( int i = 0; i < poly.contours(); ++i ) {
 	cout << "contour " << i << endl;
-	for (int j = 0; j < poly->contour[i].num_vertices; j++ ) {
-	    cout << poly->contour[i].vertex[j].x << ","
-		 << poly->contour[i].vertex[j].y << endl;
+	for (int j = 0; j < poly.contour_size( i ); ++j ) {
+	    // cout << poly->contour[i].vertex[j].x << ","
+	    //      << poly->contour[i].vertex[j].y << endl;
 	}
     }
 }
@@ -168,62 +170,53 @@ void FGClipper::merge_slivers(gpc_polygon *poly) {
 
 // Do actually clipping work
 bool FGClipper::clip_all(const point2d& min, const point2d& max) {
-    gpc_polygon accum, result_union, tmp;
-    gpc_polygon *result_diff, *remains;
-    gpcpoly_iterator current, last;
+    FGPolygon accum, result_union, tmp;
+    FGPolygon result_diff, remains;
+    // gpcpoly_iterator current, last;
 
     FG_LOG( FG_CLIPPER, FG_INFO, "Running master clipper" );
 
-    accum.num_contours = 0;
+    accum.erase();
 
     cout << "  (" << min.x << "," << min.y << ") (" 
 	 << max.x << "," << max.y << ")" << endl;
 
     // set up clipping tile
-    v_list.vertex[0].x = min.x;
-    v_list.vertex[0].y = min.y;
-
-    v_list.vertex[1].x = max.x;
-    v_list.vertex[1].y = min.y;
-
-    v_list.vertex[2].x = max.x;
-    v_list.vertex[2].y = max.y;
-
-    v_list.vertex[3].x = min.x;
-    v_list.vertex[3].y = max.y;
-
-    v_list.num_vertices = 4;
-
-    polys_in.safety_base.num_contours = 0;
-    polys_in.safety_base.contour = NULL;
-    gpc_add_contour( &polys_in.safety_base, &v_list, 0 );
+    polys_in.safety_base.erase();
+    polys_in.safety_base.add_node( 0, Point3D(min.x, min.y, 0.0) );
+    polys_in.safety_base.add_node( 0, Point3D(max.x, min.y, 0.0) );
+    polys_in.safety_base.add_node( 0, Point3D(max.x, max.y, 0.0) );
+    polys_in.safety_base.add_node( 0, Point3D(min.x, max.y, 0.0) );
 
     // int count = 0;
     // process polygons in priority order
     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 	cout << "num polys of type (" << i << ") = " 
 	     << polys_in.polys[i].size() << endl;
-	current = polys_in.polys[i].begin();
-	last = polys_in.polys[i].end();
-	for ( ; current != last; ++current ) {
+	// current = polys_in.polys[i].begin();
+	// last = polys_in.polys[i].end();
+	// for ( ; current != last; ++current ) {
+	for( int j = 0; j < (int)polys_in.polys[i].size(); ++j ) {
+	    FGPolygon current = polys_in.polys[i][j];
 	    FG_LOG( FG_CLIPPER, FG_DEBUG, get_area_name( (AreaType)i ) 
-		    << " = " << (*current)->contour->num_vertices );
+		    << " = " << current.contours() );
 
 #ifdef EXTRA_SAFETY_CLIP
 	    // clip to base tile
-	    gpc_polygon_clip(GPC_INT, *current, &polys_in.safety_base, &tmp);
+	    tmp = polygon_int( current, polys_in.safety_base );
 #else
-	    tmp = *current;
+	    tmp = current;
 #endif
 
 	    // clip current polygon against previous higher priority
 	    // stuff
-	    result_diff = new gpc_polygon;
-	    result_diff->num_contours = 0;
-	    result_diff->contour = NULL;
 
-	    if ( accum.num_contours == 0 ) {
-		*result_diff = tmp;
+	    // result_diff = new gpc_polygon;
+	    // result_diff->num_contours = 0;
+	    // result_diff->contour = NULL;
+
+	    if ( accum.contours() == 0 ) {
+		result_diff = tmp;
 		result_union = tmp;
 	    } else {
    		// cout << "DIFF: tmp.num_contours = " << tmp.num_contours
@@ -239,8 +232,8 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
 		// gpc_write_polygon(ofp, 1, &accum);
 		// fclose(ofp);
 
-		gpc_polygon_clip(GPC_DIFF, &tmp, &accum, result_diff);
-		gpc_polygon_clip(GPC_UNION, &tmp, &accum, &result_union);
+		result_diff = polygon_diff( tmp, accum);
+		result_union = polygon_union( tmp, accum);
 	    }
 
 	    /*
@@ -264,7 +257,7 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
 	      */
 
 	    // only add to output list if the clip left us with a polygon
-	    if ( result_diff->num_contours > 0 ) {
+	    if ( result_diff.contours() > 0 ) {
 		// merge any slivers with larger neighboring polygons
 		merge_slivers(result_diff);
 
@@ -283,15 +276,16 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
     // finally, what ever is left over goes to ocean
 
     // clip to accum against original base tile
-    remains = new gpc_polygon;
-    remains->num_contours = 0;
-    remains->contour = NULL;
-    gpc_polygon_clip(GPC_DIFF, &polys_in.safety_base, &accum, 
-		     remains);
-    if ( remains->num_contours > 0 ) {
+    // remains = new gpc_polygon;
+    // remains->num_contours = 0;
+    // remains->contour = NULL;
+    remains = polygon_diff( polys_in.safety_base, accum );
+
+    if ( remains.contours() > 0 ) {
 	polys_clipped.polys[(int)OceanArea].push_back(remains);
     }
 
+#if 0
     FILE *ofp;
 
     // tmp output accum
@@ -307,6 +301,7 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
 	gpc_write_polygon(ofp, 1, remains);
 	fclose(ofp);
     }
+#endif
 
     return true;
 }
diff --git a/Tools/Construct/Clipper/clipper.hxx b/Tools/Construct/Clipper/clipper.hxx
index e00f5789e..df0b62cfc 100644
--- a/Tools/Construct/Clipper/clipper.hxx
+++ b/Tools/Construct/Clipper/clipper.hxx
@@ -35,7 +35,7 @@
 
 #include <Include/compiler.h>
 
-
+#if 0
 // include Generic Polygon Clipping Library
 //
 //    http://www.cs.man.ac.uk/aig/staff/alan/software/
@@ -43,6 +43,9 @@
 extern "C" {
 #include <gpc.h>
 }
+#endif
+
+#include <Triangulate/polygon.hxx>
 
 #include STL_STRING
 #include <vector>
@@ -51,14 +54,14 @@ FG_USING_STD(string);
 FG_USING_STD(vector);
 
 
-typedef vector < gpc_polygon * > gpcpoly_container;
-typedef gpcpoly_container::iterator gpcpoly_iterator;
-typedef gpcpoly_container::const_iterator const_gpcpoly_iterator;
+// typedef vector < gpc_polygon * > gpcpoly_container;
+// typedef gpcpoly_container::iterator gpcpoly_iterator;
+// typedef gpcpoly_container::const_iterator const_gpcpoly_iterator;
 
 
 #define FG_MAX_AREA_TYPES 20
 #define EXTRA_SAFETY_CLIP
-#define FG_MAX_VERTICES 100000
+// #define FG_MAX_VERTICES 100000
 
 
 class point2d {
@@ -67,10 +70,10 @@ public:
 };
 
 
-class FGgpcPolyList {
+class FGPolyList {
 public:
-    gpcpoly_container polys[FG_MAX_AREA_TYPES];
-    gpc_polygon safety_base;
+    poly_list polys[FG_MAX_AREA_TYPES];
+    FGPolygon safety_base;
 };
 
 
@@ -78,9 +81,9 @@ class FGClipper {
 
 private:
 
-    gpc_vertex_list v_list;
+    // gpc_vertex_list v_list;
     // static gpc_polygon poly;
-    FGgpcPolyList polys_in, polys_clipped;
+    FGPolyList polys_in, polys_clipped;
 
 public:
 
@@ -98,13 +101,13 @@ public:
 
     // merge any slivers in the specified polygon with larger
     // neighboring polygons of higher priorigy
-    void merge_slivers(gpc_polygon *poly);
+    void merge_slivers( FGPolygon& poly);
     
     // Do actually clipping work
     bool clip_all(const point2d& min, const point2d& max);
 
     // return output poly list
-    inline FGgpcPolyList get_polys_clipped() const { return polys_clipped; }
+    inline FGPolyList get_polys_clipped() const { return polys_clipped; }
 };
 
 
diff --git a/Tools/Construct/Main/construct.hxx b/Tools/Construct/Main/construct.hxx
index 007dfeba7..44fbd2a26 100644
--- a/Tools/Construct/Main/construct.hxx
+++ b/Tools/Construct/Main/construct.hxx
@@ -69,7 +69,7 @@ private:
     FGBucket bucket;
 
     // clipped polygons (gpc format)
-    FGgpcPolyList clipped_polys;
+    FGPolyList clipped_polys;
 
     // raw node list (after triangulation)
     FGTriNodes tri_nodes;
@@ -125,8 +125,8 @@ public:
     inline void set_bucket( const FGBucket b ) { bucket = b; } 
 
     // clipped polygons
-    inline FGgpcPolyList get_clipped_polys() const { return clipped_polys; }
-    inline void set_clipped_polys( FGgpcPolyList p ) { clipped_polys = p; }
+    inline FGPolyList get_clipped_polys() const { return clipped_polys; }
+    inline void set_clipped_polys( FGPolyList p ) { clipped_polys = p; }
 
     // node list (after triangulation)
     inline FGTriNodes get_tri_nodes() const { return tri_nodes; }
diff --git a/Tools/Construct/Main/main.cxx b/Tools/Construct/Main/main.cxx
index 5923625e5..11c124ce8 100644
--- a/Tools/Construct/Main/main.cxx
+++ b/Tools/Construct/Main/main.cxx
@@ -164,7 +164,7 @@ void first_triangulate( FGConstruct& c, const FGArray& array,
 
     point_list corner_list = array.get_corner_node_list();
     point_list fit_list = array.get_fit_node_list();
-    FGgpcPolyList gpc_polys = c.get_clipped_polys();
+    FGPolyList gpc_polys = c.get_clipped_polys();
 
     cout << "ready to build node list and polygons" << endl;
     t.build( corner_list, fit_list, gpc_polys );
@@ -604,5 +604,3 @@ main(int argc, char **argv) {
 
     cout << "[Finished successfully]" << endl;
 }
-
-
diff --git a/Tools/Construct/Makefile.am b/Tools/Construct/Makefile.am
index 1456e2932..a2fd3255a 100644
--- a/Tools/Construct/Makefile.am
+++ b/Tools/Construct/Makefile.am
@@ -1,10 +1,10 @@
 SUBDIRS = \
 	Array \
+	Triangulate \
 	Clipper \
 	Combine \
 	GenOutput \
 	Match \
 	Parallel \
-	Triangulate \
 	Main
 
diff --git a/Tools/Construct/Triangulate/polygon.cxx b/Tools/Construct/Triangulate/polygon.cxx
index abad32cd0..3e25e89e5 100644
--- a/Tools/Construct/Triangulate/polygon.cxx
+++ b/Tools/Construct/Triangulate/polygon.cxx
@@ -107,24 +107,22 @@ void FGPolygon::calc_point_inside( const int contour,
     // min.y() starts greater than the biggest possible lat (degrees)
     min.sety( 100.0 );
 
-    int_list_iterator current, last;
+    point_list_iterator current, last;
     current = poly[contour].begin();
     last = poly[contour].end();
 
-    int counter = 0;
-    for ( ; current != last; ++current ) {
-	tmp = trinodes.get_node( *current );
+    for ( int i = 0; i < contour_size( contour ); ++i ) {
+	tmp = get_pt( contour, i );
 	if ( tmp.y() < min.y() ) {
 	    min = tmp;
-	    min_index = *current;
-	    min_node_index = counter;
+	    min_index = trinodes.find( min );
+	    min_node_index = i;
 
 	    // cout << "min index = " << *current 
 	    //      << " value = " << min_y << endl;
 	} else {
 	    // cout << "  index = " << *current << endl;
 	}
-	++counter;
     }
 
     cout << "min node index = " << min_node_index << endl;
@@ -136,17 +134,17 @@ void FGPolygon::calc_point_inside( const int contour,
     // fabs(slope)
 
     if ( min_node_index == 0 ) {
-	p1_index = poly[contour][1];
-	p2_index = poly[contour][poly[contour].size() - 1];
+	p1 = poly[contour][1];
+	p2 = poly[contour][poly[contour].size() - 1];
     } else if ( min_node_index == (int)(poly[contour].size()) - 1 ) {
-	p1_index = poly[contour][0];
-	p2_index = poly[contour][poly[contour].size() - 2];
+	p1 = poly[contour][0];
+	p2 = poly[contour][poly[contour].size() - 2];
     } else {
-	p1_index = poly[contour][min_node_index - 1];
-	p2_index = poly[contour][min_node_index + 1];
+	p1 = poly[contour][min_node_index - 1];
+	p2 = poly[contour][min_node_index + 1];
     }
-    p1 = trinodes.get_node( p1_index );
-    p2 = trinodes.get_node( p2_index );
+    p1_index = trinodes.find( p1 );
+    p2_index = trinodes.find( p2 );
 
     double s1 = fabs( slope(min, p1) );
     double s2 = fabs( slope(min, p2) );
@@ -175,10 +173,10 @@ void FGPolygon::calc_point_inside( const int contour,
 	for ( int j = 0; j < (int)(poly[i].size() - 1); ++j ) {
 	    // cout << "  p1 = " << poly[i][j] << " p2 = " 
 	    //      << poly[i][j+1] << endl;
-	    p1_index = poly[i][j];
-	    p2_index = poly[i][j+1];
-	    p1 = trinodes.get_node( p1_index );
-	    p2 = trinodes.get_node( p2_index );
+	    p1 = poly[i][j];
+	    p2 = poly[i][j+1];
+	    p1_index = trinodes.find( p1 );
+	    p2_index = trinodes.find( p2 );
 	
 	    if ( intersects(p1, p2, m.x(), &result) ) {
 		// cout << "intersection = " << result << endl;
@@ -191,10 +189,10 @@ void FGPolygon::calc_point_inside( const int contour,
 	}
 	// cout << "  p1 = " << poly[i][0] << " p2 = " 
 	//      << poly[i][poly[i].size() - 1] << endl;
-	p1_index = poly[i][0];
-	p2_index = poly[i][poly[i].size() - 1];
-	p1 = trinodes.get_node( p1_index );
-	p2 = trinodes.get_node( p2_index );
+	p1 = poly[i][0];
+	p2 = poly[i][poly[i].size() - 1];
+	p1_index = trinodes.find( p1 );
+	p2_index = trinodes.find( p2 );
 	if ( intersects(p1, p2, m.x(), &result) ) {
 	    // cout << "intersection = " << result << endl;
 	    if ( ( result.y() < p3.y() ) &&
@@ -220,34 +218,140 @@ void FGPolygon::calc_point_inside( const int contour,
 }
 
 
+//
 // wrapper functions for gpc polygon clip routines
+//
 
-// Difference
-FGPolygon polygon_diff(	const FGPolygon& subject, const FGPolygon& clip ) {
-    FGPolygon result;
-
-    gpc_polygon *poly = new gpc_polygon;
-    poly->num_contours = 0;
-    poly->contour = NULL;
-
+// Make a gpc_poly from an FGPolygon
+void make_gpc_poly( const FGPolygon& in, gpc_polygon *out ) {
     gpc_vertex_list v_list;
     v_list.num_vertices = 0;
     v_list.vertex = new gpc_vertex[FG_MAX_VERTICES];
 
-    // free allocated memory
-    gpc_free_polygon( poly );
+    cout << "making a gpc_poly" << endl;
+    cout << "  input contours = " << in.contours() << endl;
+
+    Point3D p;
+    // build the gpc_polygon structures
+    for ( int i = 0; i < in.contours(); ++i ) {
+	cout << "    contour " << i << " = " << in.contour_size( i ) << endl;
+	for ( int j = 0; j < in.contour_size( i ); ++j ) {
+	    p = in.get_pt( i, j );
+	    v_list.vertex[j].x = p.x();
+	    v_list.vertex[j].y = p.y();
+	}
+	v_list.num_vertices = in.contour_size( i );
+	gpc_add_contour( out, &v_list, in.get_hole_flag( i ) );
+    }
+
+    // free alocated memory
     delete v_list.vertex;
+}
+
+
+// Set operation type
+typedef enum {
+    POLY_DIFF,			// Difference
+    POLY_INT,			// Intersection
+    POLY_XOR,			// Exclusive or
+    POLY_UNION			// Union
+} clip_op;
+
+
+// Generic clipping routine
+FGPolygon polygon_clip( clip_op poly_op, const FGPolygon& subject, 
+			const FGPolygon& clip )
+{
+    FGPolygon result;
+
+    gpc_polygon *gpc_subject = new gpc_polygon;
+    gpc_subject->num_contours = 0;
+    gpc_subject->contour = NULL;
+    gpc_subject->hole = NULL;
+    make_gpc_poly( subject, gpc_subject );
+
+    gpc_polygon *gpc_clip = new gpc_polygon;
+    gpc_clip->num_contours = 0;
+    gpc_clip->contour = NULL;
+    gpc_clip->hole = NULL;
+    make_gpc_poly( clip, gpc_clip );
+
+    gpc_polygon *gpc_result = new gpc_polygon;
+    gpc_result->num_contours = 0;
+    gpc_result->contour = NULL;
+    gpc_result->hole = NULL;
+
+    gpc_op op;
+    if ( poly_op == POLY_DIFF ) {
+	op = GPC_DIFF;
+    } else if ( poly_op == POLY_INT ) {
+	op = GPC_INT;
+    } else if ( poly_op == POLY_XOR ) {
+	op = GPC_XOR;
+    } else if ( poly_op == POLY_UNION ) {
+	op = GPC_UNION;
+    } else {
+	cout << "Unknown polygon op, exiting." << endl;
+	exit(-1);
+    }
+
+    gpc_polygon_clip( op, gpc_subject, gpc_clip, gpc_result );
+
+    for ( int i = 0; i < gpc_result->num_contours; ++i ) {
+	// cout << "  processing contour = " << i << ", nodes = " 
+	//      << gpc_result->contour[i].num_vertices << ", hole = "
+	//      << gpc_result->hole[i] << endl;
+	
+	// sprintf(junkn, "g.%d", junkc++);
+	// junkfp = fopen(junkn, "w");
+
+	for ( int j = 0; j < gpc_result->contour[i].num_vertices; j++ ) {
+	    Point3D p( gpc_result->contour[i].vertex[j].x,
+		       gpc_result->contour[i].vertex[j].y,
+		       0 );
+	    // junkp = in_nodes.get_node( index );
+	    // fprintf(junkfp, "%.4f %.4f\n", junkp.x(), junkp.y());
+	    result.add_node(i, p);
+	    // cout << "  - " << index << endl;
+	}
+	// fprintf(junkfp, "%.4f %.4f\n", 
+	//    gpc_result->contour[i].vertex[0].x, 
+	//    gpc_result->contour[i].vertex[0].y);
+	// fclose(junkfp);
+
+	result.set_hole_flag( i, gpc_result->hole[i] );
+    }
+
+    // free allocated memory
+    gpc_free_polygon( gpc_subject );
+    gpc_free_polygon( gpc_clip );
+    gpc_free_polygon( gpc_result );
 
     return result;
 }
 
+
+// Difference
+FGPolygon polygon_diff(	const FGPolygon& subject, const FGPolygon& clip ) {
+    return polygon_clip( POLY_DIFF, subject, clip );
+}
+
 // Intersection
-FGPolygon polygon_int( const FGPolygon& subject, const FGPolygon& clip );
+FGPolygon polygon_int( const FGPolygon& subject, const FGPolygon& clip ) {
+    return polygon_clip( POLY_INT, subject, clip );
+}
+
 
 // Exclusive or
-FGPolygon polygon_xor( const FGPolygon& subject, const FGPolygon& clip );
+FGPolygon polygon_xor( const FGPolygon& subject, const FGPolygon& clip ) {
+    return polygon_clip( POLY_XOR, subject, clip );
+}
+
 
 // Union
-FGPolygon polygon_union( const FGPolygon& subject, const FGPolygon& clip );
+FGPolygon polygon_union( const FGPolygon& subject, const FGPolygon& clip ) {
+    return polygon_clip( POLY_UNION, subject, clip );
+}
+
 
 
diff --git a/Tools/Construct/Triangulate/polygon.hxx b/Tools/Construct/Triangulate/polygon.hxx
index ed1a6d515..cf5469710 100644
--- a/Tools/Construct/Triangulate/polygon.hxx
+++ b/Tools/Construct/Triangulate/polygon.hxx
@@ -44,7 +44,7 @@ FG_USING_STD(vector);
 #define FG_MAX_VERTICES 100000
 
 
-typedef vector < int_list > polytype;
+typedef vector < point_list > polytype;
 typedef polytype::iterator polytype_iterator;
 typedef polytype::const_iterator const_polytype_iterator;
 
@@ -64,10 +64,10 @@ public:
     ~FGPolygon( void );
 
     // Add the specified node (index) to the polygon
-    inline void add_node( int contour, int index ) {
+    inline void add_node( int contour, Point3D p ) {
 	if ( contour >= (int)poly.size() ) {
 	    // extend polygon
-	    int_list empty_contour;
+	    point_list empty_contour;
 	    empty_contour.clear();
 	    for ( int i = 0; i < contour - (int)poly.size() + 1; ++i ) {
 		poly.push_back( empty_contour );
@@ -75,7 +75,7 @@ public:
 		hole_list.push_back( 0 );
 	    }
 	}
-	poly[contour].push_back( index );
+	poly[contour].push_back( p );
     }
 
     // return size
@@ -85,7 +85,7 @@ public:
     }
 
     // return the ith polygon point index from the specified contour
-    inline int get_pt_index( int contour, int i ) const { 
+    inline Point3D get_pt( int contour, int i ) const { 
 	return poly[contour][i];
     }
 
diff --git a/Tools/Construct/Triangulate/triangle.cxx b/Tools/Construct/Triangulate/triangle.cxx
index 0379eeac3..d10babca7 100644
--- a/Tools/Construct/Triangulate/triangle.cxx
+++ b/Tools/Construct/Triangulate/triangle.cxx
@@ -38,7 +38,7 @@ FGTriangle::~FGTriangle( void ) {
 int 
 FGTriangle::build( const point_list& corner_list,
 		   const point_list& fit_list, 
-		   const FGgpcPolyList& gpc_polys )
+		   const FGPolyList& gpc_polys )
 {
     int debug_counter = 0;
     FGPolygon poly;
@@ -66,8 +66,8 @@ FGTriangle::build( const point_list& corner_list,
     }
 
     // next process the polygons
-    gpc_polygon *gpc_poly;
-    const_gpcpoly_iterator current, last;
+    FGPolygon gpc_poly;
+    const_poly_list_iterator current, last;
 
     // process polygons in priority order
     cout << "prepairing node list and polygons" << endl;
@@ -82,9 +82,9 @@ FGTriangle::build( const point_list& corner_list,
 	for ( ; current != last; ++current ) {
 	    gpc_poly = *current;
 	    cout << "processing a polygon, contours = " 
-		 << gpc_poly->num_contours << endl;
+		 << gpc_poly.contours() << endl;
 
-	    if (gpc_poly->num_contours <= 0 ) {
+	    if (gpc_poly.contours() <= 0 ) {
 		cout << "FATAL ERROR! no contours in this polygon" << endl;
 		exit(-1);
 	    }
@@ -93,22 +93,20 @@ FGTriangle::build( const point_list& corner_list,
 
 	    int j;
 
-	    for ( j = 0; j < gpc_poly->num_contours; ++j ) {
+	    for ( j = 0; j < gpc_poly.contours(); ++j ) {
 		cout << "  processing contour = " << j << ", nodes = " 
-		     << gpc_poly->contour[j].num_vertices << ", hole = "
-		     << gpc_poly->hole[j] << endl;
+		     << gpc_poly.contour_size( j ) << ", hole = "
+		     << gpc_poly.get_hole_flag( j ) << endl;
 
 		// sprintf(junkn, "g.%d", junkc++);
 		// junkfp = fopen(junkn, "w");
 
-		for ( int k = 0; k < gpc_poly->contour[j].num_vertices; k++ ) {
-		    Point3D p( gpc_poly->contour[j].vertex[k].x,
-			       gpc_poly->contour[j].vertex[k].y,
-			       0 );
+		for ( int k = 0; k < gpc_poly.contour_size( j ); k++ ) {
+		    Point3D p = gpc_poly.get_pt( j, k );
 		    index = in_nodes.unique_add( p );
 		    // junkp = in_nodes.get_node( index );
 		    // fprintf(junkfp, "%.4f %.4f\n", junkp.x(), junkp.y());
-		    poly.add_node(j, index);
+		    poly.add_node(j, p);
 		    // cout << "  - " << index << endl;
 		}
 		// fprintf(junkfp, "%.4f %.4f\n", 
@@ -116,10 +114,10 @@ FGTriangle::build( const point_list& corner_list,
 		//    gpc_poly->contour[j].vertex[0].y);
 		// fclose(junkfp);
 
-		poly.set_hole_flag( j, gpc_poly->hole[j] );
+		poly.set_hole_flag( j, gpc_poly.get_hole_flag( j ) );
 	    }
 
-	    for ( j = 0; j < gpc_poly->num_contours; ++j ) {
+	    for ( j = 0; j < gpc_poly.contours(); ++j ) {
 		poly.calc_point_inside( j, in_nodes );
 	    }
 
@@ -129,15 +127,12 @@ FGTriangle::build( const point_list& corner_list,
 		char pname[256];
 		sprintf(pname, "poly%02d-%02d-%02d", i, debug_counter, j);
 		FILE *fp = fopen( pname, "w" );
-		int index;
 		Point3D point;
 		for ( int k = 0; k < poly.contour_size( j ); ++k ) {
-		    index = poly.get_pt_index( j, k );
-		    point = in_nodes.get_node( index );
+		    point = poly.get_pt( j, k );
 		    fprintf( fp, "%.6f %.6f\n", point.x(), point.y() );
 		}
-		index = poly.get_pt_index( j, 0 );
-		point = in_nodes.get_node( index );
+		point = poly.get_pt( j, 0 );
 		fprintf( fp, "%.6f %.6f\n", point.x(), point.y() );
 		fclose(fp);
 
@@ -174,6 +169,7 @@ FGTriangle::build( const point_list& corner_list,
     // that is used by the "Triangle" lib.
 
     int i1, i2;
+    Point3D p1, p2;
     point_list node_list = in_nodes.get_node_list();
     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 	// cout << "area type = " << i << endl;
@@ -186,14 +182,18 @@ FGTriangle::build( const point_list& corner_list,
 	    poly = *tp_current;
 
 	    for ( int j = 0; j < (int)poly.contours(); ++j) {
-		for ( int k = 0; k < (int)(poly.contour_size(j)) - 1; ++k ) {
-		    i1 = poly.get_pt_index( j, k );
-		    i2 = poly.get_pt_index( j, k + 1 );
+		for ( int k = 0; k < (int)(poly.contour_size(j) - 1); ++k ) {
+		    p1 = poly.get_pt( j, k );
+		    p2 = poly.get_pt( j, k + 1 );
+		    i1 = in_nodes.find( p1 );
+		    i2 = in_nodes.find( p2 );
 		    // calc_line_params(i1, i2, &m, &b);
 		    in_segs.unique_divide_and_add( node_list, FGTriSeg(i1, i2) );
 		}
-		i1 = poly.get_pt_index( j, 0 );
-		i2 = poly.get_pt_index( j, poly.contour_size(j) - 1 );
+		p1 = poly.get_pt( j, 0 );
+		p2 = poly.get_pt( j, poly.contour_size(j) - 1 );
+		i1 = in_nodes.find( p1 );
+		i2 = in_nodes.find( p2 );
 		// calc_line_params(i1, i2, &m, &b);
 		in_segs.unique_divide_and_add( node_list, FGTriSeg(i1, i2) );
 	    }
diff --git a/Tools/Construct/Triangulate/triangle.hxx b/Tools/Construct/Triangulate/triangle.hxx
index d8df11a43..7c03175fa 100644
--- a/Tools/Construct/Triangulate/triangle.hxx
+++ b/Tools/Construct/Triangulate/triangle.hxx
@@ -79,7 +79,7 @@ public:
     // populate this class based on the specified gpc_polys list
     int build( const point_list& corner_list, 
 	       const point_list& fit_list,
-	       const FGgpcPolyList& gpc_polys );
+	       const FGPolyList& gpc_polys );
 
     // populate this class based on the specified gpc_polys list
     int rebuild( FGConstruct& c );
diff --git a/Tools/Construct/Triangulate/trinodes.cxx b/Tools/Construct/Triangulate/trinodes.cxx
index 9ff38f430..422584517 100644
--- a/Tools/Construct/Triangulate/trinodes.cxx
+++ b/Tools/Construct/Triangulate/trinodes.cxx
@@ -97,3 +97,27 @@ int FGTriNodes::course_add( const Point3D& p ) {
 }
 
 
+// Find the index of the specified point (compair to the same
+// tolerance as unique_add().  Returns -1 if not found.
+int FGTriNodes::find( const Point3D& p ) const {
+    const_point_list_iterator current, last;
+    int counter = 0;
+
+    // cout << p.x() << "," << p.y() << endl;
+
+    // see if point already exists
+    current = node_list.begin();
+    last = node_list.end();
+    for ( ; current != last; ++current ) {
+	if ( close_enough(p, *current) ) {
+	    // cout << "found an existing match!" << endl;
+	    return counter;
+	}
+	
+	++counter;
+    }
+
+    return -1;
+}
+
+
diff --git a/Tools/Construct/Triangulate/trinodes.hxx b/Tools/Construct/Triangulate/trinodes.hxx
index b39fcf1a0..23187892e 100644
--- a/Tools/Construct/Triangulate/trinodes.hxx
+++ b/Tools/Construct/Triangulate/trinodes.hxx
@@ -49,7 +49,7 @@ private:
 
     // return true of the two points are "close enough" as defined by
     // FG_PROXIMITY_EPSILON
-    bool close_enough( const Point3D& p, const Point3D& p );
+    bool close_enough( const Point3D& p, const Point3D& p ) const;
 
     // return true of the two points are "close enough" as defined by
     // FG_COURSE_EPSILON
@@ -76,7 +76,11 @@ public:
     // Use a course proximity check
     int course_add( const Point3D& p );
 
-    // return the master node list
+    // Find the index of the specified point (compair to the same
+    // tolerance as unique_add().  Returns -1 if not found.
+    int find( const Point3D& p ) const;
+
+     // return the master node list
     inline point_list get_node_list() const { return node_list; }
 
     // return the ith point
@@ -89,7 +93,9 @@ public:
 
 // return true of the two points are "close enough" as defined by
 // FG_PROXIMITY_EPSILON
-inline bool FGTriNodes::close_enough( const Point3D& p1, const Point3D& p2 ) {
+inline bool FGTriNodes::close_enough( const Point3D& p1, const Point3D& p2 )
+    const
+{
     if ( ( fabs(p1.x() - p2.x()) < FG_PROXIMITY_EPSILON ) &&
 	 ( fabs(p1.y() - p2.y()) < FG_PROXIMITY_EPSILON ) ) {
 	return true;