Changes to try to track down some "point in a contour" bugs. The terrain

was stressing this code more that airports (which was where the code was
developed.)

src/BuildTiles/Triangulate/triangle.cxx

@@ -56,7 +56,6 @@ FGTriangle::build( const point_list& corner_list,
 
     // Point3D junkp;
     // int junkc = 0;
-    // char junkn[256];
     // FILE *junkfp;
 
     // traverse the dem corner and fit lists and gpc_polys building a
@@ -82,7 +81,8 @@ FGTriangle::build( const point_list& corner_list,
     for ( i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 	polylist[i].clear();
 
-	cout << "area type = " << i << endl;
+	cout << "area type = " << i << " polys = " << gpc_polys.polys[i].size() 
+	     << endl;
 	debug_counter = 0;
 	current = gpc_polys.polys[i].begin();
 	last = gpc_polys.polys[i].end();
@@ -103,8 +103,9 @@ FGTriangle::build( const point_list& corner_list,
 		     << gpc_poly.contour_size( j ) << ", hole = "
 		     << gpc_poly.get_hole_flag( j ) << endl;
 
-		// sprintf(junkn, "g.%d", junkc++);
-		// junkfp = fopen(junkn, "w");
+		char junkn[256];
+		sprintf(junkn, "c%d", j);
+		gpc_poly.write_contour( j, junkn );
 
 		for ( int k = 0; k < gpc_poly.contour_size( j ); k++ ) {
 		    Point3D p = gpc_poly.get_pt( j, k );
@@ -119,16 +120,26 @@ FGTriangle::build( const point_list& corner_list,
 		// fclose(junkfp);
 	    }
 
+	    /* if ( i == OceanArea ) {
+		cout << "temporary exit point" << endl;
+		exit(-1);
+	    } */
+
 	    // for each contour, calculate a point inside (but not
 	    // also inside any interior contours
 
 	    // new way
 
 	    // try to make sure our polygons aren't goofy
+	    /*
 	    gpc_poly = remove_dups( gpc_poly );
 	    gpc_poly = reduce_degeneracy( gpc_poly );
 	    gpc_poly = remove_dups( gpc_poly );
 	    gpc_poly = remove_bad_contours( gpc_poly );
+	    */
+
+	    cout << "after sanity checks, contours = " 
+		 << gpc_poly.contours() << endl;
 
 	    calc_points_inside( gpc_poly );
 

src/Lib/Geometry/poly_support.cxx

@@ -718,6 +718,7 @@ static Point3D point_inside_hole( point_list contour ) {
 // consideration
 static Point3D point_inside_contour( FGContourNode *node, const FGPolygon &p ) {
     int contour_num;
+    int i;
 
     FGPolygon hole_polys;
     hole_polys.erase();
@@ -726,7 +727,7 @@ static Point3D point_inside_contour( FGContourNode *node, const FGPolygon &p ) {
     hole_pts.clear();
 
     // build list of hole points
-    for ( int i = 0; i < node->get_num_kids(); ++i ) {
+    for ( i = 0; i < node->get_num_kids(); ++i ) {
         if ( node->get_kid( i ) != NULL ) {
 	    contour_num = node->get_kid(i)->get_contour_num();
 	    hole_pts.push_back( p.get_point_inside( contour_num ) );
@@ -743,22 +744,34 @@ static Point3D point_inside_contour( FGContourNode *node, const FGPolygon &p ) {
 	exit(-1);
     }
 
-    FGTriEle t = elelist[0];
+    // find the largest triangle in the group
+    double max_area = 0.0;
+    int biggest = 0;
+    for ( i = 0; i < (int)elelist.size(); ++i ) {
+	FGTriEle t = elelist[i];
+	Point3D p1 = out_pts[ t.get_n1() ];
+	Point3D p2 = out_pts[ t.get_n2() ];
+	Point3D p3 = out_pts[ t.get_n3() ];
+	double area = fabs(0.5 * ( p1.x() * p2.y() - p2.x() * p1.y() +
+				   p2.x() * p3.y() - p3.x() * p2.y() +
+				   p3.x() * p1.y() - p1.x() * p3.y() ));
+	if ( area > max_area ) {
+	    max_area = area;
+	    biggest = i;
+	}
+    }
+
+    // find center point of largest triangle
+    cout << "biggest = " << biggest + 1 << " out of " << elelist.size() << endl;
+    FGTriEle t = elelist[biggest];
     contour_num = node->get_contour_num();
     Point3D p1 = out_pts[ t.get_n1() ];
     Point3D p2 = out_pts[ t.get_n2() ];
     Point3D p3 = out_pts[ t.get_n3() ];
-    cout << "  " << p1 << endl << "  " << p2 << endl << "  " << p3 << endl;
+    cout << "hole tri = " << p1 << endl << "  " << p2 << endl << "  " << p3 << endl;
 
-#if 0
-    // old
-    Point3D m1 = ( p1 + p2 ) / 2;
-    Point3D m2 = ( p1 + p3 ) / 2;
-    Point3D center = ( m1 + m2 ) / 2;
-#else
     // new
     Point3D center = ( p1 + p2 + p3 ) / 3;
-#endif
 
     return center;
 
@@ -807,11 +820,12 @@ static void build_contour_tree( FGContourNode *node,
 
     // see if we are building on a hole or note
     bool flag;
-    if ( node->get_contour_num() > 0 ) {
+    if ( node->get_contour_num() >= 0 ) {
 	flag = p.get_hole_flag( node->get_contour_num() );
     } else {
 	flag = true;
     }
+    cout << "  hole flag = " << flag << endl;
 
     // add all remaining hole/non-hole contours as children of the
     // current node if they are inside of it.
@@ -848,7 +862,12 @@ static void build_contour_tree( FGContourNode *node,
 	for ( int j = 0; j < node->get_num_kids(); ++j ) {
 	    if ( i != j ) {
 		if ( (node->get_kid(i) != NULL)&&(node->get_kid(j) != NULL) ) {
-		    if ( p.is_inside( i, j ) ) {
+		    int A = node->get_kid( i ) -> get_contour_num();
+		    int B = node->get_kid( j ) -> get_contour_num();
+		    if ( p.is_inside( A, B ) ) {
+			// p.write_contour( i, "a" );
+			// p.write_contour( j, "b" );
+			// exit(-1);
 		        // need to remove contour j from the kid list
 		        avail[ node->get_kid( i ) -> get_contour_num() ] = 1;
 		        node->remove_kid( i );
@@ -1147,6 +1166,8 @@ FGPolygon remove_bad_contours( const FGPolygon &poly ) {
 	    // good
 	    int flag = poly.get_hole_flag( i );
 	    result.add_contour( contour, flag );
+	} else {
+	    cout << "tossing a bad contour" << endl;
 	}
     }
 

src/Lib/Polygon/polygon.cxx

@@ -85,7 +85,7 @@ static double calc_angle(point2d a, point2d b, point2d c) {
 
 double FGPolygon::area_contour( const int contour ) const {
     // 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
+    // where i=k is defined as i=0
 
     point_list c = poly[contour];
     int size = c.size();
@@ -180,7 +180,7 @@ void FGPolygon::shift( double lon, double lat ) {
 
 
 // output
-void FGPolygon::write( const string& file ) {
+void FGPolygon::write( const string& file ) const {
     FILE *fp = fopen( file.c_str(), "w" );
     
     for ( int i = 0; i < (int)poly.size(); ++i ) {
@@ -194,6 +194,18 @@ void FGPolygon::write( const string& file ) {
 }
 
 
+// output
+void FGPolygon::write_contour( const int contour, const string& file ) const {
+    FILE *fp = fopen( file.c_str(), "w" );
+    
+    for ( int j = 0; j < (int)poly[contour].size(); ++j ) {
+	fprintf(fp, "%.6f %.6f\n", poly[contour][j].x(), poly[contour][j].y());
+    }
+
+    fclose(fp);
+}
+
+
 //
 // wrapper functions for gpc polygon clip routines
 //

src/Lib/Polygon/polygon.hxx

@@ -162,7 +162,10 @@ public:
     bool is_inside( int a, int b ) const;
 
     // output
-    void write( const string& file );
+    void write( const string& file ) const;
+
+    // output
+    void write_contour( const int contour, const string& file ) const;
 };