From c5faeae3570063fef818e0f6fddae5da3d5be29f Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Sun, 8 Feb 2015 06:52:39 -0500 Subject: [PATCH] some fixes for interpolating added points to the airport holes. fixes most gaps when allowing roads on airport landclass add zero area triangle check when triangulating without extra points. --- src/BuildTiles/Main/tgconstruct_cleanup.cxx | 4 +- src/Lib/terragear/tg_contour.cxx | 144 ++++++++++++++++++++ src/Lib/terragear/tg_contour.hxx | 3 + src/Lib/terragear/tg_nodes.cxx | 45 +++++- src/Lib/terragear/tg_nodes.hxx | 32 +++-- src/Lib/terragear/tg_polygon.cxx | 17 +++ src/Lib/terragear/tg_polygon.hxx | 3 + src/Lib/terragear/tg_polygon_clean.cxx | 12 +- src/Lib/terragear/tg_polygon_clip.cxx | 12 +- src/Lib/terragear/tg_polygon_tesselate.cxx | 7 +- 10 files changed, 255 insertions(+), 24 deletions(-) diff --git a/src/BuildTiles/Main/tgconstruct_cleanup.cxx b/src/BuildTiles/Main/tgconstruct_cleanup.cxx index 1d4b3f92..255fc35b 100644 --- a/src/BuildTiles/Main/tgconstruct_cleanup.cxx +++ b/src/BuildTiles/Main/tgconstruct_cleanup.cxx @@ -34,7 +34,7 @@ void TGConstruct::FixTJunctions( void ) { int before, after; - std::vector points; + std::vector points; tgRectangle bb; // traverse each poly, and add intermediate nodes @@ -42,7 +42,7 @@ void TGConstruct::FixTJunctions( void ) { for( unsigned int j = 0; j < polys_clipped.area_size(i); ++j ) { tgPolygon current = polys_clipped.get_poly(i, j); bb = current.GetBoundingBox(); - nodes.get_geod_inside( bb.getMin(), bb.getMax(), points ); + nodes.get_nodes_inside( bb.getMin(), bb.getMax(), points ); before = current.TotalNodes(); current = tgPolygon::AddColinearNodes( current, points ); diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx index 9898a12d..2a006078 100644 --- a/src/Lib/terragear/tg_contour.cxx +++ b/src/Lib/terragear/tg_contour.cxx @@ -6,6 +6,7 @@ #include "tg_accumulator.hxx" #include "tg_contour.hxx" #include "tg_polygon.hxx" +#include "tg_unique_tgnode.hxx" tgContour tgContour::Snap( const tgContour& subject, double snap ) { @@ -621,6 +622,87 @@ static bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, return found_node; } +static bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, + const std::vector& nodes, TGNode*& result, + double bbEpsilon, double errEpsilon ) +{ + bool found_node = false; + double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min; + + SGGeod p0 = start; + SGGeod p1 = end; + + double xdist = fabs(p0.getLongitudeDeg() - p1.getLongitudeDeg()); + double ydist = fabs(p0.getLatitudeDeg() - p1.getLatitudeDeg()); + + x_err_min = xdist + 1.0; + y_err_min = ydist + 1.0; + + if ( xdist > ydist ) { + // sort these in a sensible order + SGGeod p_min, p_max; + if ( p0.getLongitudeDeg() < p1.getLongitudeDeg() ) { + p_min = p0; + p_max = p1; + } else { + p_min = p1; + p_max = p0; + } + + m = (p_min.getLatitudeDeg() - p_max.getLatitudeDeg()) / (p_min.getLongitudeDeg() - p_max.getLongitudeDeg()); + b = p_max.getLatitudeDeg() - m * p_max.getLongitudeDeg(); + + for ( int i = 0; i < (int)nodes.size(); ++i ) { + // cout << i << endl; + SGGeod current = nodes[i]->GetPosition(); + + if ( (current.getLongitudeDeg() > (p_min.getLongitudeDeg() + (bbEpsilon))) && (current.getLongitudeDeg() < (p_max.getLongitudeDeg() - (bbEpsilon))) ) { + y_err = fabs(current.getLatitudeDeg() - (m * current.getLongitudeDeg() + b)); + + if ( y_err < errEpsilon ) { + found_node = true; + if ( y_err < y_err_min ) { + result = nodes[i]; + y_err_min = y_err; + } + } + } + } + } else { + // sort these in a sensible order + SGGeod p_min, p_max; + if ( p0.getLatitudeDeg() < p1.getLatitudeDeg() ) { + p_min = p0; + p_max = p1; + } else { + p_min = p1; + p_max = p0; + } + + m1 = (p_min.getLongitudeDeg() - p_max.getLongitudeDeg()) / (p_min.getLatitudeDeg() - p_max.getLatitudeDeg()); + b1 = p_max.getLongitudeDeg() - m1 * p_max.getLatitudeDeg(); + + for ( int i = 0; i < (int)nodes.size(); ++i ) { + SGGeod current = nodes[i]->GetPosition(); + + if ( (current.getLatitudeDeg() > (p_min.getLatitudeDeg() + (bbEpsilon))) && (current.getLatitudeDeg() < (p_max.getLatitudeDeg() - (bbEpsilon))) ) { + + x_err = fabs(current.getLongitudeDeg() - (m1 * current.getLatitudeDeg() + b1)); + + if ( x_err < errEpsilon ) { + found_node = true; + if ( x_err < x_err_min ) { + result = nodes[i]; + x_err_min = x_err; + } + } + } + } + } + + return found_node; +} + static void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, std::vector& nodes, tgContour& result, double bbEpsilon, double errEpsilon ) { SGGeod new_pt; @@ -639,6 +721,37 @@ static void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, std::vecto } } +extern SGGeod InterpolateElevation( const SGGeod& dst_node, const SGGeod& start, const SGGeod& end ); + +static void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, bool preserve3d, std::vector& nodes, tgContour& result, double bbEpsilon, double errEpsilon ) +{ + TGNode* new_pt; + SGGeod new_geode; + + SG_LOG(SG_GENERAL, SG_BULK, " " << p0 << " <==> " << p1 ); + + bool found_extra = FindIntermediateNode( p0, p1, nodes, new_pt, bbEpsilon, errEpsilon ); + + if ( found_extra ) { + if ( preserve3d ) { + // interpolate the new nodes elevation based on p0, p1 + new_geode = InterpolateElevation( new_pt->GetPosition(), p0, p1 ); + + SG_LOG(SG_GENERAL, SG_ALERT, "INTERPOLATE ELVATION between " << p0 << " and " << p1 << " returned elvation " << new_geode.getElevationM() ); + + new_pt->SetElevation( new_geode.getElevationM() ); + new_pt->SetFixedPosition( true ); + } + + AddIntermediateNodes( p0, new_pt->GetPosition(), preserve3d, nodes, result, bbEpsilon, errEpsilon ); + + result.AddNode( new_pt->GetPosition() ); + SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt->GetPosition() ); + + AddIntermediateNodes( new_pt->GetPosition(), p1, preserve3d, nodes, result, bbEpsilon, errEpsilon ); + } +} + tgContour tgContour::AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ) { SGGeod p0, p1; @@ -702,6 +815,37 @@ tgContour tgContour::AddColinearNodes( const tgContour& subject, std::vector& nodes ) +{ + SGGeod p0, p1; + tgContour result; + + for ( unsigned int n = 0; n < subject.GetSize()-1; n++ ) { + p0 = subject.GetNode( n ); + p1 = subject.GetNode( n+1 ); + + // add start of segment + result.AddNode( p0 ); + + // add intermediate points + AddIntermediateNodes( p0, p1, preserve3d, nodes, result, SG_EPSILON*10, SG_EPSILON*4 ); + } + + p0 = subject.GetNode( subject.GetSize() - 1 ); + p1 = subject.GetNode( 0 ); + + // add start of segment + result.AddNode( p0 ); + + // add intermediate points + AddIntermediateNodes( p0, p1, preserve3d, nodes, result, SG_EPSILON*10, SG_EPSILON*4 ); + + // maintain original hole flag setting + result.SetHole( subject.GetHole() ); + + return result; +} + // this is the opposite of FindColinearNodes - it takes a single SGGeode, // and tries to find the line segment the point is colinear with bool tgContour::FindColinearLine( const tgContour& subject, const SGGeod& node, SGGeod& start, SGGeod& end ) diff --git a/src/Lib/terragear/tg_contour.hxx b/src/Lib/terragear/tg_contour.hxx index ef286b39..def48cf4 100644 --- a/src/Lib/terragear/tg_contour.hxx +++ b/src/Lib/terragear/tg_contour.hxx @@ -14,6 +14,8 @@ #include "clipper.hpp" /* forward declarations */ +class TGNode; + class tgPolygon; typedef std::vector tgpolygon_list; @@ -112,6 +114,7 @@ public: static bool IsInside( const tgContour& inside, const tgContour& outside ); static tgContour AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ); static tgContour AddColinearNodes( const tgContour& subject, std::vector& nodes ); + static tgContour AddColinearNodes( const tgContour& subject, bool preserve3d, std::vector& nodes ); static bool FindColinearLine( const tgContour& subject, const SGGeod& node, SGGeod& start, SGGeod& end ); // conversions diff --git a/src/Lib/terragear/tg_nodes.cxx b/src/Lib/terragear/tg_nodes.cxx index 1280889f..ab66041e 100644 --- a/src/Lib/terragear/tg_nodes.cxx +++ b/src/Lib/terragear/tg_nodes.cxx @@ -45,6 +45,19 @@ bool TGNodes::get_geod_inside( const SGGeod& min, const SGGeod& max, std::vector return true; } +bool TGNodes::get_nodes_inside( const SGGeod& min, const SGGeod& max, std::vector& points ) const { + points.clear(); + for ( unsigned int i = 0; i < tg_node_list.size(); i++ ) { + SGGeod const& pt = tg_node_list[i].GetPosition(); + + if ( IsAlmostWithin( pt, min, max ) ) { + points.push_back( &tg_node_list[i] ); + } + } + + return true; +} + bool TGNodes::get_geod_edge( const SGBucket& b, std::vector& north, std::vector& south, std::vector& east, std::vector& west ) const { double north_compare = b.get_center_lat() + 0.5 * b.get_height(); double south_compare = b.get_center_lat() - 0.5 * b.get_height(); @@ -100,7 +113,7 @@ void TGNodes::init_spacial_query( void ) // generate the tuple Point pt( tg_node_list[i].GetPosition().getLongitudeDeg(), tg_node_list[i].GetPosition().getLatitudeDeg() ); double e( tg_node_list[i].GetPosition().getElevationM() ); - Point_and_Elevation pande(pt, e); + Point_and_Elevation pande( pt, e, &tg_node_list[i] ); // and insert into tree tg_kd_tree.insert( pande ); @@ -142,6 +155,36 @@ bool TGNodes::get_geod_inside( const SGGeod& min, const SGGeod& max, std::vector return true; } +bool TGNodes::get_nodes_inside( const SGGeod& min, const SGGeod& max, std::vector& points ) const { + points.clear(); + + // Have we generated the k-d tree? + if ( !kd_tree_valid ) { + SG_LOG(SG_GENERAL, SG_ALERT, "get_geod_inside called with invalid kdtree" ); + exit(0); + return false; + } + + // define an exact rectangulat range query (fuzziness=0) + Point ll( min.getLongitudeDeg() - fgPoint3_Epsilon, min.getLatitudeDeg() - fgPoint3_Epsilon ); + Point ur( max.getLongitudeDeg() + fgPoint3_Epsilon, max.getLatitudeDeg() + fgPoint3_Epsilon ); + Fuzzy_bb exact_bb(ll, ur); + + // list of tuples as a result + std::list result; + std::list::iterator it; + + // perform the query + tg_kd_tree.search(std::back_inserter( result ), exact_bb); + + // and convert the tuples back into SGGeod + for ( it = result.begin(); it != result.end(); it++ ) { + points.push_back( boost::get<2>(*it) ); + } + + return true; +} + // This query finds all nodes along the tile borders (north, south, east and west) bool TGNodes::get_geod_edge( const SGBucket& b, std::vector& north, std::vector& south, std::vector& east, std::vector& west ) const { double north_compare = b.get_center_lat() + 0.5 * b.get_height(); diff --git a/src/Lib/terragear/tg_nodes.hxx b/src/Lib/terragear/tg_nodes.hxx index a841b499..a4c9c6fd 100644 --- a/src/Lib/terragear/tg_nodes.hxx +++ b/src/Lib/terragear/tg_nodes.hxx @@ -15,17 +15,27 @@ #include /* Just use two dimensional lookup - elevation is a trait */ #include +#include +#include +#include + +#include "tg_unique_tgnode.hxx" + +#define FG_PROXIMITY_EPSILON 0.000001 +#define FG_COURSE_EPSILON 0.0001 + + typedef CGAL::Simple_cartesian Kernel; typedef Kernel::Point_2 Point; -typedef boost::tuple Point_and_Elevation; +typedef boost::tuple Point_and_Elevation; //definition of the property map #ifdef CGAL_OLD struct My_point_property_map{ - typedef Point value_type; - typedef const value_type& reference; - typedef const Point_and_Elevation& key_type; - typedef boost::readable_property_map_tag category; + typedef Point value_type; + typedef const value_type& reference; + typedef const Point_and_Elevation& key_type; + typedef boost::readable_property_map_tag category; }; #endif @@ -42,16 +52,6 @@ typedef CGAL::Fuzzy_iso_box Fuzzy_bb; typedef CGAL::Kd_tree Tree; -#include -#include -#include - -#include "tg_unique_tgnode.hxx" - -#define FG_PROXIMITY_EPSILON 0.000001 -#define FG_COURSE_EPSILON 0.0001 - - /* This class handles ALL of the nodes in a tile : 3d nodes in elevation data, 2d nodes generated from landclass, etc) */ class TGNodes { public: @@ -123,6 +123,8 @@ public: // Find all the nodes within a bounding box bool get_geod_inside( const SGGeod& min, const SGGeod& max, std::vector& points ) const; + bool get_nodes_inside( const SGGeod& min, const SGGeod& max, std::vector& points ) const; + // Find a;; the nodes on the tile edges bool get_geod_edge( const SGBucket& b, std::vector& north, std::vector& south, std::vector& east, std::vector& west ) const; diff --git a/src/Lib/terragear/tg_polygon.cxx b/src/Lib/terragear/tg_polygon.cxx index f5c8fd27..535f7ef3 100644 --- a/src/Lib/terragear/tg_polygon.cxx +++ b/src/Lib/terragear/tg_polygon.cxx @@ -121,6 +121,7 @@ tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, std::vector& nodes ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + result.AddContour( tgContour::AddColinearNodes( subject.GetContour(c), subject.GetPreserve3D(), nodes ) ); + } + + return result; +} + // this is the opposite of FindColinearNodes - it takes a single SGGeode, // and tries to find the line segment the point is colinear with bool tgPolygon::FindColinearLine( const tgPolygon& subject, SGGeod& node, SGGeod& start, SGGeod& end ) diff --git a/src/Lib/terragear/tg_polygon.hxx b/src/Lib/terragear/tg_polygon.hxx index 4793a9cc..f24f1948 100644 --- a/src/Lib/terragear/tg_polygon.hxx +++ b/src/Lib/terragear/tg_polygon.hxx @@ -86,6 +86,7 @@ void clipper_to_shapefile( ClipperLib::Paths polys, char* datasource ); // Forward Declaration: class tgPolygon; class tgChoppedPolygons; +class TGNode; typedef std::vector tgpolygon_list; typedef tgpolygon_list::iterator tgpolygon_list_iterator; @@ -413,6 +414,8 @@ public: static tgPolygon AddColinearNodes( const tgPolygon& subject, std::vector& nodes ); static bool FindColinearLine( const tgPolygon& subject, SGGeod& node, SGGeod& start, SGGeod& end ); + static tgPolygon AddColinearNodes( const tgPolygon& subject, std::vector& nodes ); + // IO void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); diff --git a/src/Lib/terragear/tg_polygon_clean.cxx b/src/Lib/terragear/tg_polygon_clean.cxx index ae2a3208..52fd0e51 100644 --- a/src/Lib/terragear/tg_polygon_clean.cxx +++ b/src/Lib/terragear/tg_polygon_clean.cxx @@ -9,6 +9,7 @@ tgPolygon tgPolygon::Snap( const tgPolygon& subject, double snap ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); for (unsigned int c = 0; c < subject.Contours(); c++) { result.AddContour( tgContour::Snap( subject.GetContour( c ), snap ) ); @@ -24,6 +25,7 @@ tgPolygon tgPolygon::RemoveDups( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); for ( unsigned int c = 0; c < subject.Contours(); c++ ) { result.AddContour( tgContour::RemoveDups( subject.GetContour( c ) ) ); @@ -39,6 +41,7 @@ tgPolygon tgPolygon::RemoveBadContours( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); for ( unsigned int c = 0; c < subject.Contours(); c++ ) { tgContour contour = subject.GetContour(c); @@ -59,7 +62,8 @@ tgPolygon tgPolygon::RemoveCycles( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); - + result.SetPreserve3D( subject.GetPreserve3D() ); + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { contours.clear(); @@ -84,6 +88,7 @@ tgPolygon tgPolygon::SplitLongEdges( const tgPolygon& subject, double dist ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); for ( unsigned c = 0; c < subject.Contours(); c++ ) { @@ -123,6 +128,7 @@ tgPolygon tgPolygon::StripHoles( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); return result; } @@ -148,6 +154,7 @@ tgPolygon tgPolygon::Simplify( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); return result; } @@ -160,6 +167,7 @@ tgPolygon tgPolygon::RemoveTinyContours( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); for ( unsigned int c = 0; c < subject.Contours(); c++ ) { tgContour contour = subject.GetContour( c ); @@ -183,6 +191,7 @@ tgPolygon tgPolygon::RemoveSpikes( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); for ( unsigned int c = 0; c < subject.Contours(); c++ ) { result.AddContour( tgContour::RemoveSpikes( subject.GetContour(c) ) ); @@ -275,6 +284,7 @@ tgcontour_list tgPolygon::MergeSlivers( tgpolygon_list& polys, tgcontour_list& s result.SetMaterial( polys[j].GetMaterial() ); result.SetTexParams( polys[j].GetTexParams() ); result.SetId( polys[j].GetId() ); + result.SetPreserve3D( polys[j].GetPreserve3D() ); polys[j] = result; done = true; } diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx index 1d68e8ce..1fbc958b 100644 --- a/src/Lib/terragear/tg_polygon_clip.cxx +++ b/src/Lib/terragear/tg_polygon_clip.cxx @@ -61,7 +61,9 @@ tgPolygon tgPolygon::Union( const tgPolygon& subject, tgPolygon& clip ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); + return result; } @@ -127,7 +129,9 @@ tgPolygon tgPolygon::Diff( const tgPolygon& subject, tgPolygon& clip ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); + return result; } @@ -164,7 +168,9 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + result.SetPreserve3D( subject.GetPreserve3D() ); + return result; } diff --git a/src/Lib/terragear/tg_polygon_tesselate.cxx b/src/Lib/terragear/tg_polygon_tesselate.cxx index 6d55d3e9..0f262d8f 100644 --- a/src/Lib/terragear/tg_polygon_tesselate.cxx +++ b/src/Lib/terragear/tg_polygon_tesselate.cxx @@ -190,8 +190,11 @@ void tgPolygon::Tesselate() SGGeod p1 = SGGeod::fromDeg( to_double(tri.vertex(1).x()), to_double(tri.vertex(1).y()) ); SGGeod p2 = SGGeod::fromDeg( to_double(tri.vertex(2).x()), to_double(tri.vertex(2).y()) ); - AddTriangle( p0, p1, p2 ); - + /* Check for Zero Area before inserting */ + if ( !SGGeod_isEqual2D( p0, p1 ) && !SGGeod_isEqual2D( p1, p2 ) && !SGGeod_isEqual2D( p0, p2 ) ) { + AddTriangle( p0, p1, p2 ); + } + ++count; } else { SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : face not in domain");