Sliver detection and elimination:

After I clip a polygon against the more important stuff, I
check the area and minimum interior angle of each of it's contours.
If both (area and min interior angle) are below some threshhold I'm
calling it a sliver.  Then I go and look for another polygon such that
result = polygon_union( poly, sliver ) yields a result that doesn't
have anymore contours than the original poly.  This means the sliver
is adjacent to poly.  So I replace poly with result (the union) and
keep going.

Tools/Construct/Clipper/clipper.cxx

@@ -152,17 +152,110 @@ bool FGClipper::load_polys(const string& path) {
 }
 
 
-// merge any slivers in the specified polygon with larger
+// remove any slivers from in polygon and move them to out polygon.
-// neighboring polygons of higher priorigy
+void FGClipper::move_slivers( FGPolygon& in, FGPolygon& out ) {
-void FGClipper::merge_slivers(FGPolygon& poly) {
+    cout << "Begin move slivers" << endl;
-    cout << "Begin merge slivers" << endl;
     // traverse each contour of the polygon and attempt to identify
     // likely slivers
-    for ( int i = 0; i < poly.contours(); ++i ) {
+
+    out.erase();
+
+    double angle_cutoff = 10.0 * DEG_TO_RAD;
+    double area_cutoff = 0.00001;
+    double min_angle;
+    double area;
+
+    point_list contour;
+    int hole_flag;
+
+    // process contours in reverse order so deleting a contour doesn't
+    // foul up our sequence
+    for ( int i = in.contours() - 1; i >= 0; --i ) {
 	cout << "contour " << i << endl;
-	for (int j = 0; j < poly.contour_size( i ); ++j ) {
+
-	    // cout << poly->contour[i].vertex[j].x << ","
+	min_angle = in.minangle_contour( i );
-	    //      << poly->contour[i].vertex[j].y << endl;
+	area = in.area_contour( i );
+
+	cout << "  min_angle (rad) = " 
+	     << min_angle << endl;
+	cout << "  min_angle (deg) = " 
+	     << min_angle * 180.0 / FG_PI << endl;
+	cout << "  area = " << area << endl;
+
+	if ( (min_angle < angle_cutoff) && (area < area_cutoff)  ) {
+	    cout << "      WE THINK IT'S A SLIVER!" << endl;
+
+	    // check if this is a hole
+	    hole_flag = in.get_hole_flag( i );
+
+	    if ( hole_flag ) {
+		// just delete/eliminate/remove sliver holes
+		in.delete_contour( i );
+	    } else {
+		// move sliver contour to out polygon
+		contour = in.get_contour( i );
+		in.delete_contour( i );
+		out.add_contour( contour, hole_flag );
+	    }
+	}
+    }
+}
+
+
+// for each sliver contour, see if a union with another polygon yields
+// a polygon with no increased contours (i.e. the sliver is adjacent
+// and can be merged.)  If so, replace the clipped polygon with the
+// new polygon that has the sliver merged in.
+void FGClipper::merge_slivers( FGPolyList& clipped, FGPolygon& slivers ) {
+    FGPolygon poly, result, sliver;
+    point_list contour;
+    int original_contours, result_contours;
+    bool done;
+
+    for ( int i = 0; i < slivers.contours(); ++i ) {
+	cout << "Merging sliver = " << i << endl;
+
+	// make the sliver polygon
+	contour = slivers.get_contour( i );
+	sliver.erase();
+	sliver.add_contour( contour, 0 );
+	done = false;
+
+	for ( int area = 0; area < FG_MAX_AREA_TYPES && !done; ++area ) {
+	    cout << "  testing area = " << area << " with " 
+		 << clipped.polys[area].size() << " polys" << endl;
+	    for ( int j = 0; 
+		  j < (int)clipped.polys[area].size() && !done;
+		  ++j )
+	    {
+		cout << "  polygon = " << j << endl;
+
+		poly = clipped.polys[area][j];
+		original_contours = poly.contours();
+		result = polygon_union( poly, sliver );
+		result_contours = result.contours();
+
+		if ( original_contours == result_contours ) {
+		    cout << "    FOUND a poly to merge the sliver with" << endl;
+		    clipped.polys[area][j] = result;
+		    done = true;
+		    // poly.write("orig");
+		    // sliver.write("sliver");
+		    // result.write("result");
+		    // cout << "press return: ";
+		    // string input;
+		    // cin >> input;
+		} else {
+		    cout << "    poly not a match" << endl;
+		    cout << "    original = " << original_contours
+			 << " result = " << result_contours << endl;
+		    cout << "    sliver = " << endl;
+		    for ( int k = 0; k < (int)contour.size(); ++k ) {
+			cout << "      " << contour[k].x() << ", "
+			     << contour[k].y() << endl;
+		    }
+		}
+	    }
 	}
     }
 }
@@ -171,7 +264,7 @@ void FGClipper::merge_slivers(FGPolygon& poly) {
 // Do actually clipping work
 bool FGClipper::clip_all(const point2d& min, const point2d& max) {
     FGPolygon accum, result_union, tmp;
-    FGPolygon result_diff, remains;
+    FGPolygon result_diff, slivers, remains;
     // gpcpoly_iterator current, last;
 
     FG_LOG( FG_CLIPPER, FG_INFO, "Running master clipper" );
@@ -258,10 +351,23 @@ 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.contours() > 0 ) {
-		// merge any slivers with larger neighboring polygons
+		// move slivers from result_diff polygon to slivers polygon
-		merge_slivers(result_diff);
+		move_slivers(result_diff, slivers);
+		cout << "  After sliver move:" << endl;
+		cout << "    result_diff = " << result_diff.contours() << endl;
+		cout << "    slivers = " << slivers.contours() << endl;
 
-		polys_clipped.polys[i].push_back(result_diff);
+		// merge any slivers with previously clipped
+		// neighboring polygons
+		if ( slivers.contours() > 0 ) {
+		    merge_slivers(polys_clipped, slivers);
+		}
+
+		// add the sliverless result polygon (from after the
+		// move_slivers) to the clipped polys list
+		if ( result_diff.contours() > 0  ) {
+		    polys_clipped.polys[i].push_back(result_diff);
+		}
 
 		// char filename[256];
 		// sprintf(filename, "next-result-%02d", count++);
@@ -282,7 +388,22 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
     remains = polygon_diff( polys_in.safety_base, accum );
 
     if ( remains.contours() > 0 ) {
-	polys_clipped.polys[(int)OceanArea].push_back(remains);
+	cout << "remains contours = " << remains.contours() << endl;
+	// move slivers from remains polygon to slivers polygon
+	move_slivers(remains, slivers);
+	cout << "  After sliver move:" << endl;
+	cout << "    remains = " << remains.contours() << endl;
+	cout << "    slivers = " << slivers.contours() << endl;
+
+	// merge any slivers with previously clipped
+	// neighboring polygons
+	if ( slivers.contours() > 0 ) {
+	    merge_slivers(polys_clipped, slivers);
+	}
+
+	if ( remains.contours() > 0 ) {
+	    polys_clipped.polys[(int)OceanArea].push_back(remains);
+	}
     }
 
 #if 0

Tools/Construct/Clipper/clipper.hxx

@@ -45,6 +45,7 @@ extern "C" {
 }
 #endif
 
+#include <Main/construct_types.hxx>
 #include <Triangulate/polygon.hxx>
 
 #include STL_STRING
@@ -64,12 +65,6 @@ FG_USING_STD(vector);
 // #define FG_MAX_VERTICES 100000
 
 
-class point2d {
-public:
-    double x, y;
-};
-
-
 class FGPolyList {
 public:
     poly_list polys[FG_MAX_AREA_TYPES];
@@ -99,9 +94,15 @@ public:
     // Load a polygon definition file
     bool load_polys(const string& path);
 
-    // merge any slivers in the specified polygon with larger
+    // remove any slivers from in polygon and move them to out
-    // neighboring polygons of higher priorigy
+    // polygon.
-    void merge_slivers( FGPolygon& poly);
+    void move_slivers( FGPolygon& in, FGPolygon& out );
+
+    // for each sliver contour, see if a union with another polygon
+    // yields a polygon with no increased contours (i.e. the sliver is
+    // adjacent and can be merged.)  If so, replace the clipped
+    // polygon with the new polygon that has the sliver merged in.
+    void merge_slivers( FGPolyList& clipped, FGPolygon& slivers );
     
     // Do actually clipping work
     bool clip_all(const point2d& min, const point2d& max);

Tools/Construct/Main/construct_types.hxx

@@ -48,6 +48,12 @@ typedef point_list::iterator point_list_iterator;
 typedef point_list::const_iterator const_point_list_iterator;
 
 
+class point2d {
+public:
+    double x, y;
+};
+
+
 #endif // _CONSTRUCT_TYPES_HXX
 
 

Tools/Construct/Triangulate/polygon.cxx

@@ -57,6 +57,33 @@ static double slope( const Point3D& p0, const Point3D& p1 ) {
 }
 
 
+// Calculate theta of angle (a, b, c)
+static double calc_angle(point2d a, point2d b, point2d c) {
+    point2d u, v;
+    double udist, vdist, uv_dot, tmp;
+
+    // u . v = ||u|| * ||v|| * cos(theta)
+
+    u.x = b.x - a.x;
+    u.y = b.y - a.y;
+    udist = sqrt( u.x * u.x + u.y * u.y );
+    // printf("udist = %.6f\n", udist);
+
+    v.x = b.x - c.x;
+    v.y = b.y - c.y;
+    vdist = sqrt( v.x * v.x + v.y * v.y );
+    // printf("vdist = %.6f\n", vdist);
+
+    uv_dot = u.x * v.x + u.y * v.y;
+    // printf("uv_dot = %.6f\n", uv_dot);
+
+    tmp = uv_dot / (udist * vdist);
+    // printf("tmp = %.6f\n", tmp);
+
+    return acos(tmp);
+}
+
+
 // Given a line segment specified by two endpoints p1 and p2, return
 // the y value of a point on the line that intersects with the
 // verticle line through x.  Return true if an intersection is found,
@@ -218,6 +245,81 @@ void FGPolygon::calc_point_inside( const int contour,
 }
 
 
+// return the perimeter of a contour (assumes simple polygons,
+// i.e. non-self intersecting.)
+double FGPolygon::area_contour( const int contour ) {
+    // area = 1/2 * sum[i = 0 to k-1][x(i)*y(i+1) - x(i+1)*y(i)]
+    // where k is defined as 0
+
+    point_list c = poly[contour];
+    int size = c.size();
+    double sum = 0.0;
+
+    for ( int i = 0; i < size; ++i ) {
+	sum += c[(i+1)%size].x() * c[i].y() - c[i].x() * c[(i+1)%size].y();
+    }
+
+    return sum / 2.0;
+}
+
+
+// return the smallest interior angle of the polygon
+double FGPolygon::minangle_contour( const int contour ) {
+    point_list c = poly[contour];
+    int size = c.size();
+    int p1_index, p2_index, p3_index;
+    point2d p1, p2, p3;
+    double angle;
+    double min_angle = 2.0 * FG_PI;
+
+    for ( int i = 0; i < size; ++i ) {
+	p1_index = i - 1;
+	if ( p1_index < 0 ) {
+	    p1_index += size;
+	}
+
+	p2_index = i;
+
+	p3_index = i + 1;
+	if ( p3_index >= size ) {
+	    p3_index -= size;
+	}
+
+	p1.x = c[p1_index].x();
+	p1.y = c[p1_index].y();
+
+	p2.x = c[p2_index].x();
+	p2.y = c[p2_index].y();
+
+	p3.x = c[p3_index].x();
+	p3.y = c[p3_index].y();
+
+	angle = calc_angle( p1, p2, p3 );
+
+	if ( angle < min_angle ) {
+	    min_angle = angle;
+	}
+    }
+
+    return min_angle;
+}
+
+
+// output
+void FGPolygon::write( const string& file ) {
+    FILE *fp = fopen( file.c_str(), "w" );
+    
+    for ( int i = 0; i < (int)poly.size(); ++i ) {
+	for ( int j = 0; j < (int)poly[i].size(); ++j ) {
+	    fprintf(fp, "%.6f %.6f\n", poly[i][j].x(), poly[i][j].y());
+	}
+	fprintf(fp, "%.6f %.6f\n", poly[i][0].x(), poly[i][0].y());
+    }
+
+    fclose(fp);
+}
+
+
 //
 // wrapper functions for gpc polygon clip routines
 //
@@ -228,13 +330,13 @@ void make_gpc_poly( const FGPolygon& in, gpc_polygon *out ) {
     v_list.num_vertices = 0;
     v_list.vertex = new gpc_vertex[FG_MAX_VERTICES];
 
-    cout << "making a gpc_poly" << endl;
+    // cout << "making a gpc_poly" << endl;
-    cout << "  input contours = " << in.contours() << 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;
+	// 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();

Tools/Construct/Triangulate/polygon.hxx

@@ -32,12 +32,14 @@
 
 #include <Include/compiler.h>
 
+#include <string>
 #include <vector>
 
 #include <Main/construct_types.hxx>
 
 #include "trinodes.hxx"
 
+FG_USING_STD(string);
 FG_USING_STD(vector);
 
 
@@ -63,6 +65,26 @@ public:
     FGPolygon( void );
     ~FGPolygon( void );
 
+    // Add a contour
+    inline void add_contour( const point_list contour, const int hole_flag ) {
+	poly.push_back( contour );
+	hole_list.push_back( hole_flag );
+    }
+
+    // Get a contour
+    inline point_list get_contour( const int i ) const {
+	return poly[i];
+    }
+
+    // Delete a contour
+    inline void delete_contour( const int i ) {
+	polytype_iterator start_poly = poly.begin();
+	poly.erase( start_poly + i );
+
+	int_list_iterator start_hole = hole_list.begin();
+	hole_list.erase( start_hole + i );
+    }
+
     // Add the specified node (index) to the polygon
     inline void add_node( int contour, Point3D p ) {
 	if ( contour >= (int)poly.size() ) {
@@ -104,8 +126,20 @@ public:
 	hole_list[contour] = flag;
     }
 
+    // erase
     inline void erase() { poly.clear(); }
 
+    // informational
+
+    // return the area of a contour (assumes simple polygons,
+    // i.e. non-self intersecting.)
+    double area_contour( const int contour );
+
+    // return the smallest interior angle of the polygon
+    double minangle_contour( const int contour );
+
+    // output
+    void write( const string& file );
 };
 
 

Tools/Construct/Triangulate/triangle.cxx

@@ -168,11 +168,12 @@ FGTriangle::build( const point_list& corner_list,
     // traverse the polygon lists and build the segment (edge) list
     // that is used by the "Triangle" lib.
 
+    cout << "building segment list" << endl;
     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;
+	cout << "area type = " << i << endl;
 	poly_list_iterator tp_current, tp_last;
 	tp_current = polylist[i].begin();
 	tp_last = polylist[i].end();
@@ -180,7 +181,8 @@ FGTriangle::build( const point_list& corner_list,
 	// process each polygon in list
 	for ( ; tp_current != tp_last; ++tp_current ) {
 	    poly = *tp_current;
-
+	    cout << "  processing a polygon with contours = " 
+		 << poly.contours() << endl;
 	    for ( int j = 0; j < (int)poly.contours(); ++j) {
 		for ( int k = 0; k < (int)(poly.contour_size(j) - 1); ++k ) {
 		    p1 = poly.get_pt( j, k );