From 4f29c32c1862b4af9b74e40fcc4789112d019689 Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Mon, 25 Mar 2013 19:25:05 -0400 Subject: [PATCH] - fixes to get full world vmap0 working - use exact construction when triangulating in CGAL --- src/BuildTiles/Main/tgconstruct_cleanup.cxx | 23 +-- src/BuildTiles/Main/tgconstruct_clip.cxx | 25 ++-- src/BuildTiles/Main/tgconstruct_poly.cxx | 32 +++-- src/BuildTiles/Main/tgconstruct_tesselate.cxx | 6 +- src/Lib/terragear/tg_contour.cxx | 133 ++++++++---------- src/Lib/terragear/tg_contour.hxx | 7 +- src/Lib/terragear/tg_polygon.cxx | 9 +- src/Lib/terragear/tg_polygon_clean.cxx | 39 +++-- src/Lib/terragear/tg_polygon_tesselate.cxx | 57 ++++---- 9 files changed, 178 insertions(+), 153 deletions(-) diff --git a/src/BuildTiles/Main/tgconstruct_cleanup.cxx b/src/BuildTiles/Main/tgconstruct_cleanup.cxx index 68211c2c..1d4b3f92 100644 --- a/src/BuildTiles/Main/tgconstruct_cleanup.cxx +++ b/src/BuildTiles/Main/tgconstruct_cleanup.cxx @@ -85,26 +85,29 @@ void TGConstruct::CleanClippedPolys() { // Clean the polys for ( unsigned int area = 0; area < area_defs.size(); area++ ) { for( unsigned int p = 0; p < polys_clipped.area_size(area); p++ ) { - // step 1 : snap - tgPolygon poly = polys_clipped.get_poly(area, p); - poly = tgPolygon::Snap(poly, gSnap); + char layer[32]; + tgPolygon poly = polys_clipped.get_poly(area, p); + + // step 1 : snap + poly = tgPolygon::Snap(poly, gSnap); if ( IsDebugShape( poly.GetId() ) ) { - tgShapefile::FromPolygon( poly, ds_name, "snapped", "" ); + sprintf(layer, "snapped_%d", poly.GetId() ); + tgShapefile::FromPolygon( poly, ds_name, layer, "poly" ); } // step 2 : remove_dups poly = tgPolygon::RemoveDups( poly ); - if ( IsDebugShape( poly.GetId() ) ) { - tgShapefile::FromPolygon( poly, ds_name, "rem_dups", "" ); + sprintf(layer, "rem_dups_%d", poly.GetId() ); + tgShapefile::FromPolygon( poly, ds_name, layer, "poly" ); } - // step 3 : remove_bad_contours - poly = tgPolygon::RemoveBadContours( poly ); - + // step 3 : remove cycles + poly = tgPolygon::RemoveCycles( poly ); if ( IsDebugShape( poly.GetId() ) ) { - tgShapefile::FromPolygon( poly, ds_name, "rem_bcs", "" ); + sprintf(layer, "rem_cycles_%d", poly.GetId() ); + tgShapefile::FromPolygon( poly, ds_name, layer, "poly" ); } polys_clipped.set_poly(area, p, poly); diff --git a/src/BuildTiles/Main/tgconstruct_clip.cxx b/src/BuildTiles/Main/tgconstruct_clip.cxx index 6fc15d92..f4ddea0f 100644 --- a/src/BuildTiles/Main/tgconstruct_clip.cxx +++ b/src/BuildTiles/Main/tgconstruct_clip.cxx @@ -139,16 +139,6 @@ bool TGConstruct::ClipLandclassPolys( void ) { clipped = accum.Diff( tmp ); - if ( debug_area || debug_shape ) { - char layer[32]; - char name[32]; - - sprintf(layer, "post_clip_%d", polys_in.get_poly( i, j ).GetId() ); - sprintf(name, "shape %d,%d", i,j); - - tgShapefile::FromPolygon( clipped, ds_name, layer, name ); - } - // only add to output list if the clip left us with a polygon if ( clipped.Contours() > 0 ) { @@ -165,11 +155,15 @@ bool TGConstruct::ClipLandclassPolys( void ) { // shape.sps.push_back( sp ); polys_clipped.add_poly( i, clipped ); -#if 0 if ( debug_area || debug_shape ) { - WriteDebugShape( "clipped", shape ); + char layer[32]; + char name[32]; + + sprintf(layer, "post_clip_%d", polys_in.get_poly( i, j ).GetId() ); + sprintf(name, "shape %d,%d", i,j); + + tgShapefile::FromPolygon( clipped, ds_name, layer, name ); } -#endif } } @@ -208,6 +202,8 @@ bool TGConstruct::ClipLandclassPolys( void ) { // finally, what ever is left over goes to ocean remains = accum.Diff( safety_base ); + remains = tgPolygon::RemoveDups( remains ); + remains = tgPolygon::RemoveCycles( remains ); if ( remains.Contours() > 0 ) { // cout << "remains contours = " << remains.contours() << endl; @@ -236,6 +232,9 @@ bool TGConstruct::ClipLandclassPolys( void ) { if ( remains.Contours() > 0 ) { remains.SetMaterial( area_defs.get_sliver_area_name() ); remains.SetTexMethod( TG_TEX_BY_GEODE, bucket.get_center_lat() ); + remains.SetId(9999); + + SG_LOG( SG_CLIPPER, SG_INFO, "Adding remains to area " << area_defs.get_sliver_area_priority() ); polys_clipped.add_poly( area_defs.get_sliver_area_priority(), remains ); } } diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx index 89e180f7..bca12f5e 100644 --- a/src/BuildTiles/Main/tgconstruct_poly.cxx +++ b/src/BuildTiles/Main/tgconstruct_poly.cxx @@ -90,27 +90,29 @@ int TGConstruct::LoadLandclassPolys( void ) { poly.SetMaterial( material ); poly.SetId( cur_poly_id++ ); - polys_in.add_poly( area, poly ); - total_polys_read++; + if ( poly.Contours() ) { + polys_in.add_poly( area, poly ); + total_polys_read++; - // add the nodes - for (unsigned int j=0; j 2) { + if ( subject.GetArea() > SG_EPSILON*SG_EPSILON ) { + // Step 1 - find the first duplicate point + for ( unsigned int i = 0; i < subject.GetSize() && !split; i++ ) { + // first check until the end of the vector + for ( unsigned int j = i + 1; j < subject.GetSize() && !split; j++ ) { + if ( SGGeod_isEqual2D( subject.GetNode(i), subject.GetNode(j) ) ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "detected a dupe: i = " << i << " j = " << j ); + split = true; + + tgContour first; + if ( j < subject.GetSize()-1 ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "first contour is " << 0 << ".." << i << "," << j+1 << ".." << subject.GetSize()-1); + + // first contour is (0 .. i) + (j+1..size()-1) + for ( unsigned int n=0; n<=i; n++) { + first.AddNode( subject.GetNode(n) ); + } + for ( unsigned int n=j+1; n<=subject.GetSize()-1; n++) { + first.AddNode( subject.GetNode(n) ); + } + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, "first contour is " << 0 << ".." << i ); + + // first contour is (0 .. i) + for ( unsigned int n=0; n<=i; n++) { + first.AddNode( subject.GetNode(n) ); + } + } + + tgContour second; + SG_LOG(SG_GENERAL, SG_DEBUG, "second contour is " << i << ".." << j-1 ); + + // second contour is (i..j-1) + for ( unsigned int n=i; n 1) { - SG_LOG(SG_GENERAL, SG_DEBUG, "detected a small cycle: i = " << i << " j = " << j << " Erasing from 0 " << " to " << j-1 ); - result.RemoveNodeRange( 0, j-1 ); - } - - found = true; - } - } - } - } - - iters++; - SG_LOG(SG_GENERAL, SG_DEBUG, "remove small cycles : after " << iters << " iterations, contour has " << result.GetSize() << " points" ); - } while( found ); - - return result; + return split; } tgContour tgContour::RemoveDups( const tgContour& subject ) diff --git a/src/Lib/terragear/tg_contour.hxx b/src/Lib/terragear/tg_contour.hxx index 9026417f..2063f02c 100644 --- a/src/Lib/terragear/tg_contour.hxx +++ b/src/Lib/terragear/tg_contour.hxx @@ -17,6 +17,11 @@ class tgPolygon; typedef std::vector tgpolygon_list; +class tgContour; +typedef std::vector tgcontour_list; +typedef tgcontour_list::iterator tgcontour_list_iterator; +typedef tgcontour_list::const_iterator const_tgcontour_list_iterator; + class tgContour { public: @@ -74,9 +79,9 @@ public: static tgContour Snap( const tgContour& subject, double snap ); static tgContour RemoveDups( const tgContour& subject ); - static tgContour RemoveCycles( const tgContour& subject ); static tgContour SplitLongEdges( const tgContour& subject, double dist ); static tgContour RemoveSpikes( const tgContour& subject ); + static bool RemoveCycles( const tgContour& subject, tgcontour_list& result ); static tgPolygon Union( const tgContour& subject, tgPolygon& clip ); static tgPolygon Diff( const tgContour& subject, tgPolygon& clip ); diff --git a/src/Lib/terragear/tg_polygon.cxx b/src/Lib/terragear/tg_polygon.cxx index 20a8d2a2..9d3d9023 100644 --- a/src/Lib/terragear/tg_polygon.cxx +++ b/src/Lib/terragear/tg_polygon.cxx @@ -57,9 +57,7 @@ tgPolygon tgPolygon::Expand( const tgPolygon& subject, double offset ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); + result.SetId( subject.GetId() ); return result; } @@ -81,7 +79,7 @@ tgPolygon tgPolygon::Expand( const SGGeod& subject, double offset ) contour.AddNode( SGGeod::fromDeg( subject.getLongitudeDeg() + dlon, subject.getLatitudeDeg() + dlat ) ); contour.AddNode( SGGeod::fromDeg( subject.getLongitudeDeg() - dlon, subject.getLatitudeDeg() + dlat ) ); contour.SetHole(false); - + result.AddContour( contour ); return result; @@ -118,7 +116,8 @@ tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, std::vector= 3 ) { @@ -51,12 +54,24 @@ tgPolygon tgPolygon::RemoveBadContours( const tgPolygon& subject ) tgPolygon tgPolygon::RemoveCycles( const tgPolygon& subject ) { tgPolygon result; + tgcontour_list contours; result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); + result.SetId( subject.GetId() ); for ( unsigned int c = 0; c < subject.Contours(); c++ ) { - result.AddContour( tgContour::RemoveCycles( subject.GetContour( c ) ) ); + contours.clear(); + + // sglog().setLogLevels( SG_ALL, SG_DEBUG ); + SG_LOG(SG_GENERAL, SG_DEBUG, " REMOVE DUPES FOR CONTOUR : " << c << " with " << subject.ContourSize(c) << " nodes "); + tgContour::RemoveCycles( subject.GetContour( c ), contours ); + SG_LOG(SG_GENERAL, SG_DEBUG, " REMOVE DUPES FOR CONTOUR : " << c << " COMPLETE"); + //sglog().setLogLevels( SG_ALL, SG_INFO ); + + for ( unsigned int i=0; i< contours.size(); i++ ) { + result.AddContour( contours[i] ); + } } return result; @@ -68,7 +83,8 @@ tgPolygon tgPolygon::SplitLongEdges( const tgPolygon& subject, double dist ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + for ( unsigned c = 0; c < subject.Contours(); c++ ) { result.AddContour( tgContour::SplitLongEdges( subject.GetContour(c), dist ) ); @@ -106,7 +122,8 @@ tgPolygon tgPolygon::StripHoles( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + return result; } @@ -130,7 +147,8 @@ tgPolygon tgPolygon::Simplify( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + return result; } @@ -141,7 +159,8 @@ tgPolygon tgPolygon::RemoveTinyContours( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { tgContour contour = subject.GetContour( c ); double area = contour.GetArea(); @@ -163,7 +182,8 @@ tgPolygon tgPolygon::RemoveSpikes( const tgPolygon& subject ) result.SetMaterial( subject.GetMaterial() ); result.SetTexParams( subject.GetTexParams() ); - + result.SetId( subject.GetId() ); + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { result.AddContour( tgContour::RemoveSpikes( subject.GetContour(c) ) ); } @@ -254,6 +274,7 @@ tgcontour_list tgPolygon::MergeSlivers( tgpolygon_list& polys, tgcontour_list& s SG_LOG(SG_GENERAL, SG_DEBUG, " FOUND a poly to merge the sliver with"); result.SetMaterial( polys[j].GetMaterial() ); result.SetTexParams( polys[j].GetTexParams() ); + result.SetId( polys[j].GetId() ); polys[j] = result; done = true; } diff --git a/src/Lib/terragear/tg_polygon_tesselate.cxx b/src/Lib/terragear/tg_polygon_tesselate.cxx index 1d7e4b0f..503e3b29 100644 --- a/src/Lib/terragear/tg_polygon_tesselate.cxx +++ b/src/Lib/terragear/tg_polygon_tesselate.cxx @@ -1,8 +1,10 @@ #include +#include -#include +#include #include #include +#include #include #include @@ -21,14 +23,15 @@ struct FaceInfo2 } }; -typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Exact_predicates_exact_constructions_kernel K; typedef CGAL::Triangulation_vertex_base_2 Vb; typedef CGAL::Triangulation_face_base_with_info_2 Fbb; typedef CGAL::Constrained_triangulation_face_base_2 Fb; typedef CGAL::Triangulation_data_structure_2 TDS; -typedef CGAL::Exact_predicates_tag Itag; +typedef CGAL::Exact_intersections_tag Itag; typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; -typedef CDT::Point Point; +typedef CGAL::Constrained_triangulation_plus_2 CDTPlus; +typedef CDTPlus::Point Point; typedef CGAL::Polygon_2 Polygon_2; typedef CGAL::Triangle_2 Triangle_2; @@ -38,17 +41,17 @@ static void tg_mark_domains(CDT& ct, CDT::Face_handle start, int index, std::lis return; } - std::list queue; + std::list queue; queue.push_back(start); while( !queue.empty() ){ - CDT::Face_handle fh = queue.front(); + CDTPlus::Face_handle fh = queue.front(); queue.pop_front(); if(fh->info().nesting_level == -1) { fh->info().nesting_level = index; for(int i = 0; i < 3; i++) { - CDT::Edge e(fh,i); - CDT::Face_handle n = fh->neighbor(i); + CDTPlus::Edge e(fh,i); + CDTPlus::Face_handle n = fh->neighbor(i); if(n->info().nesting_level == -1) { if(ct.is_constrained(e)) border.push_back(e); else queue.push_back(n); @@ -64,32 +67,32 @@ static void tg_mark_domains(CDT& ct, CDT::Face_handle start, int index, std::lis //level of 0. Then we recursively consider the non-explored facets incident //to constrained edges bounding the former set and increase the nesting level by 1. //Facets in the domain are those with an odd nesting level. -static void tg_mark_domains(CDT& cdt) +static void tg_mark_domains(CDTPlus& cdt) { - for(CDT::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){ + for(CDTPlus::All_faces_iterator it = cdt.all_faces_begin(); it != cdt.all_faces_end(); ++it){ it->info().nesting_level = -1; } int index = 0; - std::list border; + std::list border; tg_mark_domains(cdt, cdt.infinite_face(), index++, border); while(! border.empty()) { - CDT::Edge e = border.front(); + CDTPlus::Edge e = border.front(); border.pop_front(); - CDT::Face_handle n = e.first->neighbor(e.second); + CDTPlus::Face_handle n = e.first->neighbor(e.second); if(n->info().nesting_level == -1) { tg_mark_domains(cdt, n, e.first->info().nesting_level+1, border); } } } -static void tg_insert_polygon(CDT& cdt,const Polygon_2& polygon) +static void tg_insert_polygon(CDTPlus& cdt,const Polygon_2& polygon) { if ( polygon.is_empty() ) return; - CDT::Vertex_handle v_prev=cdt.insert(*CGAL::cpp0x::prev(polygon.vertices_end())); + CDTPlus::Vertex_handle v_prev=cdt.insert(*CGAL::cpp0x::prev(polygon.vertices_end())); for (Polygon_2::Vertex_iterator vit=polygon.vertices_begin(); vit!=polygon.vertices_end();++vit) { - CDT::Vertex_handle vh=cdt.insert(*vit); + CDTPlus::Vertex_handle vh=cdt.insert(*vit); cdt.insert_constraint(vh,v_prev); v_prev=vh; } @@ -97,7 +100,7 @@ static void tg_insert_polygon(CDT& cdt,const Polygon_2& polygon) void tgPolygon::Tesselate( const std::vector& extra ) { - CDT cdt; + CDTPlus cdt; SG_LOG( SG_GENERAL, SG_DEBUG, "Tess with extra" ); @@ -125,16 +128,17 @@ void tgPolygon::Tesselate( const std::vector& extra ) tg_insert_polygon(cdt, poly); } + tg_mark_domains( cdt ); int count=0; - for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end(); ++fit) { + for (CDTPlus::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end(); ++fit) { if ( fit->info().in_domain() ) { Triangle_2 tri = cdt.triangle(fit); - SGGeod p0 = SGGeod::fromDeg( tri.vertex(0).x(), tri.vertex(0).y() ); - SGGeod p1 = SGGeod::fromDeg( tri.vertex(1).x(), tri.vertex(1).y() ); - SGGeod p2 = SGGeod::fromDeg( tri.vertex(2).x(), tri.vertex(2).y() ); + SGGeod p0 = SGGeod::fromDeg( to_double(tri.vertex(0).x()), to_double(tri.vertex(0).y()) ); + 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 ); @@ -146,7 +150,7 @@ void tgPolygon::Tesselate( const std::vector& extra ) void tgPolygon::Tesselate() { - CDT cdt; + CDTPlus cdt; SG_LOG( SG_GENERAL, SG_DEBUG, "Tess" ); @@ -166,18 +170,19 @@ void tgPolygon::Tesselate() tg_insert_polygon(cdt, poly); } + assert(cdt.is_valid()); tg_mark_domains( cdt ); int count=0; - for (CDT::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end(); ++fit) { + for (CDTPlus::Finite_faces_iterator fit=cdt.finite_faces_begin(); fit!=cdt.finite_faces_end(); ++fit) { if ( fit->info().in_domain() ) { SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : face in domain"); Triangle_2 tri = cdt.triangle(fit); - SGGeod p0 = SGGeod::fromDeg( tri.vertex(0).x(), tri.vertex(0).y() ); - SGGeod p1 = SGGeod::fromDeg( tri.vertex(1).x(), tri.vertex(1).y() ); - SGGeod p2 = SGGeod::fromDeg( tri.vertex(2).x(), tri.vertex(2).y() ); + SGGeod p0 = SGGeod::fromDeg( to_double(tri.vertex(0).x()), to_double(tri.vertex(0).y()) ); + 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 );