From 59a35fbe8a4cf86827b709310724dad9b7d9a4c0 Mon Sep 17 00:00:00 2001 From: curt Date: Wed, 22 Nov 2000 20:41:22 +0000 Subject: [PATCH] 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 | 53 ++++++------- src/BuildTiles/GenOutput/genobj.cxx | 2 + src/BuildTiles/Main/main.cxx | 111 ++++++++++++++++++++++++---- src/BuildTiles/Match/match.cxx | 20 ++--- src/BuildTiles/Match/match.hxx | 4 + 5 files changed, 140 insertions(+), 50 deletions(-) diff --git a/src/BuildTiles/Clipper/clipper.cxx b/src/BuildTiles/Clipper/clipper.cxx index 18925856..3d198e89 100644 --- a/src/BuildTiles/Clipper/clipper.cxx +++ b/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 " - << clipped.polys[area].size() << " polys" << endl; + // 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; + // 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() << ", " - << contour[k].y() << endl; + // cout << " " << contour[k].x() << ", " + // << 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 << " result_diff = " << result_diff.contours() << endl; - cout << " slivers = " << slivers.contours() << endl; + // cout << " After sliver move:" << endl; + // cout << " result_diff = " << result_diff.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 << " remains = " << remains.contours() << endl; - cout << " slivers = " << slivers.contours() << endl; + // cout << " After sliver move:" << endl; + // cout << " remains = " << remains.contours() << endl; + // cout << " slivers = " << slivers.contours() << endl; // merge any slivers with previously clipped // neighboring polygons diff --git a/src/BuildTiles/GenOutput/genobj.cxx b/src/BuildTiles/GenOutput/genobj.cxx index 15616073..9e7ca73e 100644 --- a/src/BuildTiles/GenOutput/genobj.cxx +++ b/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(); diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index 28062d04..6b990248 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -46,6 +46,7 @@ #include #include +#include #include #include #include @@ -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(); - v2[0] = p3.x() - p1.x(); v2[1] = p3.y() - p1.y(); v2[2] = p3.z() - p1.z(); + // do some sanity checking. With the introduction of landuse + // 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 ); - sgNormalizeVec3( normal ); + bool degenerate = false; + 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; - face_normals.push_back( calc_normal( c, i ) ); + Point3D p = 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]; - - int_list_iterator current = tri_list.begin(); - int_list_iterator last = tri_list.end(); + double total_area = 0.0; Point3D average( 0.0 ); // for each triangle that shares this node - for ( ; current != last; ++current ) { - normal = face_normals[ *current ]; + for ( int j = 0; j < (int)tri_list.size(); ++j ) { + 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(); - // cout << "average = " << average << endl; + average /= total_area; + 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 ); diff --git a/src/BuildTiles/Match/match.cxx b/src/BuildTiles/Match/match.cxx index 3b7693bb..fcb91478 100644 --- a/src/BuildTiles/Match/match.cxx +++ b/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 ); } diff --git a/src/BuildTiles/Match/match.hxx b/src/BuildTiles/Match/match.hxx index aeaac5d4..42242906 100644 --- a/src/BuildTiles/Match/match.hxx +++ b/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