Workaround for point-in-polygon-calculation: Make "sure" that the point is sufficiently far away from the polygon border so that JRSTriangle can actually detect it is inside of the respective triangle.

src/Lib/Geometry/poly_support.cxx

@@ -4,6 +4,7 @@
 // Written by Curtis Olson, started October 1999.
 //
 // Copyright (C) 1999  Curtis L. Olson  - http://www.flightgear.org/~curt
+// Copyright (C) 2007, 2008  Ralf Gerlich - http://www.custom-scenery.org/
 //
 // This program is free software; you can redistribute it and/or
 // modify it under the terms of the GNU General Public License as
@@ -443,9 +444,98 @@ static void write_tree_element( TGContourNode *node, TGPolygon& p, int hole=0) {
     }
 }
 
+/*
+ * Disabled my older (and more robust) version of calc_point_inside() as
+ * TriangleJRS hat problems with the points produced.
+ *
+ * With the disabled version of calc_point_inside() polygons are sometimes
+ * attributed the wrong texture, which is most prominently visible when
+ * sea/ocean turns into land.
+ *
+ * The points generated by both the old and the new version are perfectly
+ * inside the polygons, however it seems that the "point-inside-triangle"-
+ * predicate of TriangleJRS is not sufficiently robust to detect this when the
+ * distance between the point and a triangle-edge is pretty small.
+ *
+ * This may be the case when a polygon has a sliver which has a high y-extent.
+ * Note that these are not removed by the cleanup after clipping, as the latter
+ * only removes polygons that explicitly are slivers, but not those in which
+ * a few points form a sliver and the rest of the polygon is pretty normal.
+ *
+ * Such slivers occur often when VMAP0-data is clipped against each other.
+ * The original VMAP0-input data is already clipped and inaccuracies on
+ * reading the data may lead to slivers occurring.
+ *
+ * The algorithm first sorts the vertices of the contour along the y-axis. An
+ * x-parallel line is laid between the two successive vertices which are
+ * farthest away from each other in y-direction.
+ *
+ * This is done to ensure that the line does not intersect any contour vertex,
+ * which spares us possibly shaky special case handling for such intersections.
+ *
+ * The x-line produces intersections with the polygon contour, which are
+ * then sorted along the x-axis. The intersections then separate sections of
+ * the x-line which are alternating between the inside and the outside of
+ * the polygon. The midpoint of the longest inside section is chosen as point
+ * in the polygon.
+ *
+ * In case of a sliver, the x-line is most probably placed through the sliver,
+ * as that has the largest y-extent. The inside-point is placed perfectly
+ * inside the sliver - as visual checks showed for the tiles in which the bug
+ * was detected - but is too near to the polygon contour for TriangleJRS to
+ * handle it correctly.
+ *
+ * The new algorithm does select the y-line and the x-parallel section
+ * differently. Instead maximizing the respective point distances, the
+ * x-line is placed between the first two contour vertices, that are more
+ * than a given limit apart from each other along the y-axis. Similarly,
+ * from the inside segments of that x-line the first segment longer than
+ * a given limit is chosen.
+ *
+ * This way it is possible to try again with a different line without
+ * performance issues, if the chosen line does not deliver a satisfactory
+ * point-in-polygon.
+ *
+ * Possible, the new algorithm is a bit faster than the old one, but it will
+ * fail on very small polygons. These should have been already removed after
+ * clipping.
+ *
+ * Another argument in favor of the new algorithm is that it will simply
+ * throw an exception if it can't deliver a proper point instead of silently
+ * generating a point TriangleJRS cannot handle.
+ *
+ * After all, fixing TriangleJRS would have been the better solution, but due
+ * to the complexity of that piece of software this wasn't feasible for now.
+ *
+ * - Ralf Gerlich
+ */
+#if 0
 // recurse the contour tree and build up the point inside list for
 // each contour/hole
 static void calc_point_inside( TGContourNode *node, TGPolygon &p ) {
+    /*
+     * Replaced the old algorithm for robustness.
+     *
+     * The old algorithm created a triangulation of the polygon and chose the
+     * midpoint from the largest triangle. In some situations, TriangleJRS
+     * choked on the polygons (while it does not in later stages of the
+     * triangulation).
+     *
+     * Therefore the old algorithm was replaced by a different and more robust
+     * one, taken from the shapefile-import-code in GRASS.
+     *
+     * The contour points are sorted in y-direction. The two successive
+     * points which are farthest away from each other in y-direction are chosen
+     * and an x-parallel line is laid between them. This line does not
+     * touch any of the contour points.
+     *
+     * Intersections of that line with the contour are calculated and sorted
+     * along the x-axis. The segments between successive pairs of intersections
+     * are alternating between inside and outside the polygon.
+     *
+     * The longest inside segment is chosen and its midpoint is selected as
+     * point inside the polygon.
+     */
     int contour_num = node->get_contour_num();
     // cout << "starting calc_point_inside() with contour = " << contour_num << endl;
 
@@ -458,19 +548,6 @@ static void calc_point_inside( TGContourNode *node, TGPolygon &p ) {
     if ( contour_num < 0 )
             return;
     
-    /*
-     * Find a line intersecting the contour and intersect it with the segments
-     * of our children. Sort the intersection points along the line. They then
-     * partition the line in IN/OUT parts. Find the longest segment and take
-     * its midpoint as point inside the contour.
-     */
-    
-    /*
-     * Try to find a line on which none of the contour points lie. For that,
-     * sort all contour points (also those of our direct children) by
-     * y-coordinate, find the two with the largest distance and take their
-     * center y-coordinate.
-     */
     point_list allpoints;
     
     collect_contour_points( node, p, allpoints );
@@ -548,6 +625,108 @@ static void calc_point_inside( TGContourNode *node, TGPolygon &p ) {
     
     p.set_point_inside( contour_num, Point3D(x, yline, -9999.0) );
 }
+#else
+// recurse the contour tree and build up the point inside list for
+// each contour/hole
+static void calc_point_inside( TGContourNode *node, TGPolygon &p ) {
+    int contour_num = node->get_contour_num();
+    // cout << "starting calc_point_inside() with contour = " << contour_num << endl;
+
+    for ( int i = 0; i < node->get_num_kids(); ++i ) {
+	if ( node->get_kid( i ) != NULL ) {
+	    calc_point_inside( node->get_kid( i ), p );
+	}
+    }
+    
+    if ( contour_num < 0 )
+            return;
+    
+    /*
+     * Find a line intersecting the contour and intersect it with the segments
+     * of our children. Sort the intersection points along the line. They then
+     * partition the line in IN/OUT parts. Find the longest segment and take
+     * its midpoint as point inside the contour.
+     */
+    
+    /*
+     * Try to find a line on which none of the contour points lie. For that,
+     * sort all contour points (also those of our direct children) by
+     * y-coordinate, find the two with the largest distance and take their
+     * center y-coordinate.
+     */
+    point_list allpoints;
+    
+    collect_contour_points( node, p, allpoints );
+    
+    for ( int i = 0; i < node->get_num_kids(); ++i ) {
+            if ( node->get_kid( i ) != NULL ) {
+                    collect_contour_points( node->get_kid( i ), p, allpoints );
+            }
+    }
+    
+    if ( allpoints.size() < 2 ) {
+            throw sg_exception("Polygon must have at least 2 contour points");
+    }
+    
+    sort(allpoints.begin(), allpoints.end(), Point3DOrdering(PY));
+
+    point_list::iterator point_it;
+    
+    point_it=allpoints.begin();
+    
+    double yline; // the y-location of the intersection line
+    
+    while ((++point_it) != allpoints.end()) {
+            double diff=point_it->y()-(point_it-1)->y();
+            if (diff<=8.0*SG_EPSILON) {
+                    continue;
+            }
+            yline=point_it->y()+diff/2.0;
+            
+            // cout << "calc_point_inside() " << allpoints.size() << " points ";
+            // copy(allpoints.begin(), allpoints.end(), ostream_iterator<Point3D>(cout, " "));
+            // cout << endl;
+        
+            // cout << "calc_point_inside() maxdiff=" << maxdiff << " yline=" << yline << endl;
+            
+            vector < double > xcuts;
+            
+            intersect_yline_with_contour( yline, node, p, xcuts );
+            
+            for ( int i = 0; i < node->get_num_kids(); ++i ) {
+                    if ( node->get_kid( i ) != NULL ) {
+                            intersect_yline_with_contour( yline, node->get_kid( i ), p, xcuts );
+                    }
+            }
+            
+            sort( xcuts.begin(), xcuts.end() );
+            
+            // cout << "calc_point_inside() " << xcuts.size() << " intersections ";
+            // copy(xcuts.begin(), xcuts.end(), ostream_iterator<double>(cout, " "));
+            // cout << endl;
+            
+            if ( xcuts.size() < 2 || (xcuts.size() % 2) != 0 ) {
+                    throw sg_exception("Geometric inconsistency in calc_point_inside()");
+            }
+            
+            for ( int i = 0; i < xcuts.size(); i+=2 ) {
+                    double x0 = xcuts[ i ];
+                    double x1 = xcuts[ i + 1 ];
+                    
+                    if ((x1-x0) > 8.0*SG_EPSILON) {
+                            double x = (x0+x1)/2.0;
+                            
+                            // cout << "calc_point_inside() found point inside x=" << x << " y=" << yline << endl;
+                            
+                            p.set_point_inside( contour_num, Point3D(x, yline, -9999.0) );
+                            
+                            return;
+                    }
+            }
+    }
+    throw sg_exception("No proper point inside found in calc_point_inside()");
+}
+#endif
 
 static void print_contour_tree( TGContourNode *node, string indent ) {
     cout << indent << node->get_contour_num() << endl;