Cleaned up some debugging output.

Added adjacent triangle area weighting for calculating vertices of nodes.
Check for super small or degenerate triangles which could blow up our face
  normal calculations.

src/BuildTiles/Clipper/clipper.cxx

@@ -109,7 +109,7 @@ bool FGClipper::load_polys(const string& path) {
 	    in >> startx;
 	    in >> starty;
 	    p = Point3D(startx, starty, 0.0);
-	    cout << "poly pt = " << p << endl;
+	    // cout << "poly pt = " << p << endl;
 	    poly.add_node( i, p );
 	    FG_LOG( FG_CLIPPER, FG_BULK, "0 = " 
 		    << startx << ", " << starty );
@@ -118,7 +118,7 @@ bool FGClipper::load_polys(const string& path) {
 		in >> x;
 		in >> y;
 		p = Point3D( x, y, 0.0 );
-		cout << "poly pt = " << p << endl;
+		// cout << "poly pt = " << p << endl;
 		poly.add_node( i, p );
 	    }
 
@@ -130,7 +130,7 @@ bool FGClipper::load_polys(const string& path) {
 		// last point same as first, discard
 	    } else {
 		p = Point3D( lastx, lasty, 0.0 );
-		cout << "poly pt = " << p << endl;
+		// cout << "poly pt = " << p << endl;
 		poly.add_node( i, p );
 	    }
 
@@ -169,10 +169,11 @@ void FGClipper::add_poly (int area, const FGPolygon &poly)
 
 // remove any slivers from in polygon and move them to out polygon.
 void FGClipper::move_slivers( FGPolygon& in, FGPolygon& out ) {
-    cout << "Begin move slivers" << endl;
     // traverse each contour of the polygon and attempt to identify
     // likely slivers
 
+    // cout << "Begin move slivers" << endl;
+
     out.erase();
 
     double angle_cutoff = 10.0 * DEG_TO_RAD;
@@ -186,28 +187,28 @@ void FGClipper::move_slivers( FGPolygon& in, FGPolygon& out ) {
     // 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;
+	// cout << "contour " << i << endl;
 
 	min_angle = in.minangle_contour( i );
 	area = in.area_contour( i );
 
-	cout << "  min_angle (rad) = " 
+	/* cout << "  min_angle (rad) = " 
 	     << min_angle << endl;
 	cout << "  min_angle (deg) = " 
 	     << min_angle * 180.0 / FG_PI << endl;
-	cout << "  area = " << area << endl;
+	cout << "  area = " << area << endl; */
 
 	if ( ((min_angle < angle_cutoff) && (area < area_cutoff)) ||
 	     ( area < area_cutoff / 10.0) )
 	{
-	    cout << "      WE THINK IT'S A SLIVER!" << endl;
+	    // 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
-		cout << "just deleting a sliver hole" << endl;
+		// cout << "just deleting a sliver hole" << endl;
 		in.delete_contour( i );
 	    } else {
 		// move sliver contour to out polygon
@@ -231,7 +232,7 @@ void FGClipper::merge_slivers( FGPolyList& clipped, FGPolygon& slivers ) {
     bool done;
 
     for ( int i = 0; i < slivers.contours(); ++i ) {
-	cout << "Merging sliver = " << i << endl;
+	// cout << "Merging sliver = " << i << endl;
 
 	// make the sliver polygon
 	contour = slivers.get_contour( i );
@@ -246,13 +247,13 @@ void FGClipper::merge_slivers( FGPolyList& clipped, FGPolygon& slivers ) {
 		continue;
 	    }
 
-	    cout << "  testing area = " << area << " with " 
+	    // cout << "  testing area = " << area << " with " 
-		 << clipped.polys[area].size() << " polys" << endl;
+	    //      << clipped.polys[area].size() << " polys" << endl;
 	    for ( int j = 0; 
 		  j < (int)clipped.polys[area].size() && !done;
 		  ++j )
 	    {
-		cout << "  polygon = " << j << endl;
+		// cout << "  polygon = " << j << endl;
 
 		poly = clipped.polys[area][j];
 		original_contours = poly.contours();
@@ -260,7 +261,7 @@ void FGClipper::merge_slivers( FGPolyList& clipped, FGPolygon& slivers ) {
 		result_contours = result.contours();
 
 		if ( original_contours == result_contours ) {
-		    cout << "    FOUND a poly to merge the sliver with" << endl;
+		    // cout << "    FOUND a poly to merge the sliver with" << endl;
 		    clipped.polys[area][j] = result;
 		    done = true;
 		    // poly.write("orig");
@@ -270,19 +271,19 @@ void FGClipper::merge_slivers( FGPolyList& clipped, FGPolygon& slivers ) {
 		    // string input;
 		    // cin >> input;
 		} else {
-		    cout << "    poly not a match" << endl;
+		    /* cout << "    poly not a match" << endl;
 		    cout << "    original = " << original_contours
 			 << " result = " << result_contours << endl;
-		    cout << "    sliver = " << endl;
+		    cout << "    sliver = " << endl; */
 		    for ( int k = 0; k < (int)contour.size(); ++k ) {
-			cout << "      " << contour[k].x() << ", "
+			// cout << "      " << contour[k].x() << ", "
-			     << contour[k].y() << endl;
+			//      << contour[k].y() << endl;
 		    }
 		}
 	    }
 	}
 	if ( !done ) {
-	    cout << "no suitable polys found for sliver merge" << endl;
+	    // cout << "no suitable polys found for sliver merge" << endl;
 	}
     }
 }
@@ -430,9 +431,9 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
 	    if ( result_diff.contours() > 0 ) {
 		// move slivers from result_diff polygon to slivers polygon
 		move_slivers(result_diff, slivers);
-		cout << "  After sliver move:" << endl;
+		// cout << "  After sliver move:" << endl;
-		cout << "    result_diff = " << result_diff.contours() << endl;
+		// cout << "    result_diff = " << result_diff.contours() << endl;
-		cout << "    slivers = " << slivers.contours() << endl;
+		// cout << "    slivers = " << slivers.contours() << endl;
 
 		// merge any slivers with previously clipped
 		// neighboring polygons
@@ -465,12 +466,12 @@ bool FGClipper::clip_all(const point2d& min, const point2d& max) {
     remains = polygon_diff( polys_in.safety_base, accum );
 
     if ( remains.contours() > 0 ) {
-	cout << "remains contours = " << remains.contours() << endl;
+	// cout << "remains contours = " << remains.contours() << endl;
 	// move slivers from remains polygon to slivers polygon
 	move_slivers(remains, slivers);
-	cout << "  After sliver move:" << endl;
+	// cout << "  After sliver move:" << endl;
-	cout << "    remains = " << remains.contours() << endl;
+	// cout << "    remains = " << remains.contours() << endl;
-	cout << "    slivers = " << slivers.contours() << endl;
+	// cout << "    slivers = " << slivers.contours() << endl;
 
 	// merge any slivers with previously clipped
 	// neighboring polygons

src/BuildTiles/GenOutput/genobj.cxx

@@ -353,6 +353,7 @@ int FGGenOutput::write( FGConstruct &c ) {
 
     // write nodes
     point_list wgs84_nodes = c.get_wgs84_nodes();
+    cout << "writing nodes = " << wgs84_nodes.size() << endl;
     fprintf(fp, "# vertex list\n");
     const_point_list_iterator w_current = wgs84_nodes.begin();
     const_point_list_iterator w_last = wgs84_nodes.end();
@@ -364,6 +365,7 @@ int FGGenOutput::write( FGConstruct &c ) {
     
     // write vertex normals
     point_list point_normals = c.get_point_normals();
+    cout << "writing normals = " << point_normals.size() << endl;
     fprintf(fp, "# vertex normal list\n");
     const_point_list_iterator n_current = point_normals.begin();
     const_point_list_iterator n_last = point_normals.end();

src/BuildTiles/Main/main.cxx

@@ -46,6 +46,7 @@
 #include <simgear/bucket/newbucket.hxx>
 #include <simgear/debug/logstream.hxx>
 
+#include <Geometry/poly_support.hxx>
 #include <Array/array.hxx>
 #include <Clipper/clipper.hxx>
 #include <GenOutput/genobj.hxx>
@@ -362,7 +363,14 @@ static void second_triangulate( FGConstruct& c, FGTriangle& t ) {
     cout << "done re building node list and polygons" << endl;
 
     cout << "ready to do second triangulation" << endl;
+
+    cout << "  (pre) nodes = " << c.get_tri_nodes().size() << endl;
+    cout << "  (pre) normals = " << c.get_point_normals().size() << endl;
+
     t.run_triangulate( c.get_angle(), 2 );
+
+    cout << "  (post) nodes = " << t.get_out_nodes().size() << endl;
+
     cout << "finished second triangulation" << endl;
 }
 
@@ -387,7 +395,8 @@ static void fix_point_heights( FGConstruct& c, const FGArray& array ) {
     for ( i = 0; i < (int)raw_nodes.size(); ++i ) {
 	z = array.interpolate_altitude( raw_nodes[i].x() * 3600.0, 
 					raw_nodes[i].y() * 3600.0 );
-	cout << "  old z = " << raw_nodes[i].z() << "  new z = " << z << endl;
+	// cout << "  old z = " << raw_nodes[i].z() << "  new z = " << z
+	//      << endl;
 	if ( raw_nodes[i].z() != z ) {
 	    cout << "    DIFFERENT" << endl;
 	}
@@ -543,6 +552,18 @@ static belongs_to_list gen_node_ele_lookup_table( FGConstruct& c ) {
 }
 
 
+// caclulate the area for the specified triangle face
+static double tri_ele_area( const FGConstruct& c, const FGTriEle tri ) {
+    point_list nodes = c.get_geod_nodes();
+
+    Point3D p1 = nodes[ tri.get_n1() ];
+    Point3D p2 = nodes[ tri.get_n2() ];
+    Point3D p3 = nodes[ tri.get_n3() ];
+
+    return triangle_area( p1, p2, p3 );
+}
+
+
 // caclulate the normal for the specified triangle face
 static Point3D calc_normal( FGConstruct& c, int i ) {
     sgVec3 v1, v2, normal;
@@ -554,11 +575,52 @@ static Point3D calc_normal( FGConstruct& c, int i ) {
     Point3D p2 = wgs84_nodes[ tri_elements[i].get_n2() ];
     Point3D p3 = wgs84_nodes[ tri_elements[i].get_n3() ];
 
-    v1[0] = p2.x() - p1.x(); v1[1] = p2.y() - p1.y(); v1[2] = p2.z() - p1.z();
+    // do some sanity checking.  With the introduction of landuse
-    v2[0] = p3.x() - p1.x(); v2[1] = p3.y() - p1.y(); v2[2] = p3.z() - p1.z();
+    // areas, we can get some long skinny triangles that blow up our
+    // "normal" calculations here.  Let's check for really small
+    // triangle areas and check if one dimension of the triangle
+    // coordinates is nearly coincident.  If so, assign the "default"
+    // normal of straight up.
 
-    sgVectorProductVec3( normal, v1, v2 );
+    bool degenerate = false;
-    sgNormalizeVec3( normal );
+    const double area_eps = 1.0e-12;
+    double area = tri_ele_area( c, tri_elements[i] );
+    // cout << "   area = " << area << endl;
+    if ( area < area_eps ) {
+	degenerate = true;
+    }
+
+    // cout << "  " << p1 << endl;
+    // cout << "  " << p2 << endl;
+    // cout << "  " << p3 << endl;
+    if ( fabs(p1.x() - p2.x()) < FG_EPSILON &&
+	 fabs(p1.x() - p3.x()) < FG_EPSILON ) {
+	degenerate = true;
+    }
+    if ( fabs(p1.y() - p2.y()) < FG_EPSILON &&
+	 fabs(p1.y() - p3.y()) < FG_EPSILON ) {
+	degenerate = true;
+    }
+    if ( fabs(p1.z() - p2.z()) < FG_EPSILON &&
+	 fabs(p1.z() - p3.z()) < FG_EPSILON ) {
+	degenerate = true;
+    }
+
+    if ( degenerate ) {
+	sgSetVec3( normal, p1.x(), p1.y(), p1.z() );
+	sgNormalizeVec3( normal );
+	cout << "Degenerate tri!" << endl;
+    } else {
+	v1[0] = p2.x() - p1.x();
+	v1[1] = p2.y() - p1.y();
+	v1[2] = p2.z() - p1.z();
+	v2[0] = p3.x() - p1.x();
+	v2[1] = p3.y() - p1.y();
+	v2[2] = p3.z() - p1.z();
+
+	sgVectorProductVec3( normal, v1, v2 );
+	sgNormalizeVec3( normal );
+    }
 
     return Point3D( normal[0], normal[1], normal[2] );
 }
@@ -574,8 +636,9 @@ static point_list gen_face_normals( FGConstruct& c ) {
 
     triele_list tri_elements = c.get_tri_elements();
     for ( int i = 0; i < (int)tri_elements.size(); i++ ) {
-	// cout << calc_normal( i ) << endl;
+	Point3D p = calc_normal(c,  i );
-	face_normals.push_back( calc_normal( c, i ) );
+	cout << p << endl;
+	face_normals.push_back( p );
     }
 
     return face_normals;
@@ -587,30 +650,32 @@ static point_list gen_point_normals( FGConstruct& c ) {
     point_list point_normals;
 
     Point3D normal;
-    cout << "caculating node normals" << endl;
+    cout << "calculating node normals" << endl;
 
     point_list wgs84_nodes = c.get_wgs84_nodes();
     belongs_to_list reverse_ele_lookup = c.get_reverse_ele_lookup();
     point_list face_normals = c.get_face_normals();
+    triele_list tri_elements = c.get_tri_elements();
 
     // for each node
     for ( int i = 0; i < (int)wgs84_nodes.size(); ++i ) {
 	int_list tri_list = reverse_ele_lookup[i];
-
+	double total_area = 0.0;
-	int_list_iterator current = tri_list.begin();
-	int_list_iterator last = tri_list.end();
 
 	Point3D average( 0.0 );
 
 	// for each triangle that shares this node
-	for ( ; current != last; ++current ) {
+	for ( int j = 0; j < (int)tri_list.size(); ++j ) {
-	    normal = face_normals[ *current ];
+	    normal = face_normals[ tri_list[j] ];
+	    double area = tri_ele_area( c, tri_elements[ tri_list[j] ] );
+	    normal *= area;	// scale normal weight relative to area
+	    total_area += area;
 	    average += normal;
 	    // cout << normal << endl;
 	}
 
-	average /= tri_list.size();
+	average /= total_area;
-	// cout << "average = " << average << endl;
+	cout << "average = " << average << endl;
 
 	point_normals.push_back( average );
     }
@@ -801,6 +866,22 @@ static void construct_tile( FGConstruct& c ) {
     c.set_tri_elements( t.get_elelist() );
     c.set_tri_segs( t.get_out_segs() );
 
+    // double check on the off chance that the triangulator was forced
+    // to introduce additional points
+    if ( c.get_tri_nodes().size() > c.get_point_normals().size() ) {
+	cout << "oops, need to add normals :-(" << endl;
+	point_list normals = c.get_point_normals();
+	int start = normals.size();
+	int end = c.get_tri_nodes().size();
+	for ( int i = start; i < end; ++i ) {
+	    cout << "adding a normal for " << c.get_tri_nodes().get_node(i)
+		 << endl;
+	    Point3D p = tgFakeNormal( c.get_tri_nodes().get_node(i) );
+	    normals.push_back( p );
+	}
+	c.set_point_normals( normals );
+    }
+
     // calculate wgs84 (cartesian) form of node list
     build_wgs_84_point_list( c, array );
 

src/BuildTiles/Match/match.cxx

@@ -53,13 +53,15 @@ void FGMatch::scan_share_file( const string& dir, const FGBucket& b,
 {
     string file = dir + "/"  + b.gen_base_path() + "/" + b.gen_index_str();
 
+    cout << "reading shared data from " << file << endl;
+
     fg_gzifstream in( file );
     if ( !in.is_open() ) {
         cout << "Cannot open file: " << file << endl;
 	return;
     }
 
-    cout << "reading shared data from " << file << endl;
+    cout << "open successful." << endl;
 
     string target;
     if ( search == SW_Corner ) {
@@ -268,7 +270,7 @@ void FGMatch::load_neighbor_shared( FGConstruct& c ) {
 
 
 // fake a normal for a point which is basically straight up
-static Point3D fake_normal( const Point3D& p ) {
+Point3D tgFakeNormal( const Point3D& p ) {
     Point3D radians = Point3D( p.x() * DEG_TO_RAD,
 			       p.y() * DEG_TO_RAD,
 			       p.z() );
@@ -276,7 +278,7 @@ static Point3D fake_normal( const Point3D& p ) {
     double len = Point3D(0.0).distance3D(cart);
     // cout << "len = " << len << endl;
     cart /= len;
-    // cout << "fake normal = " << cart << endl;
+    cout << "new fake normal = " << cart << endl;
 
     return cart;
 }
@@ -303,19 +305,19 @@ void FGMatch::split_tile( FGConstruct& c ) {
     // defaults "just in case"
     if ( ! sw_flag ) {
 	sw_node = Point3D( min.x, min.y, 0.0 );
-	sw_normal = fake_normal( sw_node );
+	sw_normal = tgFakeNormal( sw_node );
     }
     if ( ! se_flag ) {
 	se_node = Point3D( max.x, min.y, 0.0 );
-	se_normal = fake_normal( se_node );
+	se_normal = tgFakeNormal( se_node );
     }
     if ( ! nw_flag ) {
 	nw_node = Point3D( min.x, max.y, 0.0 );
- 	nw_normal = fake_normal( nw_node );
+ 	nw_normal = tgFakeNormal( nw_node );
     }
     if ( ! ne_flag ) {
 	ne_node = Point3D( max.x, max.y, 0.0 );
- 	ne_normal = fake_normal( ne_node );
+ 	ne_normal = tgFakeNormal( ne_node );
    }
 
     // separate nodes and normals into components
@@ -738,7 +740,7 @@ void FGMatch::assemble_tile( FGConstruct& c ) {
 	if ( n1 >= (int)new_normals.size() ) {
 	    cout << "Adding a segment resulted in a new node, faking a normal" 
 		 << endl;
-	    Point3D fake = fake_normal( p1 );
+	    Point3D fake = tgFakeNormal( p1 );
 	    insert_normal( new_normals, fake, n1 );
 	}
 
@@ -746,7 +748,7 @@ void FGMatch::assemble_tile( FGConstruct& c ) {
 	if ( n2 >= (int)new_normals.size() ) {
 	    cout << "Adding a segment resulted in a new node, faking a normal" 
 		 << endl;
-	    Point3D fake = fake_normal( p2 );
+	    Point3D fake = tgFakeNormal( p2 );
 	    insert_normal( new_normals, fake, n2 );
 	}
 

src/BuildTiles/Match/match.hxx

@@ -106,4 +106,8 @@ public:
 };
 
 
+// fake a normal for a point which is basically straight up
+Point3D tgFakeNormal( const Point3D& p );
+
+
 #endif // _MATCH_HXX