From 10db5bfbffcb3cc4bdbf9dcb81bdbdde1fea5761 Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Sun, 16 Dec 2012 10:46:25 -0500 Subject: [PATCH] - renamed Polygon library to terragear - split the new implementation into respective class files (tg_misc is a catch all...) --- src/Airports/GenAirports850/CMakeLists.txt | 2 +- src/Airports/GenAirports850/airport.cxx | 17 +- src/Airports/GenAirports850/apt_math.hxx | 2 +- src/Airports/GenAirports850/closedpoly.cxx | 7 +- src/Airports/GenAirports850/closedpoly.hxx | 2 +- src/Airports/GenAirports850/elevations.hxx | 2 +- src/Airports/GenAirports850/helipad.hxx | 4 +- src/Airports/GenAirports850/linearfeature.hxx | 4 +- .../GenAirports850/linked_objects.hxx | 2 +- src/Airports/GenAirports850/main.cxx | 3 - src/Airports/GenAirports850/object.hxx | 3 +- src/Airports/GenAirports850/runway.cxx | 2 +- src/Airports/GenAirports850/runway.hxx | 4 +- src/Airports/GenAirports850/rwy_gen.cxx | 10 +- src/Airports/GenAirports850/scheduler.hxx | 2 +- src/Airports/GenAirports850/taxiway.cxx | 6 +- src/Airports/GenAirports850/taxiway.hxx | 4 +- src/BuildTiles/Main/CMakeLists.txt | 4 +- src/BuildTiles/Main/cliptst.cxx | 12 +- src/BuildTiles/Main/main.cxx | 3 - src/BuildTiles/Main/tgconstruct.hxx | 2 +- src/BuildTiles/Main/tgconstruct_cleanup.cxx | 8 +- src/BuildTiles/Main/tgconstruct_clip.cxx | 17 +- src/BuildTiles/Main/tgconstruct_output.cxx | 4 +- src/BuildTiles/Main/tgconstruct_tesselate.cxx | 4 +- src/BuildTiles/Main/tglandclass.hxx | 2 +- src/Lib/CMakeLists.txt | 2 +- src/Lib/Polygon/CMakeLists.txt | 21 - src/Lib/Polygon/polygon.cxx | 2966 ----------------- src/Lib/terragear/CMakeLists.txt | 35 + src/Lib/{Polygon => terragear}/TNT/jama_qr.h | 0 .../{Polygon => terragear}/TNT/tnt_array1d.h | 0 .../{Polygon => terragear}/TNT/tnt_array2d.h | 0 .../{Polygon => terragear}/TNT/tnt_i_refvec.h | 0 .../TNT/tnt_math_utils.h | 0 src/Lib/{Polygon => terragear}/clipper.cpp | 0 src/Lib/{Polygon => terragear}/clipper.hpp | 0 src/Lib/terragear/tg_accumulator.cxx | 131 + src/Lib/terragear/tg_accumulator.hxx | 25 + src/Lib/terragear/tg_chopper.cxx | 274 ++ src/Lib/terragear/tg_chopper.hxx | 27 + src/Lib/terragear/tg_contour.cxx | 855 +++++ src/Lib/terragear/tg_contour.hxx | 112 + src/Lib/terragear/tg_light.hxx | 95 + src/Lib/terragear/tg_misc.cxx | 269 ++ src/Lib/terragear/tg_misc.hxx | 33 + src/Lib/{Polygon => terragear}/tg_nodes.cxx | 0 src/Lib/{Polygon => terragear}/tg_nodes.hxx | 2 +- src/Lib/terragear/tg_polygon.cxx | 519 +++ .../polygon.hxx => terragear/tg_polygon.hxx} | 280 +- src/Lib/terragear/tg_polygon_clean.cxx | 272 ++ src/Lib/terragear/tg_polygon_clip.cxx | 175 + src/Lib/terragear/tg_polygon_tesselate.cxx | 194 ++ .../tg_rectangle.cxx} | 2 +- .../tg_rectangle.hxx} | 0 src/Lib/terragear/tg_shapefile.cxx | 208 ++ src/Lib/terragear/tg_shapefile.hxx | 20 + .../surface.cxx => terragear/tg_surface.cxx} | 2 +- .../surface.hxx => terragear/tg_surface.hxx} | 2 +- .../{Polygon => terragear}/tg_unique_geod.hxx | 0 .../tg_unique_tgnode.hxx | 0 .../tg_unique_vec2f.hxx | 0 .../tg_unique_vec3d.hxx | 0 .../tg_unique_vec3f.hxx | 0 src/Prep/OGRDecode/CMakeLists.txt | 2 +- src/Prep/OGRDecode/ogr-decode.cxx | 8 +- src/Utils/poly2ogr/CMakeLists.txt | 3 +- src/Utils/poly2ogr/poly2ogr.cxx | 2 +- 68 files changed, 3350 insertions(+), 3318 deletions(-) delete mode 100644 src/Lib/Polygon/CMakeLists.txt delete mode 100644 src/Lib/Polygon/polygon.cxx create mode 100644 src/Lib/terragear/CMakeLists.txt rename src/Lib/{Polygon => terragear}/TNT/jama_qr.h (100%) rename src/Lib/{Polygon => terragear}/TNT/tnt_array1d.h (100%) rename src/Lib/{Polygon => terragear}/TNT/tnt_array2d.h (100%) rename src/Lib/{Polygon => terragear}/TNT/tnt_i_refvec.h (100%) rename src/Lib/{Polygon => terragear}/TNT/tnt_math_utils.h (100%) rename src/Lib/{Polygon => terragear}/clipper.cpp (100%) rename src/Lib/{Polygon => terragear}/clipper.hpp (100%) create mode 100644 src/Lib/terragear/tg_accumulator.cxx create mode 100644 src/Lib/terragear/tg_accumulator.hxx create mode 100644 src/Lib/terragear/tg_chopper.cxx create mode 100644 src/Lib/terragear/tg_chopper.hxx create mode 100644 src/Lib/terragear/tg_contour.cxx create mode 100644 src/Lib/terragear/tg_contour.hxx create mode 100644 src/Lib/terragear/tg_light.hxx create mode 100644 src/Lib/terragear/tg_misc.cxx create mode 100644 src/Lib/terragear/tg_misc.hxx rename src/Lib/{Polygon => terragear}/tg_nodes.cxx (100%) rename src/Lib/{Polygon => terragear}/tg_nodes.hxx (99%) create mode 100644 src/Lib/terragear/tg_polygon.cxx rename src/Lib/{Polygon/polygon.hxx => terragear/tg_polygon.hxx} (64%) create mode 100644 src/Lib/terragear/tg_polygon_clean.cxx create mode 100644 src/Lib/terragear/tg_polygon_clip.cxx create mode 100644 src/Lib/terragear/tg_polygon_tesselate.cxx rename src/Lib/{Polygon/rectangle.cxx => terragear/tg_rectangle.cxx} (98%) rename src/Lib/{Polygon/rectangle.hxx => terragear/tg_rectangle.hxx} (100%) create mode 100644 src/Lib/terragear/tg_shapefile.cxx create mode 100644 src/Lib/terragear/tg_shapefile.hxx rename src/Lib/{Polygon/surface.cxx => terragear/tg_surface.cxx} (99%) rename src/Lib/{Polygon/surface.hxx => terragear/tg_surface.hxx} (99%) rename src/Lib/{Polygon => terragear}/tg_unique_geod.hxx (100%) rename src/Lib/{Polygon => terragear}/tg_unique_tgnode.hxx (100%) rename src/Lib/{Polygon => terragear}/tg_unique_vec2f.hxx (100%) rename src/Lib/{Polygon => terragear}/tg_unique_vec3d.hxx (100%) rename src/Lib/{Polygon => terragear}/tg_unique_vec3f.hxx (100%) diff --git a/src/Airports/GenAirports850/CMakeLists.txt b/src/Airports/GenAirports850/CMakeLists.txt index fd02c58c..3c2ac3da 100644 --- a/src/Airports/GenAirports850/CMakeLists.txt +++ b/src/Airports/GenAirports850/CMakeLists.txt @@ -25,7 +25,7 @@ add_executable(genapts850 ) target_link_libraries(genapts850 - Polygon + terragear Array ${GDAL_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} diff --git a/src/Airports/GenAirports850/airport.cxx b/src/Airports/GenAirports850/airport.cxx index dcfaa1a2..6b9099a4 100644 --- a/src/Airports/GenAirports850/airport.cxx +++ b/src/Airports/GenAirports850/airport.cxx @@ -12,10 +12,11 @@ #include #include -#include -#include -#include -#include +#include +#include +#include +#include +#include #include "airport.hxx" #include "beznode.hxx" @@ -627,7 +628,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src ) for ( unsigned int k = 0; k < rwy_polys.size(); ++k ) { tgPolygon poly = rwy_polys[k]; - poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes.get_list() ); + poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes ); GENAPT_LOG(SG_GENERAL, SG_DEBUG, "total size after add nodes = " << poly.TotalNodes()); rwy_polys[k] = poly; } @@ -636,7 +637,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src ) for ( unsigned int k = 0; k < pvmt_polys.size(); ++k ) { tgPolygon poly = pvmt_polys[k]; - poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes.get_list() ); + poly = tgPolygon::AddColinearNodes( poly, tmp_pvmt_nodes ); GENAPT_LOG(SG_GENERAL, SG_DEBUG, "total size after add nodes = " << poly.TotalNodes()); pvmt_polys[k] = poly; } @@ -645,7 +646,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src ) for ( unsigned int k = 0; k < line_polys.size(); ++k ) { tgPolygon poly = line_polys[k]; - poly = tgPolygon::AddColinearNodes( poly, tmp_feat_nodes.get_list() ); + poly = tgPolygon::AddColinearNodes( poly, tmp_feat_nodes ); GENAPT_LOG(SG_GENERAL, SG_DEBUG, "total size after add nodes = " << poly.TotalNodes()); line_polys[k] = poly; } @@ -673,7 +674,7 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src ) log_time = time(0); GENAPT_LOG( SG_GENERAL, SG_ALERT, "Finished cleaning polys for " << icao << " at " << DebugTimeToString(log_time) ); - base_poly = tgPolygon::AddColinearNodes( base_poly, tmp_pvmt_nodes.get_list() ); + base_poly = tgPolygon::AddColinearNodes( base_poly, tmp_pvmt_nodes ); base_poly = tgPolygon::Snap( base_poly, gSnap ); // Finally find slivers in base diff --git a/src/Airports/GenAirports850/apt_math.hxx b/src/Airports/GenAirports850/apt_math.hxx index e1c27ed3..d443705c 100644 --- a/src/Airports/GenAirports850/apt_math.hxx +++ b/src/Airports/GenAirports850/apt_math.hxx @@ -1,7 +1,7 @@ #ifndef _APT_MATH_HXX_ #define _APT_MATH_HXX_ -#include +#include tgContour gen_wgs84_area( SGGeod origin, double length_m, diff --git a/src/Airports/GenAirports850/closedpoly.cxx b/src/Airports/GenAirports850/closedpoly.cxx index 622fb3ea..f4504b59 100644 --- a/src/Airports/GenAirports850/closedpoly.cxx +++ b/src/Airports/GenAirports850/closedpoly.cxx @@ -2,7 +2,8 @@ #include -#include +#include +#include #include "global.hxx" #include "beznode.hxx" @@ -455,7 +456,7 @@ int ClosedPoly::BuildBtg( tgpolygon_list& rwy_polys, tgcontour_list& slivers, tg if ( is_pavement && pre_tess.Contours() ) { if( shapefile_name.size() ) { - tgPolygon::ToShapefile( pre_tess, "./airport_dbg", std::string("preclip"), shapefile_name ); + tgShapefile::FromPolygon( pre_tess, "./airport_dbg", std::string("preclip"), shapefile_name ); accum.ToShapefiles( "./airport_dbg", "accum" ); } @@ -463,7 +464,7 @@ int ClosedPoly::BuildBtg( tgpolygon_list& rwy_polys, tgcontour_list& slivers, tg if ( clipped.Contours() ) { if( shapefile_name.size() ) { - tgPolygon::ToShapefile( clipped, "./airport_dbg", std::string("postclip"), shapefile_name ); + tgShapefile::FromPolygon( clipped, "./airport_dbg", std::string("postclip"), shapefile_name ); } tgPolygon::RemoveSlivers( clipped, slivers ); diff --git a/src/Airports/GenAirports850/closedpoly.hxx b/src/Airports/GenAirports850/closedpoly.hxx index a61884e0..b31f2f09 100644 --- a/src/Airports/GenAirports850/closedpoly.hxx +++ b/src/Airports/GenAirports850/closedpoly.hxx @@ -1,7 +1,7 @@ #ifndef _BEZPOLY_H_ #define _BEZPOLY_H_ -#include +#include #include "beznode.hxx" #include "linearfeature.hxx" diff --git a/src/Airports/GenAirports850/elevations.hxx b/src/Airports/GenAirports850/elevations.hxx index e992a218..1da27119 100644 --- a/src/Airports/GenAirports850/elevations.hxx +++ b/src/Airports/GenAirports850/elevations.hxx @@ -20,7 +20,7 @@ // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // -#include +#include // lookup node elevations for each point in the SGGeod list. Returns // average of all points. Doesn't modify the original list. diff --git a/src/Airports/GenAirports850/helipad.hxx b/src/Airports/GenAirports850/helipad.hxx index 22e0976a..654734d4 100644 --- a/src/Airports/GenAirports850/helipad.hxx +++ b/src/Airports/GenAirports850/helipad.hxx @@ -16,7 +16,9 @@ #ifndef _HELIPAD_HXX #define _HELIPAD_HXX -#include +#include +#include +#include class Helipad { diff --git a/src/Airports/GenAirports850/linearfeature.hxx b/src/Airports/GenAirports850/linearfeature.hxx index 61131597..34675acc 100644 --- a/src/Airports/GenAirports850/linearfeature.hxx +++ b/src/Airports/GenAirports850/linearfeature.hxx @@ -1,7 +1,9 @@ #ifndef _LINEARFEATURE_H_ #define _LINEARFEATURE_H_ -#include +#include +#include +#include #include "beznode.hxx" diff --git a/src/Airports/GenAirports850/linked_objects.hxx b/src/Airports/GenAirports850/linked_objects.hxx index 68d099d5..c7ddf071 100644 --- a/src/Airports/GenAirports850/linked_objects.hxx +++ b/src/Airports/GenAirports850/linked_objects.hxx @@ -1,7 +1,7 @@ #ifndef _LINKED_OBJECTS_H_ #define _LINKED_OBJECTS_H_ -#include +#include class Windsock { diff --git a/src/Airports/GenAirports850/main.cxx b/src/Airports/GenAirports850/main.cxx index 75103312..a2238743 100644 --- a/src/Airports/GenAirports850/main.cxx +++ b/src/Airports/GenAirports850/main.cxx @@ -301,9 +301,6 @@ int main(int argc, char **argv) exit(-1); } - // Initialize shapefile support (for debugging) - tgShapefile::Init(); - sg_gzifstream in( input_file ); if ( !in.is_open() ) { diff --git a/src/Airports/GenAirports850/object.hxx b/src/Airports/GenAirports850/object.hxx index f2ec27d6..361bd1b9 100644 --- a/src/Airports/GenAirports850/object.hxx +++ b/src/Airports/GenAirports850/object.hxx @@ -1,7 +1,8 @@ #ifndef _OBJECT_H_ #define _OBJECT_H_ -#include +#include +#include class LightingObj { diff --git a/src/Airports/GenAirports850/runway.cxx b/src/Airports/GenAirports850/runway.cxx index ef4a339b..6987f3ec 100644 --- a/src/Airports/GenAirports850/runway.cxx +++ b/src/Airports/GenAirports850/runway.cxx @@ -5,7 +5,7 @@ #include #include -#include +#include #include "global.hxx" #include "apt_math.hxx" diff --git a/src/Airports/GenAirports850/runway.hxx b/src/Airports/GenAirports850/runway.hxx index 60fea423..05424505 100644 --- a/src/Airports/GenAirports850/runway.hxx +++ b/src/Airports/GenAirports850/runway.hxx @@ -1,7 +1,9 @@ #ifndef _RUNWAY_H_ #define _RUNWAY_H_ -#include +#include +#include +#include #include "apt_math.hxx" diff --git a/src/Airports/GenAirports850/rwy_gen.cxx b/src/Airports/GenAirports850/rwy_gen.cxx index ee5ffe0e..89032a2a 100644 --- a/src/Airports/GenAirports850/rwy_gen.cxx +++ b/src/Airports/GenAirports850/rwy_gen.cxx @@ -19,6 +19,8 @@ // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // +#include + #include #include #include @@ -26,12 +28,12 @@ #include +#include + #include "global.hxx" #include "beznode.hxx" #include "runway.hxx" -#include - using std::string; struct sections { @@ -301,14 +303,14 @@ void Runway::gen_runway_section( const tgPolygon& runway, GENAPT_LOG(SG_GENERAL, SG_DEBUG, section ); if( shapefile_name.size() ) { - tgPolygon::ToShapefile( section, "./airport_dbg", std::string("preclip"), shapefile_name ); + tgShapefile::FromPolygon( section, "./airport_dbg", std::string("preclip"), shapefile_name ); } // Clip the new polygon against what ever has already been created. tgPolygon clipped = accum.Diff( section ); if( shapefile_name.size() ) { - tgPolygon::ToShapefile( clipped, "./airport_dbg", std::string("postclip"), shapefile_name ); + tgShapefile::FromPolygon( clipped, "./airport_dbg", std::string("postclip"), shapefile_name ); } tgPolygon::RemoveSlivers( clipped, slivers ); diff --git a/src/Airports/GenAirports850/scheduler.hxx b/src/Airports/GenAirports850/scheduler.hxx index c68917dd..7a7870b2 100644 --- a/src/Airports/GenAirports850/scheduler.hxx +++ b/src/Airports/GenAirports850/scheduler.hxx @@ -10,7 +10,7 @@ #include #include #include -#include +#include #include "airport.hxx" #define P_STATE_INIT (0) diff --git a/src/Airports/GenAirports850/taxiway.cxx b/src/Airports/GenAirports850/taxiway.cxx index 88b2a925..6ea1b860 100644 --- a/src/Airports/GenAirports850/taxiway.cxx +++ b/src/Airports/GenAirports850/taxiway.cxx @@ -4,6 +4,8 @@ #include #include +#include + #include "global.hxx" #include "apt_math.hxx" #include "beznode.hxx" @@ -135,7 +137,7 @@ int Taxiway::BuildBtg( tgpolygon_list& rwy_polys, tglightcontour_list& rwy_light tgPolygon taxi_poly; taxi_poly.AddContour( taxi_contour ); - tgPolygon::ToShapefile( taxi_poly, "./airport_dbg", std::string("preclip"), shapefile_name ); + tgShapefile::FromPolygon( taxi_poly, "./airport_dbg", std::string("preclip"), shapefile_name ); accum.ToShapefiles( "./airport_dbg", "accum" ); } @@ -151,7 +153,7 @@ int Taxiway::BuildBtg( tgpolygon_list& rwy_polys, tglightcontour_list& rwy_light rwy_polys.push_back( split ); if( shapefile_name.size() ) { - tgPolygon::ToShapefile( split, "./airport_dbg", std::string("postclip"), shapefile_name ); + tgShapefile::FromPolygon( split, "./airport_dbg", std::string("postclip"), shapefile_name ); } accum.Add( taxi_contour ); diff --git a/src/Airports/GenAirports850/taxiway.hxx b/src/Airports/GenAirports850/taxiway.hxx index 50da361d..7da527a6 100644 --- a/src/Airports/GenAirports850/taxiway.hxx +++ b/src/Airports/GenAirports850/taxiway.hxx @@ -1,7 +1,9 @@ #ifndef _TAXIWAY_H_ #define _TAXIWAY_H_ -#include +#include +#include +#include #include "apt_math.hxx" diff --git a/src/BuildTiles/Main/CMakeLists.txt b/src/BuildTiles/Main/CMakeLists.txt index 54318b3d..7afacdd3 100644 --- a/src/BuildTiles/Main/CMakeLists.txt +++ b/src/BuildTiles/Main/CMakeLists.txt @@ -27,7 +27,7 @@ set_target_properties(tg-construct PROPERTIES "DEFAULT_USGS_MAPFILE=\"${PKGDATADIR}/usgsmap.txt\";DEFAULT_PRIORITIES_FILE=\"${PKGDATADIR}/default_priorities.txt\"" ) target_link_libraries(tg-construct - Polygon + terragear Array landcover ${Boost_LIBRARIES} ${GDAL_LIBRARY} @@ -43,7 +43,7 @@ add_executable(cliptst cliptst.cxx) target_link_libraries(cliptst - Polygon + terragear ${GDAL_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} diff --git a/src/BuildTiles/Main/cliptst.cxx b/src/BuildTiles/Main/cliptst.cxx index 3374a58f..1e31f0c1 100644 --- a/src/BuildTiles/Main/cliptst.cxx +++ b/src/BuildTiles/Main/cliptst.cxx @@ -12,8 +12,9 @@ #include #include -#include -#include +#include +#include +#include //--------------------------------------------------------------------------- @@ -259,9 +260,8 @@ int main(int argc, char* argv[]) } if (show_svg) { - tgShapefile::Init(); - clipper_to_shapefile( subject, "./clptst_subject" ); - clipper_to_shapefile( clip, "./clptst_clip" ); + tgShapefile::FromClipper( subject, "./clptst_subject", "subject", "subject" ); + tgShapefile::FromClipper( clip, "./clptst_clip", "clip", "clip" ); } ClipType clipType; @@ -315,7 +315,7 @@ int main(int argc, char* argv[]) if( show_svg ) { cout << "Generating shapefile\n"; - clipper_to_shapefile( solution, "./clptst_solution" ); + tgShapefile::FromClipper( solution, "./clptst_solution", "solution", "solution" ); } } else cout << "failed.\n\n"; diff --git a/src/BuildTiles/Main/main.cxx b/src/BuildTiles/Main/main.cxx index c6362fcf..db26b8e1 100644 --- a/src/BuildTiles/Main/main.cxx +++ b/src/BuildTiles/Main/main.cxx @@ -71,9 +71,6 @@ int main(int argc, char **argv) { sglog().setLogLevels( SG_ALL, SG_INFO ); - // Initialize shapefile support (for debugging) - tgShapefile::Init(); - // // Parse the command-line arguments. // diff --git a/src/BuildTiles/Main/tgconstruct.hxx b/src/BuildTiles/Main/tgconstruct.hxx index 3ad9bcb4..8f4ad27b 100644 --- a/src/BuildTiles/Main/tgconstruct.hxx +++ b/src/BuildTiles/Main/tgconstruct.hxx @@ -34,7 +34,7 @@ #define TG_MAX_AREA_TYPES 128 #include -#include +#include #include #include "tglandclass.hxx" diff --git a/src/BuildTiles/Main/tgconstruct_cleanup.cxx b/src/BuildTiles/Main/tgconstruct_cleanup.cxx index 351afcb2..78f763ce 100644 --- a/src/BuildTiles/Main/tgconstruct_cleanup.cxx +++ b/src/BuildTiles/Main/tgconstruct_cleanup.cxx @@ -28,6 +28,8 @@ #include #include +#include + #include "tgconstruct.hxx" void TGConstruct::FixTJunctions( void ) { @@ -86,21 +88,21 @@ void TGConstruct::CleanClippedPolys() { poly = tgPolygon::Snap(poly, gSnap); if ( IsDebugShape( poly.GetId() ) ) { - tgPolygon::ToShapefile( poly, ds_name, "snapped", "" ); + tgShapefile::FromPolygon( poly, ds_name, "snapped", "" ); } // step 2 : remove_dups poly = tgPolygon::RemoveDups( poly ); if ( IsDebugShape( poly.GetId() ) ) { - tgPolygon::ToShapefile( poly, ds_name, "rem_dups", "" ); + tgShapefile::FromPolygon( poly, ds_name, "rem_dups", "" ); } // step 3 : remove_bad_contours poly = tgPolygon::RemoveBadContours( poly ); if ( IsDebugShape( poly.GetId() ) ) { - tgPolygon::ToShapefile( poly, ds_name, "rem_bcs", "" ); + tgShapefile::FromPolygon( poly, ds_name, "rem_bcs", "" ); } polys_clipped.set_poly(area, p, poly); diff --git a/src/BuildTiles/Main/tgconstruct_clip.cxx b/src/BuildTiles/Main/tgconstruct_clip.cxx index 2b2d8b4d..578e48de 100644 --- a/src/BuildTiles/Main/tgconstruct_clip.cxx +++ b/src/BuildTiles/Main/tgconstruct_clip.cxx @@ -26,6 +26,9 @@ #include +#include +#include + #include "tgconstruct.hxx" using std::string; @@ -97,9 +100,9 @@ bool TGConstruct::ClipLandclassPolys( void ) { // Dump the masks if ( debug_all || debug_shapes.size() || debug_areas.size() ) { - tgPolygon::ToShapefile( land_mask, ds_name, "land_mask", "" ); - tgPolygon::ToShapefile( water_mask, ds_name, "water_mask", "" ); - tgPolygon::ToShapefile( island_mask, ds_name, "island_mask", "" ); + tgShapefile::FromPolygon( land_mask, ds_name, "land_mask", "" ); + tgShapefile::FromPolygon( water_mask, ds_name, "water_mask", "" ); + tgShapefile::FromPolygon( island_mask, ds_name, "island_mask", "" ); } // process polygons in priority order @@ -130,7 +133,7 @@ bool TGConstruct::ClipLandclassPolys( void ) { sprintf(layer, "pre_clip_%d", polys_in.get_poly( i, j ).GetId() ); sprintf(name, "shape %d,%d", i,j); - tgPolygon::ToShapefile( tmp, ds_name, layer, name ); + tgShapefile::FromPolygon( tmp, ds_name, layer, name ); sprintf(layer, "pre_clip_accum_%d_%d", accum_idx, polys_in.get_poly( i, j ).GetId() ); accum.ToShapefiles( ds_name, layer ); @@ -145,7 +148,7 @@ bool TGConstruct::ClipLandclassPolys( void ) { sprintf(layer, "post_clip_%d", polys_in.get_poly( i, j ).GetId() ); sprintf(name, "shape %d,%d", i,j); - tgPolygon::ToShapefile( clipped, ds_name, layer, name ); + tgShapefile::FromPolygon( clipped, ds_name, layer, name ); } // only add to output list if the clip left us with a polygon @@ -190,7 +193,7 @@ bool TGConstruct::ClipLandclassPolys( void ) { char name[32]; for ( unsigned int i=0; i #include -#include -#include +#include +#include #include "tgconstruct.hxx" diff --git a/src/BuildTiles/Main/tgconstruct_tesselate.cxx b/src/BuildTiles/Main/tgconstruct_tesselate.cxx index 78a0a1fb..d809559e 100644 --- a/src/BuildTiles/Main/tgconstruct_tesselate.cxx +++ b/src/BuildTiles/Main/tgconstruct_tesselate.cxx @@ -27,6 +27,8 @@ #include +#include + #include "tgconstruct.hxx" void TGConstruct::TesselatePolys( void ) @@ -40,7 +42,7 @@ void TGConstruct::TesselatePolys( void ) tgPolygon poly = polys_clipped.get_poly(area, p ); if ( IsDebugShape( poly.GetId() ) ) { - tgPolygon::ToShapefile( poly, ds_name, "preteselate", "" ); + tgShapefile::FromPolygon( poly, ds_name, "preteselate", "" ); } tgRectangle rect = poly.GetBoundingBox(); diff --git a/src/BuildTiles/Main/tglandclass.hxx b/src/BuildTiles/Main/tglandclass.hxx index 39f9913d..ec895f7c 100644 --- a/src/BuildTiles/Main/tglandclass.hxx +++ b/src/BuildTiles/Main/tglandclass.hxx @@ -36,7 +36,7 @@ #include #include -#include +#include #define TG_MAX_AREA_TYPES 128 diff --git a/src/Lib/CMakeLists.txt b/src/Lib/CMakeLists.txt index fd831558..0543f73f 100644 --- a/src/Lib/CMakeLists.txt +++ b/src/Lib/CMakeLists.txt @@ -3,5 +3,5 @@ include_directories(${PROJECT_SOURCE_DIR}/src/Lib) add_subdirectory(Array) add_subdirectory(DEM) add_subdirectory(HGT) -add_subdirectory(Polygon) add_subdirectory(landcover) +add_subdirectory(terragear) diff --git a/src/Lib/Polygon/CMakeLists.txt b/src/Lib/Polygon/CMakeLists.txt deleted file mode 100644 index 4d5ba863..00000000 --- a/src/Lib/Polygon/CMakeLists.txt +++ /dev/null @@ -1,21 +0,0 @@ -include_directories(${GDAL_INCLUDE_DIR}) - -include( ${CGAL_USE_FILE} ) - -add_library(Polygon STATIC - clipper.cpp - clipper.hpp - polygon.cxx - polygon.hxx - rectangle.cxx - rectangle.hxx - surface.cxx - surface.hxx - tg_nodes.cxx - tg_nodes.hxx - tg_unique_geod.hxx - tg_unique_tgnode.hxx - tg_unique_vec2f.hxx - tg_unique_vec3d.hxx - tg_unique_vec3f.hxx -) \ No newline at end of file diff --git a/src/Lib/Polygon/polygon.cxx b/src/Lib/Polygon/polygon.cxx deleted file mode 100644 index d11618c7..00000000 --- a/src/Lib/Polygon/polygon.cxx +++ /dev/null @@ -1,2966 +0,0 @@ -// polygon.cxx -- polygon (with holes) management class -// -// Written by Curtis Olson, started March 1999. -// -// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program; if not, write to the Free Software -// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. -// - -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "polygon.hxx" - -#ifdef _MSC_VER -# define LONG_LONG_MAX LLONG_MAX -# define LONG_LONG_MIN LLONG_MIN -#endif - -///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// -// NEW IMPLEMENTATIONS -///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - -#include -#include - -#include -#include -#include -#include -#include -#include - - - -// TODO : Where does this belong -// Calculate theta of angle (a, b, c) -double tgPolygonCalcAngle(SGVec2d a, SGVec2d b, SGVec2d c) { - SGVec2d u, v; - double udist, vdist, uv_dot; - - // u . v = ||u|| * ||v|| * cos(theta) - - u = b-a; - udist = dist(a,b); - - v = b-c; - vdist = dist(b,c); - - uv_dot = dot(u,v); - - return acos(uv_dot / (udist * vdist)); -} - - - - - -// EXACT ( polygon list from contour ) -typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; -typedef Kernel::Point_2 Point_2; -typedef CGAL::Segment_2 Segment_2; - - -// INEXACT ( triangulate ) - -/* determining if a face is within the reulting poly */ -struct FaceInfo2 -{ - FaceInfo2() {} - int nesting_level; - - bool in_domain(){ - return nesting_level%2 == 1; - } -}; - -typedef CGAL::Exact_predicates_inexact_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::Constrained_Delaunay_triangulation_2 CDT; -typedef CDT::Point Point; -typedef CGAL::Polygon_2 Polygon_2; -typedef CGAL::Triangle_2 Triangle_2; - - -// GEOD HELPERS - maybe some math??? - -const double isEqual2D_Epsilon = 0.000001; - -static SGGeod SGGeod_snap( const SGGeod& in, double grid ) -{ - return SGGeod::fromDegM( grid * SGMisc::round( in.getLongitudeDeg()/grid ), - grid * SGMisc::round( in.getLatitudeDeg() /grid ), - grid * SGMisc::round( in.getElevationM() /grid ) ); -} - -static bool SGGeod_isEqual2D( const SGGeod& g0, const SGGeod& g1 ) -{ - return ( (fabs( g0.getLongitudeDeg() - g1.getLongitudeDeg() ) < isEqual2D_Epsilon) && - (fabs( g0.getLatitudeDeg() - g1.getLatitudeDeg() ) < isEqual2D_Epsilon ) ); -} - -static SGVec2d SGGeod_ToSGVec2d( const SGGeod& p ) -{ - return SGVec2d( p.getLongitudeDeg(), p.getLatitudeDeg() ); -} - -double SGGeod_CalculateTheta( const SGGeod& p0, const SGGeod& p1, const SGGeod& p2 ) -{ - SGVec2d v0, v1, v2; - SGVec2d u, v; - double udist, vdist, uv_dot; - - v0 = SGGeod_ToSGVec2d( p0 ); - v1 = SGGeod_ToSGVec2d( p1 ); - v2 = SGGeod_ToSGVec2d( p2 ); - - u = v1 - v0; - udist = norm(u); - - v = v1 - v2; - vdist = norm(v); - - uv_dot = dot(u, v); - - return acos( uv_dot / (udist * vdist) ); -} - -bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, - const std::vector& nodes, SGGeod& 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]; - - 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 = current; - 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]; - - 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 = current; - x_err_min = x_err; - } - } - } - } - } - - return found_node; -} - -void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, std::vector& nodes, tgContour& result, double bbEpsilon, double errEpsilon ) -{ - SGGeod new_pt; - - SG_LOG(SG_GENERAL, SG_BULK, " " << p0 << " <==> " << p1 ); - - bool found_extra = FindIntermediateNode( p0, p1, nodes, new_pt, bbEpsilon, errEpsilon ); - - if ( found_extra ) { - AddIntermediateNodes( p0, new_pt, nodes, result, bbEpsilon, errEpsilon ); - - result.AddNode( new_pt ); - SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt); - - AddIntermediateNodes( new_pt, p1, nodes, result, bbEpsilon, errEpsilon ); - } -} - -#define CLIPPER_FIXEDPT (10000000000000000) -#define CLIPPER_FIXED1M ( 90090) - -static ClipperLib::IntPoint SGGeod_ToClipper( const SGGeod& p ) -{ - ClipperLib::long64 x, y; - - x = (ClipperLib::long64)( p.getLongitudeDeg() * CLIPPER_FIXEDPT ); - y = (ClipperLib::long64)( p.getLatitudeDeg() * CLIPPER_FIXEDPT ); - - return ClipperLib::IntPoint( x, y ); -} - -static SGGeod SGGeod_FromClipper( const ClipperLib::IntPoint& p ) -{ - double lon, lat; - - lon = (double)( ((double)p.X) / (double)CLIPPER_FIXEDPT ); - lat = (double)( ((double)p.Y) / (double)CLIPPER_FIXEDPT ); - - return SGGeod::fromDeg( lon, lat ); -} - -static double Dist_ToClipper( double dist ) -{ - return ( dist * ( CLIPPER_FIXEDPT / CLIPPER_FIXED1M ) ); -} - -static tgRectangle BoundingBox_FromClipper( const ClipperLib::Polygons& subject ) -{ - ClipperLib::IntPoint min_pt, max_pt; - SGGeod min, max; - - min_pt.X = min_pt.Y = LONG_LONG_MAX; - max_pt.X = max_pt.Y = LONG_LONG_MIN; - - // for each polygon, we need to check the orientation, to set the hole flag... - for (unsigned int i=0; i max_pt.X ) { - max_pt.X = subject[i][j].X; - } - if ( subject[i][j].Y > max_pt.Y ) { - max_pt.Y = subject[i][j].Y; - } - } - } - - min = SGGeod_FromClipper( min_pt ); - max = SGGeod_FromClipper( max_pt ); - - return tgRectangle( min, max ); -} - - -// tgContour static functions - -tgContour tgContour::Snap( const tgContour& subject, double snap ) -{ - tgContour result; - SGGeod pt; - - for (unsigned int i = 0; i < subject.GetSize(); i++) { - pt = SGGeod_snap( subject.GetNode(i), snap ); - result.AddNode(pt); - } - result.SetHole( subject.GetHole() ); - - return result; -} - -double tgContour::GetMinimumAngle( void ) const -{ - unsigned int p1_index, p2_index, p3_index; - SGVec2d p1, p2, p3; - double angle; - double min_angle = 2.0 * SGD_PI; - unsigned int size = node_list.size(); - - SG_LOG(SG_GENERAL, SG_DEBUG, " tgContour::GetMinimumAngle() : contour size is " << size ); - - for ( unsigned int i = 0; i < size; i++ ) { - if ( i == 0) { - p1_index = size -1; - } else { - p1_index = i - 1; - } - - p2_index = i; - - if ( i == size ) { - p3_index = 0; - } else { - p3_index = i + 1; - } - - SG_LOG(SG_GENERAL, SG_DEBUG, " tgContour::GetMinimumAngle() iteration " << i << " of " << size << " p1 " << p1_index << " p2 " << p2_index << " p3 " << p3_index ); - - p1 = SGGeod_ToSGVec2d( node_list[p1_index] ); - p2 = SGGeod_ToSGVec2d( node_list[p2_index] ); - p3 = SGGeod_ToSGVec2d( node_list[p3_index] ); - - SG_LOG(SG_GENERAL, SG_DEBUG, " tgContour::GetMinimumAngle() iteration " << i << " of " << size << " p1 " << p1 << " p2 " << p2 << " p3 " << p3 ); - - angle = tgPolygonCalcAngle( p1, p2, p3 ); - - SG_LOG(SG_GENERAL, SG_DEBUG, " tgContour::GetMinimumAngle() iteration " << i << " of " << size << "angle is " << angle ); - - if ( angle < min_angle ) { - min_angle = angle; - } - } - - return min_angle; -} - -double tgContour::GetArea( void ) const -{ - double area = 0.0; - SGVec2d a, b; - unsigned int i, j; - - if ( node_list.size() ) { - j = node_list.size() - 1; - for (i=0; i 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; -} - -tgContour tgContour::RemoveDups( const tgContour& subject ) -{ - tgContour result; - - int iters = 0; - bool found; - - for ( unsigned int i = 0; i < subject.GetSize(); i++ ) - { - result.AddNode( subject.GetNode( i ) ); - } - result.SetHole( subject.GetHole() ); - - SG_LOG( SG_GENERAL, SG_DEBUG, "remove contour dups : original contour has " << result.GetSize() << " points" ); - - do - { - SG_LOG( SG_GENERAL, SG_DEBUG, "remove_contour_dups: start new iteration" ); - found = false; - - // Step 1 - find a neighboring duplicate points - SGGeod cur, next; - unsigned int cur_loc, next_loc; - - for ( unsigned int i = 0; i < result.GetSize() && !found; i++ ) { - if (i == result.GetSize() - 1 ) { - cur_loc = i; - cur = result.GetNode(i); - - next_loc = 0; - next = result.GetNode(0); - - SG_LOG( SG_GENERAL, SG_DEBUG, " cur is last point: " << cur << " next is first point: " << next ); - } else { - cur_loc = i; - cur = result.GetNode(i); - - next_loc = i+1; - next = result.GetNode(i+1); - SG_LOG( SG_GENERAL, SG_DEBUG, " cur is: " << cur << " next is : " << next ); - } - - if ( SGGeod_isEqual2D( cur, next ) ) { - // keep the point with higher Z - if ( cur.getElevationM() < next.getElevationM() ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_dups: erasing " << cur ); - result.RemoveNodeAt( cur_loc ); - // result.RemoveNodeAt( result.begin()+cur ); - } else { - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_dups: erasing " << next ); - result.RemoveNodeAt( next_loc ); - } - found = true; - } - } - - iters++; - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_dups : after " << iters << " iterations, contour has " << result.GetSize() << " points" ); - } while( found ); - - return result; - -} - -tgContour tgContour::SplitLongEdges( const tgContour& subject, double max_len ) -{ - SGGeod p0, p1; - double dist; - tgContour result; - - for ( unsigned i = 0; i < subject.GetSize() - 1; i++ ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "point = " << i); - - p0 = subject.GetNode( i ); - p1 = subject.GetNode( i + 1 ); - - SG_LOG(SG_GENERAL, SG_DEBUG, " " << p0 << " - " << p1); - - if ( fabs(p0.getLatitudeDeg()) < (90.0 - SG_EPSILON) || - fabs(p1.getLatitudeDeg()) < (90.0 - SG_EPSILON) ) - { - dist = SGGeodesy::distanceM( p0, p1 ); - SG_LOG(SG_GENERAL, SG_DEBUG, "distance = " << dist); - - if ( dist > max_len ) { - unsigned int segments = (int)(dist / max_len) + 1; - SG_LOG(SG_GENERAL, SG_DEBUG, "segments = " << segments); - - double dx = (p1.getLongitudeDeg() - p0.getLongitudeDeg()) / segments; - double dy = (p1.getLatitudeDeg() - p0.getLatitudeDeg()) / segments; - - for ( unsigned int j = 0; j < segments; j++ ) { - SGGeod tmp = SGGeod::fromDeg( p0.getLongitudeDeg() + dx * j, p0.getLatitudeDeg() + dy * j ); - SG_LOG(SG_GENERAL, SG_DEBUG, tmp); - result.AddNode( tmp ); - } - } else { - SG_LOG(SG_GENERAL, SG_DEBUG, p0); - result.AddNode( p0 ); - } - } else { - SG_LOG(SG_GENERAL, SG_DEBUG, p0); - result.AddNode( p0 ); - } - - // end of segment is beginning of next segment - } - - p0 = subject.GetNode( subject.GetSize() - 1 ); - p1 = subject.GetNode( 0 ); - - dist = SGGeodesy::distanceM( p0, p1 ); - SG_LOG(SG_GENERAL, SG_DEBUG, "distance = " << dist); - - if ( dist > max_len ) { - unsigned int segments = (int)(dist / max_len) + 1; - SG_LOG(SG_GENERAL, SG_DEBUG, "segments = " << segments); - - double dx = (p1.getLongitudeDeg() - p0.getLongitudeDeg()) / segments; - double dy = (p1.getLatitudeDeg() - p0.getLatitudeDeg()) / segments; - - for ( unsigned int i = 0; i < segments; i++ ) { - SGGeod tmp = SGGeod::fromDeg( p0.getLongitudeDeg() + dx * i, p0.getLatitudeDeg() + dy * i ); - SG_LOG(SG_GENERAL, SG_DEBUG, tmp); - result.AddNode( tmp ); - } - } else { - SG_LOG(SG_GENERAL, SG_DEBUG, p0); - result.AddNode( p0 ); - } - - // maintain original hole flag setting - result.SetHole( subject.GetHole() ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "split_long_edges() complete"); - - return result; -} - -tgContour tgContour::RemoveSpikes( const tgContour& subject ) -{ - tgContour result; - int iters = 0; - double theta; - bool found; - SGGeod cur, prev, next; - - for ( unsigned int i = 0; i < subject.GetSize(); i++ ) - { - result.AddNode( subject.GetNode( i ) ); - } - - SG_LOG(SG_GENERAL, SG_DEBUG, "remove contour spikes : original contour has " << result.GetSize() << " points" ); - - do - { - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_spikes: start new iteration"); - found = false; - - // Step 1 - find a duplicate point - for ( unsigned int i = 0; i < result.GetSize() && !found; i++ ) { - if (i == 0) { - SG_LOG(SG_GENERAL, SG_DEBUG, " cur is first point: " << i << ": " << result.GetNode(i) ); - cur = result.GetNode( 0 ); - prev = result.GetNode( result.GetSize()-1 ); - next = result.GetNode( 1 ); - } else if ( i == result.GetSize()-1 ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " cur is last point: " << i << ": " << result.GetNode(i) ); - cur = result.GetNode( i ); - prev = result.GetNode( i-1 ); - next = result.GetNode( 0 ); - } else { - SG_LOG(SG_GENERAL, SG_DEBUG, " cur is: " << i << ": " << result.GetNode(i) ); - cur = result.GetNode( i ); - prev = result.GetNode( i-1 ); - next = result.GetNode( i+1 ); - } - - theta = SGMiscd::rad2deg( SGGeod_CalculateTheta(prev, cur, next) ); - - if ( abs(theta) < 0.1 ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_spikes: (theta is " << theta << ") erasing " << i << " prev is " << prev << " cur is " << cur << " next is " << next ); - result.RemoveNodeAt( i ); - found = true; - } - } - - iters++; - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_spikes : after " << iters << " iterations, contour has " << result.GetSize() << " points" ); - } while( found ); - - return result; -} - -ClipperLib::Polygon tgContour::ToClipper( const tgContour& subject ) -{ - ClipperLib::Polygon contour; - - for ( unsigned int i=0; i::infinity(); - double miny = std::numeric_limits::infinity(); - double maxx = -std::numeric_limits::infinity(); - double maxy = -std::numeric_limits::infinity(); - - for (unsigned int i = 0; i < node_list.size(); i++) { - SGGeod pt = GetNode(i); - if ( pt.getLongitudeDeg() < minx ) { minx = pt.getLongitudeDeg(); } - if ( pt.getLongitudeDeg() > maxx ) { maxx = pt.getLongitudeDeg(); } - if ( pt.getLatitudeDeg() < miny ) { miny = pt.getLatitudeDeg(); } - if ( pt.getLatitudeDeg() > maxy ) { maxy = pt.getLatitudeDeg(); } - } - - min = SGGeod::fromDeg( minx, miny ); - max = SGGeod::fromDeg( maxx, maxy ); - - return tgRectangle( min, max ); -} - -tgPolygon tgContour::Union( const tgContour& subject, tgPolygon& clip ) -{ - tgPolygon result; - UniqueSGGeodSet all_nodes; - - /* before diff - gather all nodes */ - for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { - all_nodes.add( subject.GetNode(i) ); - } - - ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); - ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); - ClipperLib::Polygons clipper_result; - - ClipperLib::Clipper c; - c.Clear(); - c.AddPolygon(clipper_subject, ClipperLib::ptSubject); - c.AddPolygons(clipper_clip, ClipperLib::ptClip); - c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); - - result = tgPolygon::FromClipper( clipper_result ); - result = tgPolygon::AddColinearNodes( result, all_nodes ); - - return result; -} - -tgPolygon tgContour::Diff( const tgContour& subject, tgPolygon& clip ) -{ - tgPolygon result; - UniqueSGGeodSet all_nodes; - - /* before diff - gather all nodes */ - for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { - all_nodes.add( subject.GetNode(i) ); - } - - ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); - ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); - ClipperLib::Polygons clipper_result; - - ClipperLib::Clipper c; - c.Clear(); - c.AddPolygon(clipper_subject, ClipperLib::ptSubject); - c.AddPolygons(clipper_clip, ClipperLib::ptClip); - c.Execute(ClipperLib::ctDifference, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); - - result = tgPolygon::FromClipper( clipper_result ); - result = tgPolygon::AddColinearNodes( result, all_nodes ); - - return result; -} - -tgPolygon tgPolygon::Union( const tgContour& subject, tgPolygon& clip ) -{ - tgPolygon result; - UniqueSGGeodSet all_nodes; - - /* before diff - gather all nodes */ - for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { - all_nodes.add( subject.GetNode(i) ); - } - - ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); - ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); - ClipperLib::Polygons clipper_result; - - ClipperLib::Clipper c; - c.Clear(); - c.AddPolygon(clipper_subject, ClipperLib::ptSubject); - c.AddPolygons(clipper_clip, ClipperLib::ptClip); - c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); - - result = tgPolygon::FromClipper( clipper_result ); - result = tgPolygon::AddColinearNodes( result, all_nodes ); - - return result; -} - -tgContour tgContour::AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ) -{ - SGGeod p0, p1; - tgContour result; - std::vector& tmp_nodes = nodes.get_list(); - - 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, tmp_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, tmp_nodes, result, SG_EPSILON*10, SG_EPSILON*4 ); - - // maintain original hole flag setting - result.SetHole( subject.GetHole() ); - - return result; -} - -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, 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, 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 ) -{ - SGGeod p0, p1; - SGGeod new_pt; - std::vector tmp_nodes; - - tmp_nodes.push_back( node ); - for ( unsigned int n = 0; n < subject.GetSize()-1; n++ ) { - p0 = subject.GetNode( n ); - p1 = subject.GetNode( n+1 ); - - // add intermediate points - bool found_extra = FindIntermediateNode( p0, p1, tmp_nodes, new_pt, SG_EPSILON*10, SG_EPSILON*4 ); - if ( found_extra ) { - start = p0; - end = p1; - return true; - } - } - - // check last segment - p0 = subject.GetNode( subject.GetSize() - 1 ); - p1 = subject.GetNode( 0 ); - - // add intermediate points - bool found_extra = FindIntermediateNode( p0, p1, tmp_nodes, new_pt, SG_EPSILON*10, SG_EPSILON*4 ); - if ( found_extra ) { - start = p0; - end = p1; - return true; - } - - return false; -} - -void tgContour::SaveToGzFile( gzFile& fp ) const -{ - // Save the nodelist - sgWriteUInt( fp, node_list.size() ); - for (unsigned int i = 0; i < node_list.size(); i++) { - sgWriteGeod( fp, node_list[i] ); - } - - // and the hole flag - sgWriteInt( fp, (int)hole ); -} - -void tgContour::LoadFromGzFile( gzFile& fp ) -{ - unsigned int count; - SGGeod node; - - // Start Clean - Erase(); - - // Load the nodelist - sgReadUInt( fp, &count ); - for (unsigned int i = 0; i < count; i++) { - sgReadGeod( fp, node ); - node_list.push_back( node ); - } - - sgReadInt( fp, (int *)&hole ); -} - -std::ostream& operator<< ( std::ostream& output, const tgContour& subject ) -{ - // Save the data - output << "NumNodes: " << subject.node_list.size() << "\n"; - - for( unsigned int n=0; n= 3 ) { - /* keeping the contour */ - result.AddContour( contour ); - } - } - - return result; -} - -tgPolygon tgPolygon::RemoveCycles( const tgPolygon& subject ) -{ - tgPolygon result; - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - for ( unsigned int c = 0; c < subject.Contours(); c++ ) { - result.AddContour( tgContour::RemoveCycles( subject.GetContour( c ) ) ); - } - - return result; -} - -tgPolygon tgPolygon::SplitLongEdges( const tgPolygon& subject, double dist ) -{ - tgPolygon result; - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - for ( unsigned c = 0; c < subject.Contours(); c++ ) - { - result.AddContour( tgContour::SplitLongEdges( subject.GetContour(c), dist ) ); - } - - return result; -} - -tgPolygon tgPolygon::StripHoles( const tgPolygon& subject ) -{ - tgPolygon result; - UniqueSGGeodSet all_nodes; - - /* before diff - gather all nodes */ - for ( unsigned int i = 0; i < subject.Contours(); ++i ) { - for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { - all_nodes.add( subject.GetNode(i, j) ); - } - } - - ClipperLib::Polygons clipper_result; - ClipperLib::Clipper c; - c.Clear(); - - for ( unsigned int i = 0; i < subject.Contours(); i++ ) { - tgContour contour = subject.GetContour( i ); - if ( !contour.GetHole() ) { - c.AddPolygon( tgContour::ToClipper( contour ), ClipperLib::ptClip ); - } - } - c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); - - result = tgPolygon::FromClipper( clipper_result ); - result = tgPolygon::AddColinearNodes( result, all_nodes ); - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - return result; -} - -tgPolygon tgPolygon::Simplify( const tgPolygon& subject ) -{ - tgPolygon result; - UniqueSGGeodSet all_nodes; - - /* before diff - gather all nodes */ - for ( unsigned int i = 0; i < subject.Contours(); ++i ) { - for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { - all_nodes.add( subject.GetNode(i, j) ); - } - } - - ClipperLib::Polygons clipper_poly = tgPolygon::ToClipper( subject ); - SimplifyPolygons( clipper_poly ); - - result = tgPolygon::FromClipper( clipper_poly ); - result = tgPolygon::AddColinearNodes( result, all_nodes ); - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - return result; -} - -tgPolygon tgPolygon::RemoveTinyContours( const tgPolygon& subject ) -{ - double min_area = SG_EPSILON*SG_EPSILON; - tgPolygon result; - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - for ( unsigned int c = 0; c < subject.Contours(); c++ ) { - tgContour contour = subject.GetContour( c ); - double area = contour.GetArea(); - - if ( area >= min_area) { - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_tiny_contours NO - " << c << " area is " << area << " requirement is " << min_area); - result.AddContour( contour ); - } else { - SG_LOG(SG_GENERAL, SG_DEBUG, "remove_tiny_contours " << c << " area is " << area << ": removing"); - } - } - - return result; -} - -tgPolygon tgPolygon::RemoveSpikes( const tgPolygon& subject ) -{ - tgPolygon result; - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - for ( unsigned int c = 0; c < subject.Contours(); c++ ) { - result.AddContour( tgContour::RemoveSpikes( subject.GetContour(c) ) ); - } - - return result; -} - -ClipperLib::Polygons tgPolygon::ToClipper( const tgPolygon& subject ) -{ - ClipperLib::Polygons result; - - for ( unsigned int i=0; i *ipoint = CGAL::object_cast >(&result)) { - // handle the point intersection case with *ipoint. - return true; - } else { - if (const CGAL::Segment_2 *iseg = CGAL::object_cast >(&result)) { - // handle the segment intersection case with *iseg. - return false; - } else { - // handle the no intersection case. - return false; - } - } -} - -tgPolygon tgPolygon::Expand( const SGGeod& subject, double offset ) -{ - tgPolygon result; - tgContour contour; - SGGeod pt; - - pt = SGGeodesy::direct( subject, 90, offset/2.0 ); - double dlon = pt.getLongitudeDeg() - subject.getLongitudeDeg(); - - pt = SGGeodesy::direct( subject, 0, offset/2.0 ); - double dlat = pt.getLatitudeDeg() - subject.getLatitudeDeg(); - - contour.AddNode( SGGeod::fromDeg( subject.getLongitudeDeg() - dlon, subject.getLatitudeDeg() - dlat ) ); - contour.AddNode( SGGeod::fromDeg( subject.getLongitudeDeg() + dlon, subject.getLatitudeDeg() - dlat ) ); - 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; -} - -tgpolygon_list tgContour::ExpandToPolygons( const tgContour& subject, double width ) -{ - int turn_dir; - - SGGeod cur_inner; - SGGeod cur_outer; - SGGeod prev_inner; - SGGeod prev_outer; - SGGeod calc_inner; - SGGeod calc_outer; - - double last_end_v = 0.0f; - - tgContour expanded; - tgPolygon segment; - tgAccumulator accum; - tgpolygon_list result; - - // generate poly and texparam lists for each line segment - for (unsigned int i = 0; i < subject.GetSize(); i++) - { - last_end_v = 0.0f; - turn_dir = 0; - - sglog().setLogLevels( SG_ALL, SG_INFO ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "makePolygonsTP: calculating offsets for segment " << i); - - // for each point on the PointsList, generate a quad from - // start to next, offset by 1/2 width from the edge - if (i == 0) - { - // first point on the list - offset heading is 90deg - cur_outer = OffsetPointFirst( subject.GetNode(i), subject.GetNode(i+1), -width/2.0f ); - cur_inner = OffsetPointFirst( subject.GetNode(i), subject.GetNode(i+1), width/2.0f ); - } - else if (i == subject.GetSize()-1) - { - // last point on the list - offset heading is 90deg - cur_outer = OffsetPointLast( subject.GetNode(i-1), subject.GetNode(i), -width/2.0f ); - cur_inner = OffsetPointLast( subject.GetNode(i-1), subject.GetNode(i), width/2.0f ); - } - else - { - // middle section - cur_outer = OffsetPointMiddle( subject.GetNode(i-1), subject.GetNode(i), subject.GetNode(i+1), -width/2.0f, turn_dir ); - cur_inner = OffsetPointMiddle( subject.GetNode(i-1), subject.GetNode(i), subject.GetNode(i+1), width/2.0f, turn_dir ); - } - - if ( i > 0 ) - { - SGGeod prev_mp = midpoint( prev_outer, prev_inner ); - SGGeod cur_mp = midpoint( cur_outer, cur_inner ); - SGGeod intersect; - double heading; - double dist; - double az2; - - SGGeodesy::inverse( prev_mp, cur_mp, heading, az2, dist ); - - expanded.Erase(); - segment.Erase(); - - expanded.AddNode( prev_inner ); - expanded.AddNode( prev_outer ); - - // we need to extend one of the points so we're sure we don't create adjacent edges - if (turn_dir == 0) - { - // turned right - offset outer - - if ( getIntersection_cgal( prev_inner, prev_outer, cur_inner, cur_outer, intersect ) ) - { - // yes - make a triangle with inner edge = 0 - expanded.AddNode( cur_outer ); - cur_inner = prev_inner; - } - else - { - expanded.AddNode( cur_outer ); - expanded.AddNode( cur_inner ); - } - } - else - { - // turned left - offset inner - - if ( getIntersection_cgal( prev_inner, prev_outer, cur_inner, cur_outer, intersect ) ) - { - // yes - make a triangle with outer edge = 0 - expanded.AddNode( cur_inner ); - cur_outer = prev_outer; - } - else - { - expanded.AddNode( cur_outer ); - expanded.AddNode( cur_inner ); - } - } - - expanded.SetHole(false); - segment.AddContour(expanded); - segment.SetTexParams( prev_inner, width, 20.0f, heading ); - segment.SetTexLimits( 0, last_end_v, 1, 1 ); - segment.SetTexMethod( TG_TEX_BY_TPS_CLIPU, -1.0, 0.0, 1.0, 0.0 ); - result.push_back( segment ); - - last_end_v = 1.0f - (fmod( (double)(dist - last_end_v), (double)1.0f )); - } - - prev_outer = cur_outer; - prev_inner = cur_inner; - } - - sglog().setLogLevels( SG_ALL, SG_INFO ); - - return result; -} - -tgPolygon tgPolygon::Union( const tgPolygon& subject, tgPolygon& clip ) -{ - tgPolygon result; - UniqueSGGeodSet all_nodes; - std::ofstream dmpfile; - - /* before union - gather all nodes */ - for ( unsigned int i = 0; i < subject.Contours(); ++i ) { - for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { - all_nodes.add( subject.GetNode(i, j) ); - } - } - - ClipperLib::Polygons clipper_subject = tgPolygon::ToClipper( subject ); - ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); - ClipperLib::Polygons clipper_result; - - if ( clipper_dump ) { - dmpfile.open ("subject.txt"); - dmpfile << clipper_subject; - dmpfile.close(); - - dmpfile.open ("clip.txt"); - dmpfile << clipper_clip; - dmpfile.close(); - } - - ClipperLib::Clipper c; - c.Clear(); - c.AddPolygons(clipper_subject, ClipperLib::ptSubject); - c.AddPolygons(clipper_clip, ClipperLib::ptClip); - c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); - - if ( clipper_dump ) { - dmpfile.open ("result.txt"); - dmpfile << clipper_result; - dmpfile.close(); - } - - result = tgPolygon::FromClipper( clipper_result ); - result = tgPolygon::AddColinearNodes( result, all_nodes ); - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - return result; -} - -tgPolygon tgPolygon::Union( const tgpolygon_list& polys ) -{ - ClipperLib::Polygons clipper_result; - ClipperLib::Clipper c; - UniqueSGGeodSet all_nodes; - tgPolygon result; - - /* before union - gather all nodes */ - for ( unsigned int i=0; i::infinity(); - double miny = std::numeric_limits::infinity(); - double maxx = -std::numeric_limits::infinity(); - double maxy = -std::numeric_limits::infinity(); - - for ( unsigned int i = 0; i < Contours(); i++ ) { - for (unsigned int j = 0; j < ContourSize(i); j++) { - SGGeod pt = GetNode(i,j); - if ( pt.getLongitudeDeg() < minx ) { minx = pt.getLongitudeDeg(); } - if ( pt.getLongitudeDeg() > maxx ) { maxx = pt.getLongitudeDeg(); } - if ( pt.getLatitudeDeg() < miny ) { miny = pt.getLatitudeDeg(); } - if ( pt.getLatitudeDeg() > maxy ) { maxy = pt.getLatitudeDeg(); } - } - } - - min = SGGeod::fromDeg( minx, miny ); - max = SGGeod::fromDeg( maxx, maxy ); - - return tgRectangle( min, max ); -} - -void clipperToShapefile( ClipperLib::Polygons polys, const std::string& path, const std::string&layer, const std::string& name ) -{ - tgPolygon poly = tgPolygon::FromClipper(polys); - tgPolygon::ToShapefile( poly, path, layer, name); -} - -// Move slivers from in polygon to out polygon. -void tgPolygon::RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers ) -{ -#if 0 - // traverse each contour of the polygon and attempt to identify - // likely slivers - SG_LOG(SG_GENERAL, SG_DEBUG, "tgPolygon::RemoveSlivers()"); - - tgPolygon result; - tgContour contour; - int i; - - double angle_cutoff = 10.0 * SGD_DEGREES_TO_RADIANS; - double area_cutoff = 0.000000001; - double min_angle; - double area; - - // process contours in reverse order so deleting a contour doesn't - // foul up our sequence - for ( i = subject.Contours() - 1; i >= 0; --i ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "contour " << i ); - - contour = subject.GetContour(i); - - SG_LOG(SG_GENERAL, SG_DEBUG, " calc min angle for contour " << i); - min_angle = contour.GetMinimumAngle(); - SG_LOG(SG_GENERAL, SG_DEBUG, " min_angle (rad) = " << min_angle ); - - area = contour.GetArea(); - SG_LOG(SG_GENERAL, SG_DEBUG, " area = " << area ); - - if ( ((min_angle < angle_cutoff) && (area < area_cutoff)) || - ( area < area_cutoff / 10.0) ) - { - if ((min_angle < angle_cutoff) && (area < area_cutoff)) - { - SG_LOG(SG_GENERAL, SG_DEBUG, " WE THINK IT'S A SLIVER! - min angle < 10 deg, and area < 10 sq meters"); - } - else - { - SG_LOG(SG_GENERAL, SG_DEBUG, " WE THINK IT'S A SLIVER! - min angle > 10 deg, but area < 1 sq meters"); - } - - // Remove the sliver from source - subject.DeleteContourAt( i ); - - // And add it to the slive list if it isn't a hole - if ( !contour.GetHole() ) { - // move sliver contour to sliver list - SG_LOG(SG_GENERAL, SG_DEBUG, " Found SLIVER!"); - - slivers.push_back( contour ); - } - } - } -#endif -} - -tgcontour_list tgPolygon::MergeSlivers( tgpolygon_list& polys, tgcontour_list& sliver_list ) { - tgPolygon poly, result; - tgContour sliver; - tgContour contour; - tgcontour_list unmerged; - unsigned int original_contours, result_contours; - bool done; - - for ( unsigned int i = 0; i < sliver_list.size(); i++ ) { - sliver = sliver_list[i]; - SG_LOG(SG_GENERAL, SG_DEBUG, "Merging sliver = " << i ); - - sliver.SetHole( false ); - - done = false; - - // try to merge the slivers with the list of clipped polys - for ( unsigned int j = 0; j < polys.size() && !done; j++ ) { - poly = polys[j]; - original_contours = poly.Contours(); - result = tgPolygon::Union( sliver, poly ); - result_contours = result.Contours(); - - if ( original_contours == result_contours ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " FOUND a poly to merge the sliver with"); - result.SetMaterial( polys[j].GetMaterial() ); - result.SetTexParams( polys[j].GetTexParams() ); - polys[j] = result; - done = true; - } - } - - if ( !done ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "couldn't merge sliver " << i ); - unmerged.push_back( sliver ); - } - } - - return unmerged; -} - -tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, UniqueSGGeodSet& nodes ) -{ - tgPolygon result; - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - for ( unsigned int c = 0; c < subject.Contours(); c++ ) { - result.AddContour( tgContour::AddColinearNodes( subject.GetContour(c), nodes ) ); - } - - return result; -} - -tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, std::vector& nodes ) -{ - tgPolygon result; - - result.SetMaterial( subject.GetMaterial() ); - result.SetTexParams( subject.GetTexParams() ); - - for ( unsigned int c = 0; c < subject.Contours(); c++ ) { - result.AddContour( tgContour::AddColinearNodes( subject.GetContour(c), 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 ) -{ - bool found = false; - - for ( unsigned int c = 0; c < subject.Contours() && !found; c++ ) { - found = tgContour::FindColinearLine( subject.GetContour(c), node, start, end ); - } - - return found; -} - -SGGeod InterpolateElevation( const SGGeod& dst_node, const SGGeod& start, const SGGeod& end ) -{ - double total_dist = SGGeodesy::distanceM( start, end ); - double inter_dist = SGGeodesy::distanceM( start, dst_node ); - double delta = inter_dist/total_dist; - - double dest_elevation = start.getElevationM() + (delta * ( end.getElevationM() - start.getElevationM() )); - - return SGGeod::fromDegM( dst_node.getLongitudeDeg(), dst_node.getLatitudeDeg(), dest_elevation ); -} - - -void tgPolygon::InheritElevations( const tgPolygon& source ) -{ - UniqueSGGeodSet src_nodes; - - // build a list of points from the source polygon - for ( unsigned int i = 0; i < source.Contours(); ++i ) { - for ( unsigned int j = 0; j < source.ContourSize(i); ++j ) { - src_nodes.add( source.GetNode( i, j ) ); - } - } - - // traverse the dest polygon and build a mirror image but with - // elevations from the source polygon - for ( unsigned int i = 0; i < contours.size(); ++i ) { - for ( unsigned int j = 0; j < contours[i].GetSize(); ++j ) { - SGGeod dst_node = GetNode(i,j); - int index = src_nodes.find( dst_node ); - if ( index >= 0 ) { - SetNode( i, j, src_nodes.get_list()[index] ); - } else { - /* node not is source - we need to find the two points to interpolate from */ - SGGeod start, end, result; - if ( FindColinearLine( source, dst_node, start, end ) ) { - dst_node = InterpolateElevation( dst_node, start, end ); - SetNode( i, j, dst_node ); - } - } - } - } -} - -void tgPolygon::Texture( void ) -{ - SGGeod p; - SGVec2f t; - double x, y; - float tx, ty; - - SG_LOG(SG_GENERAL, SG_DEBUG, "Texture Poly with material " << material << " method " << tp.method << " tpref " << tp.ref << " heading " << tp.heading ); - - switch( tp.method ) { - case TG_TEX_BY_GEODE: - { - // The Simgear General texture coordinate routine takes a fan. - // Simgear could probably use a new function that just takes a Geod vector - // For now, just create an identity fan... - std::vector< int > node_idxs; - for (int i = 0; i < 3; i++) { - node_idxs.push_back(i); - } - - for ( unsigned int i = 0; i < triangles.size(); i++ ) { - std::vector< SGVec2f > tc_list; - std::vector< SGGeod > nodes; - - nodes = triangles[i].GetNodeList(); - tc_list = sgCalcTexCoords( tp.center_lat, nodes, node_idxs ); - triangles[i].SetTexCoordList( tc_list ); - } - } - break; - - case TG_TEX_BY_TPS_NOCLIP: - case TG_TEX_BY_TPS_CLIPU: - case TG_TEX_BY_TPS_CLIPV: - case TG_TEX_BY_TPS_CLIPUV: - { - for ( unsigned int i = 0; i < triangles.size(); i++ ) { - for ( unsigned int j = 0; j < 3; j++ ) { - p = triangles[i].GetNode( j ); - SG_LOG(SG_GENERAL, SG_DEBUG, "point = " << p); - - // - // 1. Calculate distance and bearing from the center of - // the poly - // - - // given alt, lat1, lon1, lat2, lon2, calculate starting - // and ending az1, az2 and distance (s). Lat, lon, and - // azimuth are in degrees. distance in meters - double az1, az2, dist; - SGGeodesy::inverse( tp.ref, p, az1, az2, dist ); - SG_LOG(SG_GENERAL, SG_DEBUG, "basic course = " << az2); - - // - // 2. Rotate this back into a coordinate system where Y - // runs the length of the poly and X runs crossways. - // - - double course = SGMiscd::normalizePeriodic(0, 360, az2 - tp.heading); - SG_LOG( SG_GENERAL, SG_DEBUG," course = " << course << " dist = " << dist ); - - // - // 3. Convert from polar to cartesian coordinates - // - - x = sin( course * SGD_DEGREES_TO_RADIANS ) * dist; - y = cos( course * SGD_DEGREES_TO_RADIANS ) * dist; - SG_LOG(SG_GENERAL, SG_DEBUG, " x = " << x << " y = " << y); - - // - // 4. Map x, y point into texture coordinates - // - float tmp; - - tmp = (float)x / (float)tp.width; - tx = tmp * (float)(tp.maxu - tp.minu) + (float)tp.minu; - SG_LOG(SG_GENERAL, SG_DEBUG, " (" << tx << ")"); - - // clip u? - if ( (tp.method == TG_TEX_BY_TPS_CLIPU) || (tp.method == TG_TEX_BY_TPS_CLIPUV) ) { - if ( tx < (float)tp.min_clipu ) { tx = (float)tp.min_clipu; } - if ( tx > (float)tp.max_clipu ) { tx = (float)tp.max_clipu; } - } - - tmp = (float)y / (float)tp.length; - ty = tmp * (float)(tp.maxv - tp.minv) + (float)tp.minv; - SG_LOG(SG_GENERAL, SG_DEBUG, " (" << ty << ")"); - - // clip v? - if ( (tp.method == TG_TEX_BY_TPS_CLIPV) || (tp.method == TG_TEX_BY_TPS_CLIPUV) ) { - if ( ty < (float)tp.min_clipv ) { ty = (float)tp.min_clipv; } - if ( ty > (float)tp.max_clipv ) { ty = (float)tp.max_clipv; } - } - - t = SGVec2f( tx, ty ); - SG_LOG(SG_GENERAL, SG_DEBUG, " (" << tx << ", " << ty << ")"); - - triangles[i].SetTexCoord( j, t ); - } - } - } - break; - } -} - -void tgPolygon::ToShapefile( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description ) -{ - void* ds_id = tgShapefile::OpenDatasource( datasource.c_str() ); - SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenDatasource returned " << (unsigned long)ds_id); - - OGRLayer* l_id = (OGRLayer *)tgShapefile::OpenLayer( ds_id, layer.c_str() ); - SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenLayer returned " << (unsigned long)l_id); - - OGRPolygon* polygon = new OGRPolygon(); - - SG_LOG(SG_GENERAL, SG_DEBUG, "subject has " << subject.Contours() << " contours "); - - for ( unsigned int i = 0; i < subject.Contours(); i++ ) { - bool skip_ring=false; - tgContour contour = subject.GetContour( i ); - - if (contour.GetSize() < 3) { - SG_LOG(SG_GENERAL, SG_DEBUG, "Polygon with less than 3 points"); - skip_ring=true; - } - - // FIXME: Current we ignore the hole-flag and instead assume - // that the first ring is not a hole and the rest - // are holes - OGRLinearRing *ring=new OGRLinearRing(); - for (unsigned int pt = 0; pt < contour.GetSize(); pt++) { - OGRPoint *point=new OGRPoint(); - - point->setX( contour.GetNode(pt).getLongitudeDeg() ); - point->setY( contour.GetNode(pt).getLatitudeDeg() ); - point->setZ( 0.0 ); - ring->addPoint(point); - } - ring->closeRings(); - - if (!skip_ring) { - polygon->addRingDirectly(ring); - } - - OGRFeature* feature = NULL; - feature = new OGRFeature( l_id->GetLayerDefn() ); - feature->SetField("ID", description.c_str()); - feature->SetGeometry(polygon); - if( l_id->CreateFeature( feature ) != OGRERR_NONE ) - { - SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile"); - } - OGRFeature::DestroyFeature(feature); - } - - // close after each write - ds_id = tgShapefile::CloseDatasource( ds_id ); -} - -tgPolygon tgPolygon::FromOGR( const OGRPolygon* subject ) -{ - OGRLinearRing const *ring = subject->getExteriorRing(); - tgContour contour; - tgPolygon result; - - for (int i = 0; i < ring->getNumPoints(); i++) { - contour.AddNode( SGGeod::fromDegM( ring->getX(i), ring->getY(i), ring->getZ(i)) ); - } - contour.SetHole( false ); - result.AddContour( contour ); - - // then add the inner rings - for ( int j = 0 ; j < subject->getNumInteriorRings(); j++ ) { - ring = subject->getInteriorRing( j ); - contour.Erase(); - - for (int i = 0; i < ring->getNumPoints(); i++) { - contour.AddNode( SGGeod::fromDegM( ring->getX(i), ring->getY(i), ring->getZ(i)) ); - } - contour.SetHole( true ); - result.AddContour( contour ); - } - result.SetTexMethod( TG_TEX_BY_GEODE ); - - return result; -} - -void tgPolygon::SaveToGzFile( gzFile& fp ) const -{ - // Save the contours - sgWriteUInt( fp, contours.size() ); - for (unsigned int i = 0; i < contours.size(); i++) { - contours[i].SaveToGzFile( fp ); - } - - // Save the triangles - sgWriteUInt( fp, triangles.size() ); - for (unsigned int i = 0; i < triangles.size(); i++) { - triangles[i].SaveToGzFile( fp ); - } - - // Save the tex params - tp.SaveToGzFile( fp ); - - // and the rest - sgWriteString( fp, material.c_str() ); - sgWriteString( fp, flag.c_str() ); - sgWriteInt( fp, (int)preserve3d ); -} - -void tgPolygon::LoadFromGzFile( gzFile& fp ) -{ - unsigned int count; - tgContour contour; - tgTriangle triangle; - char *strbuff; - - // Start clean - Erase(); - - // Load the contours - sgReadUInt( fp, &count ); - for (unsigned int i = 0; i < count; i++) { - contour.LoadFromGzFile( fp ); - AddContour(contour); - } - - // load the triangles - sgReadUInt( fp, &count ); - for (unsigned int i = 0; i < count; i++) { - triangle.LoadFromGzFile( fp ); - AddTriangle(triangle); - } - - // Load the tex params - tp.LoadFromGzFile( fp ); - - // and the rest - sgReadString( fp, &strbuff ); - if ( strbuff ) { - material = strbuff; - delete strbuff; - } - - sgReadString( fp, &strbuff ); - if ( strbuff ) { - flag = strbuff; - delete strbuff; - } - - sgReadInt( fp, (int *)&preserve3d ); -} - -#if 0 -////////////////////////////// CHOP //////////////////////////////// - -// initialize the unique polygon index counter stored in path -static long int poly_index = 0; -static std::string poly_path; - -bool tgPolygon::ChopIdxInit( const std::string& path ) -{ - poly_path = path; - FILE *fp = fopen( poly_path.c_str(), "r" ); - - if ( fp == NULL ) { - SG_LOG(SG_GENERAL, SG_DEBUG, "Warning: cannot open " << path); - poly_index = 0; - return false; - } - - fscanf( fp, "%ld", &poly_index ); - fclose( fp ); - - return true; -} -#endif - -/************************ TESSELATION ***********************************/ - -void tg_mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list& border ) -{ - if(start->info().nesting_level != -1) { - return; - } - - std::list queue; - queue.push_back(start); - - while( !queue.empty() ){ - CDT::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); - if(n->info().nesting_level == -1) { - if(ct.is_constrained(e)) border.push_back(e); - else queue.push_back(n); - } - } - } - } -} - -//explore set of facets connected with non constrained edges, -//and attribute to each such set a nesting level. -//We start from facets incident to the infinite vertex, with a nesting -//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. -void tg_mark_domains(CDT& cdt) -{ - for(CDT::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; - tg_mark_domains(cdt, cdt.infinite_face(), index++, border); - while(! border.empty()) { - CDT::Edge e = border.front(); - border.pop_front(); - CDT::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); - } - } -} - -void tg_insert_polygon(CDT& cdt,const Polygon_2& polygon) -{ - if ( polygon.is_empty() ) return; - - CDT::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); - cdt.insert_constraint(vh,v_prev); - v_prev=vh; - } -} - -void tgPolygon::Tesselate( const std::vector& extra ) -{ - CDT cdt; - - SG_LOG( SG_GENERAL, SG_DEBUG, "Tess with extra" ); - - // Bail right away if polygon is empty - if ( contours.size() != 0 ) { - // First, convert the extra points to cgal Points - std::vector points; - points.reserve(extra.size()); - for (unsigned int n = 0; n < extra.size(); n++) { - points.push_back( Point(extra[n].getLongitudeDeg(), extra[n].getLatitudeDeg() ) ); - } - - cdt.insert(points.begin(), points.end()); - - // then insert each polygon as a constraint into the triangulation - for ( unsigned int c = 0; c < contours.size(); c++ ) { - tgContour contour = contours[c]; - Polygon_2 poly; - - for (unsigned int n = 0; n < contour.GetSize(); n++ ) { - SGGeod node = contour.GetNode(n); - poly.push_back( Point( node.getLongitudeDeg(), node.getLatitudeDeg() ) ); - } - - 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) { - 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() ); - - AddTriangle( p0, p1, p2 ); - - ++count; - } - } - } -} - -void tgPolygon::Tesselate() -{ - CDT cdt; - - SG_LOG( SG_GENERAL, SG_DEBUG, "Tess" ); - - // Bail right away if polygon is empty - if ( contours.size() != 0 ) { - // insert each polygon as a constraint into the triangulation - for ( unsigned int c = 0; c < contours.size(); c++ ) { - tgContour contour = contours[c]; - Polygon_2 poly; - - for (unsigned int n = 0; n < contour.GetSize(); n++ ) { - SGGeod node = contour.GetNode(n); - SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : Adding GEOD " << node); - poly.push_back( Point( node.getLongitudeDeg(), node.getLatitudeDeg() ) ); - } - - 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) { - 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() ); - - AddTriangle( p0, p1, p2 ); - - ++count; - } else { - SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : face not in domain"); - } - } - } - else - { - SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : no contours" ); - } -} - -// Friends for serialization -std::ostream& operator<< ( std::ostream& output, const tgPolygon& subject ) -{ - // Save the data - output << "NumContours: " << subject.contours.size() << "\n"; - - for( unsigned int c=0; c= -89.0) && (c.getLatitudeDeg() < 89.0) ) { - min = SGGeod::fromDeg( c.getLongitudeDeg() - span/2.0, c.getLatitudeDeg() - SG_HALF_BUCKET_SPAN ); - max = SGGeod::fromDeg( c.getLongitudeDeg() + span/2.0, c.getLatitudeDeg() + SG_HALF_BUCKET_SPAN ); - } else if ( c.getLatitudeDeg() < -89.0) { - min = SGGeod::fromDeg( -90.0, -180.0 ); - max = SGGeod::fromDeg( -89.0, 180.0 ); - } else if ( c.getLatitudeDeg() >= 89.0) { - min = SGGeod::fromDeg( 89.0, -180.0 ); - max = SGGeod::fromDeg( 90.0, 180.0 ); - } else { - SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << c.getLatitudeDeg() ); - } - - SG_LOG( SG_GENERAL, SG_DEBUG, " (" << min << ") (" << max << ")" ); - - // set up clipping tile - base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), min.getLatitudeDeg()) ); - base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), min.getLatitudeDeg()) ); - base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), max.getLatitudeDeg()) ); - base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), max.getLatitudeDeg()) ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "shape contours = " << subject.Contours() ); - for ( unsigned int ii = 0; ii < subject.Contours(); ii++ ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << subject.GetContour(ii).GetHole() ); - } - - result = tgPolygon::Intersect( subject, base ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "result contours = " << result.Contours() ); - for ( unsigned int ii = 0; ii < result.Contours(); ii++ ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << result.GetContour(ii).GetHole() ); - } - - if ( subject.GetPreserve3D() ) { - result.InheritElevations( subject ); - } - - if ( result.Contours() > 0 ) { - result.SetPreserve3D( subject.GetPreserve3D() ); - result.SetTexParams( subject.GetTexParams() ); - if ( subject.GetTexMethod() == TG_TEX_BY_GEODE ) { - // need to set center latitude for geodetic texturing - result.SetTexMethod( TG_TEX_BY_GEODE, b.get_center_lat() ); - } - result.SetFlag(type); - - lock.lock(); - bp_map[b.gen_index()].push_back( result ); - lock.unlock(); - } -} - -void tgChopper::Add( const tgPolygon& subject, const std::string& type ) -{ - tgRectangle bb; - SGGeod p; - - // bail out immediately if polygon is empty - if ( subject.Contours() == 0 ) - return; - - bb = subject.GetBoundingBox(); - - SG_LOG( SG_GENERAL, SG_DEBUG, " min = " << bb.getMin() << " max = " << bb.getMax() ); - - // find buckets for min, and max points of convex hull. - // note to self: self, you should think about checking for - // polygons that span the date line - SGBucket b_min( bb.getMin() ); - SGBucket b_max( bb.getMax() ); - SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket min = " << b_min ); - SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket max = " << b_max ); - - if ( b_min == b_max ) { - // shape entirely contained in a single bucket, write and bail - Clip( subject, type, b_min ); - return; - } - - SGBucket b_cur; - int dx, dy; - - sgBucketDiff(b_min, b_max, &dx, &dy); - SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" ); - SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy ); - - if ( (dx > 2880) || (dy > 1440) ) - throw sg_exception("something is really wrong in split_polygon()!!!!"); - - if ( dy <= 1 ) { - // we are down to at most two rows, write each column and then bail - double min_center_lat = b_min.get_center_lat(); - double min_center_lon = b_min.get_center_lon(); - for ( int j = 0; j <= dy; ++j ) { - for ( int i = 0; i <= dx; ++i ) { - b_cur = sgBucketOffset(min_center_lon, min_center_lat, i, j); - Clip( subject, type, b_cur ); - } - } - return; - } - - // we have two or more rows left, split in half (along a - // horizontal dividing line) and recurse with each half - - // find mid point (integer math) - int mid = (dy + 1) / 2 - 1; - - // determine horizontal clip line - SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, mid); - double clip_line = b_clip.get_center_lat(); - if ( (clip_line >= -90.0 + SG_HALF_BUCKET_SPAN) - && (clip_line < 90.0 - SG_HALF_BUCKET_SPAN) ) - clip_line += SG_HALF_BUCKET_SPAN; - else if ( clip_line < -89.0 ) - clip_line = -89.0; - else if ( clip_line >= 89.0 ) - clip_line = 90.0; - else { - SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << clip_line ); - } - - { - // - // Crop bottom area (hopefully by putting this in it's own - // scope we can shorten the life of some really large data - // structures to reduce memory use) - // - - SG_LOG( SG_GENERAL, SG_DEBUG, "Generating bottom half (" << bb.getMin().getLatitudeDeg() << "-" << clip_line << ")" ); - - tgPolygon bottom, bottom_clip; - - bottom.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMin().getLatitudeDeg()) ); - bottom.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMin().getLatitudeDeg()) ); - bottom.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) ); - bottom.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) ); - - bottom_clip = tgPolygon::Intersect( subject, bottom ); - - // the texparam should be constant over each clipped poly. - // when they are reassembled, we want the texture map to - // be seamless - Add( bottom_clip, type ); - } - - { - // - // Crop top area (hopefully by putting this in it's own scope - // we can shorten the life of some really large data - // structures to reduce memory use) - // - - SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half (" << clip_line << "-" << bb.getMax().getLatitudeDeg() << ")" ); - - tgPolygon top, top_clip; - - top.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) ); - top.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) ); - top.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMax().getLatitudeDeg()) ); - top.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMax().getLatitudeDeg()) ); - - top_clip = tgPolygon::Intersect( subject, top ); - - if ( top_clip.TotalNodes() == subject.TotalNodes() ) { - SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half - total nodes is the same after clip" << subject.TotalNodes() ); - exit(0); - } - - Add( top_clip, type ); - } -} - -long int tgChopper::GenerateIndex( std::string path ) -{ - std::string index_file = path + "/chop.idx"; - long int index = 0; - - //Open or create the named mutex - boost::interprocess::named_mutex mutex(boost::interprocess::open_or_create, "tgChopper_index2"); - { -// SG_LOG(SG_GENERAL, SG_ALERT, "getting lock"); - boost::interprocess::scoped_lock lock(mutex); -// SG_LOG(SG_GENERAL, SG_ALERT, " - got it"); - - /* first try to read the file */ - FILE *fp = fopen( index_file.c_str(), "r+" ); - if ( fp == NULL ) { - /* doesn't exist - create it */ - fp = fopen( index_file.c_str(), "w" ); - if ( fp == NULL ) { - SG_LOG(SG_GENERAL, SG_ALERT, "Error cannot open Index file " << index_file << " for writing"); - boost::interprocess::named_mutex::remove("tgChopper_index2"); - exit( 0 ); - } - } else { - fread( (void*)&index, sizeof(long int), 1, fp ); -// SG_LOG(SG_GENERAL, SG_ALERT, " SUCCESS READING INDEX FILE - READ INDEX " << index ); - } - - index++; - - rewind( fp ); - fwrite( (void*)&index, sizeof(long int), 1, fp ); - fclose( fp ); - } - - boost::interprocess::named_mutex::remove("tgChopper_index2"); - - return index; -} - -void tgChopper::Save( void ) -{ - // traverse the bucket list - bucket_polys_map_interator it; - char tile_name[16]; - char poly_ext[16]; - - for (it=bp_map.begin(); it != bp_map.end(); it++) { - SGBucket b( (*it).first ); - tgpolygon_list const& polys = (*it).second; - - std::string path = root_path + "/" + b.gen_base_path(); - sprintf( tile_name, "%ld", b.gen_index() ); - - std::string polyfile = path + "/" + tile_name; - - SGPath sgp( polyfile ); - sgp.create_dir( 0755 ); - - long int poly_index = GenerateIndex( path ); - - sprintf( poly_ext, "%ld", poly_index ); - polyfile = polyfile + "." + poly_ext; - - gzFile fp; - if ( (fp = gzopen( polyfile.c_str(), "wb9" )) == NULL ) { - SG_LOG( SG_GENERAL, SG_INFO, "ERROR: opening " << polyfile.c_str() << " for writing!" ); - return; - } - - /* Write polys to the file */ - sgWriteUInt( fp, polys.size() ); - for ( unsigned int i=0; iGetDriverByName(format_name); - if (!ogrdriver) { - SG_LOG(SG_GENERAL, SG_ALERT, "Unknown datasource format driver: " << format_name); - exit(1); - } - - datasource = ogrdriver->Open(datasource_name, TRUE); - - if (!datasource) { - datasource = ogrdriver->CreateDataSource(datasource_name, NULL); - } - - if (!datasource) { - SG_LOG(SG_GENERAL, SG_ALERT, "Unable to open or create datasource: " << datasource_name); - exit(1); - } - - return (void*)datasource; -} - -void* tgShapefile::OpenLayer( void* ds_id, const char* layer_name ) { - OGRDataSource* datasource = (OGRDataSource *)ds_id; - OGRLayer* layer; - - OGRSpatialReference srs; - srs.SetWellKnownGeogCS("WGS84"); - - layer = datasource->GetLayerByName(layer_name); - - if (!layer) { - layer = datasource->CreateLayer( layer_name, &srs, wkbPolygon25D, NULL); - - OGRFieldDefn descriptionField("ID", OFTString); - descriptionField.SetWidth(128); - - if( layer->CreateField( &descriptionField ) != OGRERR_NONE ) { - SG_LOG(SG_GENERAL, SG_ALERT, "Creation of field 'Description' failed"); - } - } - - if (!layer) { - SG_LOG(SG_GENERAL, SG_ALERT, "Creation of layer '" << layer_name << "' failed"); - return NULL; - } - - return (void*)layer; -} - -void tgShapefile::CloseLayer( void* l_id ) -{ - //OGRLayer::DestroyLayer( layer ); -} - -void* tgShapefile::CloseDatasource( void* ds_id ) -{ - OGRDataSource* datasource = (OGRDataSource *)ds_id; - - OGRDataSource::DestroyDataSource( datasource ); - - return (void *)-1; -} - -void clipper_to_shapefile( ClipperLib::Polygons polys, char* ds ) -{ -#if 0 - ClipperLib::Polygons contour; - TGPolygon tgcontour; - char layer[32]; - - void* ds_id = tgShapefile::OpenDatasource( ds ); - - for (unsigned int i = 0; i < polys.size(); ++i) { - if ( Orientation( polys[i] ) ) { - sprintf( layer, "%04d_boundary", i ); - } else { - sprintf( layer, "%04d_hole", i ); - } - - void* l_id = tgShapefile::OpenLayer( ds_id, layer ); - contour.clear(); - contour.push_back( polys[i] ); - - tgcontour.erase(); - make_tg_poly_from_clipper( contour, &tgcontour ); - - tgShapefile::CreateFeature( ds_id, l_id, tgcontour, "contour" ); - } - - // close after each write - ds_id = tgShapefile::CloseDatasource( ds_id ); -#endif -} diff --git a/src/Lib/terragear/CMakeLists.txt b/src/Lib/terragear/CMakeLists.txt new file mode 100644 index 00000000..f9c9a257 --- /dev/null +++ b/src/Lib/terragear/CMakeLists.txt @@ -0,0 +1,35 @@ +include_directories(${GDAL_INCLUDE_DIR}) + +include( ${CGAL_USE_FILE} ) + +add_library(terragear STATIC + clipper.cpp + clipper.hpp + tg_accumulator.cxx + tg_accumulator.hxx + tg_chopper.cxx + tg_chopper.hxx + tg_contour.cxx + tg_contour.hxx + tg_light.hxx + tg_misc.cxx + tg_misc.hxx + tg_nodes.cxx + tg_nodes.hxx + tg_polygon.cxx + tg_polygon.hxx + tg_polygon_clean.cxx + tg_polygon_clip.cxx + tg_polygon_tesselate.cxx + tg_rectangle.cxx + tg_rectangle.hxx + tg_shapefile.cxx + tg_shapefile.hxx + tg_surface.cxx + tg_surface.hxx + tg_unique_geod.hxx + tg_unique_tgnode.hxx + tg_unique_vec2f.hxx + tg_unique_vec3d.hxx + tg_unique_vec3f.hxx +) \ No newline at end of file diff --git a/src/Lib/Polygon/TNT/jama_qr.h b/src/Lib/terragear/TNT/jama_qr.h similarity index 100% rename from src/Lib/Polygon/TNT/jama_qr.h rename to src/Lib/terragear/TNT/jama_qr.h diff --git a/src/Lib/Polygon/TNT/tnt_array1d.h b/src/Lib/terragear/TNT/tnt_array1d.h similarity index 100% rename from src/Lib/Polygon/TNT/tnt_array1d.h rename to src/Lib/terragear/TNT/tnt_array1d.h diff --git a/src/Lib/Polygon/TNT/tnt_array2d.h b/src/Lib/terragear/TNT/tnt_array2d.h similarity index 100% rename from src/Lib/Polygon/TNT/tnt_array2d.h rename to src/Lib/terragear/TNT/tnt_array2d.h diff --git a/src/Lib/Polygon/TNT/tnt_i_refvec.h b/src/Lib/terragear/TNT/tnt_i_refvec.h similarity index 100% rename from src/Lib/Polygon/TNT/tnt_i_refvec.h rename to src/Lib/terragear/TNT/tnt_i_refvec.h diff --git a/src/Lib/Polygon/TNT/tnt_math_utils.h b/src/Lib/terragear/TNT/tnt_math_utils.h similarity index 100% rename from src/Lib/Polygon/TNT/tnt_math_utils.h rename to src/Lib/terragear/TNT/tnt_math_utils.h diff --git a/src/Lib/Polygon/clipper.cpp b/src/Lib/terragear/clipper.cpp similarity index 100% rename from src/Lib/Polygon/clipper.cpp rename to src/Lib/terragear/clipper.cpp diff --git a/src/Lib/Polygon/clipper.hpp b/src/Lib/terragear/clipper.hpp similarity index 100% rename from src/Lib/Polygon/clipper.hpp rename to src/Lib/terragear/clipper.hpp diff --git a/src/Lib/terragear/tg_accumulator.cxx b/src/Lib/terragear/tg_accumulator.cxx new file mode 100644 index 00000000..f08933ef --- /dev/null +++ b/src/Lib/terragear/tg_accumulator.cxx @@ -0,0 +1,131 @@ +#include + +#include "tg_accumulator.hxx" +#include "tg_shapefile.hxx" +#include "tg_misc.hxx" + +tgPolygon tgAccumulator::Diff( const tgContour& subject ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { + all_nodes.add( subject.GetNode(i) ); + } + + unsigned int num_hits = 0; + tgRectangle box1 = subject.GetBoundingBox(); + + ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); + ClipperLib::Polygons clipper_result; + + ClipperLib::Clipper c; + c.Clear(); + + c.AddPolygon(clipper_subject, ClipperLib::ptSubject); + + // clip result against all polygons in the accum that intersect our bb + for (unsigned int i=0; i < accum.size(); i++) { + tgRectangle box2 = BoundingBox_FromClipper( accum[i] ); + + if ( box2.intersects(box1) ) + { + c.AddPolygons(accum[i], ClipperLib::ptClip); + num_hits++; + } + } + + if (num_hits) { + if ( !c.Execute(ClipperLib::ctDifference, clipper_result, ClipperLib::pftNonZero, ClipperLib::pftNonZero) ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Diff With Accumulator returned FALSE" ); + exit(-1); + } + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + } else { + result.AddContour( subject ); + } + + return result; +} + +void tgAccumulator::Add( const tgContour& subject ) +{ + tgPolygon poly; + poly.AddContour( subject ); + + ClipperLib::Polygons clipper_subject = tgPolygon::ToClipper( poly ); + accum.push_back( clipper_subject ); +} + +void tgAccumulator::ToShapefiles( const std::string& path, const std::string& layer_prefix ) +{ + char shapefile[16]; + char layer[16]; + + for (unsigned int i=0; i < accum.size(); i++) { + sprintf( layer, "%s_%d", layer_prefix.c_str(), i ); + sprintf( shapefile, "accum_%d", i ); + tgShapefile::FromClipper( accum[i], path, layer, std::string(shapefile) ); + } +} + +tgPolygon tgAccumulator::Diff( const tgPolygon& subject ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.Contours(); ++i ) { + for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { + all_nodes.add( subject.GetNode(i, j) ); + } + } + + unsigned int num_hits = 0; + tgRectangle box1 = subject.GetBoundingBox(); + + ClipperLib::Polygons clipper_subject = tgPolygon::ToClipper( subject ); + ClipperLib::Polygons clipper_result; + + ClipperLib::Clipper c; + c.Clear(); + + c.AddPolygons(clipper_subject, ClipperLib::ptSubject); + + // clip result against all polygons in the accum that intersect our bb + for (unsigned int i=0; i < accum.size(); i++) { + tgRectangle box2 = BoundingBox_FromClipper( accum[i] ); + + if ( box2.intersects(box1) ) + { + c.AddPolygons(accum[i], ClipperLib::ptClip); + num_hits++; + } + } + + if (num_hits) { + if ( !c.Execute(ClipperLib::ctDifference, clipper_result, ClipperLib::pftNonZero, ClipperLib::pftNonZero) ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Diff With Accumulator returned FALSE" ); + exit(-1); + } + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + // Make sure we keep texturing info + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + } else { + result = subject; + } + + return result; +} + +void tgAccumulator::Add( const tgPolygon& subject ) +{ + ClipperLib::Polygons clipper_subject = tgPolygon::ToClipper( subject ); + accum.push_back( clipper_subject ); +} diff --git a/src/Lib/terragear/tg_accumulator.hxx b/src/Lib/terragear/tg_accumulator.hxx new file mode 100644 index 00000000..d80192ed --- /dev/null +++ b/src/Lib/terragear/tg_accumulator.hxx @@ -0,0 +1,25 @@ +#ifndef _TGACCUMULATOR_HXX +#define _TGACCUMULATOR_HXX + +#include "tg_polygon.hxx" +#include "tg_contour.hxx" +#include "clipper.hpp" + +class tgAccumulator +{ +public: + tgPolygon Diff( const tgContour& subject ); + tgPolygon Diff( const tgPolygon& subject ); + + void Add( const tgContour& subject ); + void Add( const tgPolygon& subject ); + + void ToShapefiles( const std::string& path, const std::string& layer ); + +private: + typedef std::vector < ClipperLib::Polygons > clipper_polygons_list; + + clipper_polygons_list accum; +}; + +#endif // _TGACCUMULATOR_HXX \ No newline at end of file diff --git a/src/Lib/terragear/tg_chopper.cxx b/src/Lib/terragear/tg_chopper.cxx new file mode 100644 index 00000000..1a6e8167 --- /dev/null +++ b/src/Lib/terragear/tg_chopper.cxx @@ -0,0 +1,274 @@ +#include +#include + +#include +#include +#include +#include + + +#include "tg_chopper.hxx" + +void tgChopper::Clip( const tgPolygon& subject, + const std::string& type, + SGBucket& b ) +{ + // p; + + SGGeod min, max; + SGGeod c = b.get_center(); + double span = b.get_width(); + tgPolygon base, result; + + // calculate bucket dimensions + if ( (c.getLatitudeDeg() >= -89.0) && (c.getLatitudeDeg() < 89.0) ) { + min = SGGeod::fromDeg( c.getLongitudeDeg() - span/2.0, c.getLatitudeDeg() - SG_HALF_BUCKET_SPAN ); + max = SGGeod::fromDeg( c.getLongitudeDeg() + span/2.0, c.getLatitudeDeg() + SG_HALF_BUCKET_SPAN ); + } else if ( c.getLatitudeDeg() < -89.0) { + min = SGGeod::fromDeg( -90.0, -180.0 ); + max = SGGeod::fromDeg( -89.0, 180.0 ); + } else if ( c.getLatitudeDeg() >= 89.0) { + min = SGGeod::fromDeg( 89.0, -180.0 ); + max = SGGeod::fromDeg( 90.0, 180.0 ); + } else { + SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << c.getLatitudeDeg() ); + } + + SG_LOG( SG_GENERAL, SG_DEBUG, " (" << min << ") (" << max << ")" ); + + // set up clipping tile + base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), min.getLatitudeDeg()) ); + base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), min.getLatitudeDeg()) ); + base.AddNode( 0, SGGeod::fromDeg( max.getLongitudeDeg(), max.getLatitudeDeg()) ); + base.AddNode( 0, SGGeod::fromDeg( min.getLongitudeDeg(), max.getLatitudeDeg()) ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "shape contours = " << subject.Contours() ); + for ( unsigned int ii = 0; ii < subject.Contours(); ii++ ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << subject.GetContour(ii).GetHole() ); + } + + result = tgPolygon::Intersect( subject, base ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "result contours = " << result.Contours() ); + for ( unsigned int ii = 0; ii < result.Contours(); ii++ ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " hole = " << result.GetContour(ii).GetHole() ); + } + + if ( subject.GetPreserve3D() ) { + result.InheritElevations( subject ); + } + + if ( result.Contours() > 0 ) { + result.SetPreserve3D( subject.GetPreserve3D() ); + result.SetTexParams( subject.GetTexParams() ); + if ( subject.GetTexMethod() == TG_TEX_BY_GEODE ) { + // need to set center latitude for geodetic texturing + result.SetTexMethod( TG_TEX_BY_GEODE, b.get_center_lat() ); + } + result.SetFlag(type); + + lock.lock(); + bp_map[b.gen_index()].push_back( result ); + lock.unlock(); + } +} + +void tgChopper::Add( const tgPolygon& subject, const std::string& type ) +{ + tgRectangle bb; + SGGeod p; + + // bail out immediately if polygon is empty + if ( subject.Contours() == 0 ) + return; + + bb = subject.GetBoundingBox(); + + SG_LOG( SG_GENERAL, SG_DEBUG, " min = " << bb.getMin() << " max = " << bb.getMax() ); + + // find buckets for min, and max points of convex hull. + // note to self: self, you should think about checking for + // polygons that span the date line + SGBucket b_min( bb.getMin() ); + SGBucket b_max( bb.getMax() ); + SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket min = " << b_min ); + SG_LOG( SG_GENERAL, SG_DEBUG, " Bucket max = " << b_max ); + + if ( b_min == b_max ) { + // shape entirely contained in a single bucket, write and bail + Clip( subject, type, b_min ); + return; + } + + SGBucket b_cur; + int dx, dy; + + sgBucketDiff(b_min, b_max, &dx, &dy); + SG_LOG( SG_GENERAL, SG_DEBUG, " polygon spans tile boundaries" ); + SG_LOG( SG_GENERAL, SG_DEBUG, " dx = " << dx << " dy = " << dy ); + + if ( (dx > 2880) || (dy > 1440) ) + throw sg_exception("something is really wrong in split_polygon()!!!!"); + + if ( dy <= 1 ) { + // we are down to at most two rows, write each column and then bail + double min_center_lat = b_min.get_center_lat(); + double min_center_lon = b_min.get_center_lon(); + for ( int j = 0; j <= dy; ++j ) { + for ( int i = 0; i <= dx; ++i ) { + b_cur = sgBucketOffset(min_center_lon, min_center_lat, i, j); + Clip( subject, type, b_cur ); + } + } + return; + } + + // we have two or more rows left, split in half (along a + // horizontal dividing line) and recurse with each half + + // find mid point (integer math) + int mid = (dy + 1) / 2 - 1; + + // determine horizontal clip line + SGBucket b_clip = sgBucketOffset( bb.getMin().getLongitudeDeg(), bb.getMin().getLatitudeDeg(), 0, mid); + double clip_line = b_clip.get_center_lat(); + if ( (clip_line >= -90.0 + SG_HALF_BUCKET_SPAN) + && (clip_line < 90.0 - SG_HALF_BUCKET_SPAN) ) + clip_line += SG_HALF_BUCKET_SPAN; + else if ( clip_line < -89.0 ) + clip_line = -89.0; + else if ( clip_line >= 89.0 ) + clip_line = 90.0; + else { + SG_LOG( SG_GENERAL, SG_ALERT, "Out of range latitude in clip_and_write_poly() = " << clip_line ); + } + + { + // + // Crop bottom area (hopefully by putting this in it's own + // scope we can shorten the life of some really large data + // structures to reduce memory use) + // + + SG_LOG( SG_GENERAL, SG_DEBUG, "Generating bottom half (" << bb.getMin().getLatitudeDeg() << "-" << clip_line << ")" ); + + tgPolygon bottom, bottom_clip; + + bottom.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMin().getLatitudeDeg()) ); + bottom.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMin().getLatitudeDeg()) ); + bottom.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) ); + bottom.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) ); + + bottom_clip = tgPolygon::Intersect( subject, bottom ); + + // the texparam should be constant over each clipped poly. + // when they are reassembled, we want the texture map to + // be seamless + Add( bottom_clip, type ); + } + + { + // + // Crop top area (hopefully by putting this in it's own scope + // we can shorten the life of some really large data + // structures to reduce memory use) + // + + SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half (" << clip_line << "-" << bb.getMax().getLatitudeDeg() << ")" ); + + tgPolygon top, top_clip; + + top.AddNode( 0, SGGeod::fromDeg(-180.0, clip_line) ); + top.AddNode( 0, SGGeod::fromDeg( 180.0, clip_line) ); + top.AddNode( 0, SGGeod::fromDeg( 180.0, bb.getMax().getLatitudeDeg()) ); + top.AddNode( 0, SGGeod::fromDeg(-180.0, bb.getMax().getLatitudeDeg()) ); + + top_clip = tgPolygon::Intersect( subject, top ); + + if ( top_clip.TotalNodes() == subject.TotalNodes() ) { + SG_LOG( SG_GENERAL, SG_DEBUG, "Generating top half - total nodes is the same after clip" << subject.TotalNodes() ); + exit(0); + } + + Add( top_clip, type ); + } +} + +long int tgChopper::GenerateIndex( std::string path ) +{ + std::string index_file = path + "/chop.idx"; + long int index = 0; + + //Open or create the named mutex + boost::interprocess::named_mutex mutex(boost::interprocess::open_or_create, "tgChopper_index2"); + { +// SG_LOG(SG_GENERAL, SG_ALERT, "getting lock"); + boost::interprocess::scoped_lock lock(mutex); +// SG_LOG(SG_GENERAL, SG_ALERT, " - got it"); + + /* first try to read the file */ + FILE *fp = fopen( index_file.c_str(), "r+" ); + if ( fp == NULL ) { + /* doesn't exist - create it */ + fp = fopen( index_file.c_str(), "w" ); + if ( fp == NULL ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Error cannot open Index file " << index_file << " for writing"); + boost::interprocess::named_mutex::remove("tgChopper_index2"); + exit( 0 ); + } + } else { + fread( (void*)&index, sizeof(long int), 1, fp ); +// SG_LOG(SG_GENERAL, SG_ALERT, " SUCCESS READING INDEX FILE - READ INDEX " << index ); + } + + index++; + + rewind( fp ); + fwrite( (void*)&index, sizeof(long int), 1, fp ); + fclose( fp ); + } + + boost::interprocess::named_mutex::remove("tgChopper_index2"); + + return index; +} + +void tgChopper::Save( void ) +{ + // traverse the bucket list + bucket_polys_map_interator it; + char tile_name[16]; + char poly_ext[16]; + + for (it=bp_map.begin(); it != bp_map.end(); it++) { + SGBucket b( (*it).first ); + tgpolygon_list const& polys = (*it).second; + + std::string path = root_path + "/" + b.gen_base_path(); + sprintf( tile_name, "%ld", b.gen_index() ); + + std::string polyfile = path + "/" + tile_name; + + SGPath sgp( polyfile ); + sgp.create_dir( 0755 ); + + long int poly_index = GenerateIndex( path ); + + sprintf( poly_ext, "%ld", poly_index ); + polyfile = polyfile + "." + poly_ext; + + gzFile fp; + if ( (fp = gzopen( polyfile.c_str(), "wb9" )) == NULL ) { + SG_LOG( SG_GENERAL, SG_INFO, "ERROR: opening " << polyfile.c_str() << " for writing!" ); + return; + } + + /* Write polys to the file */ + sgWriteUInt( fp, polys.size() ); + for ( unsigned int i=0; i + +#include "tg_polygon.hxx" + +// for ogr-decode : generate a bunch of polygons, mapped by bucket id +typedef std::map bucket_polys_map; +typedef bucket_polys_map::iterator bucket_polys_map_interator; + +class tgChopper +{ +public: + tgChopper( const std::string& path ) { + root_path = path; + } + + void Add( const tgPolygon& poly, const std::string& type ); + void Save( void ); + +private: + long int GenerateIndex( std::string path ); + void Clip( const tgPolygon& subject, const std::string& type, SGBucket& b ); + void Chop( const tgPolygon& subject, const std::string& type ); + + std::string root_path; + bucket_polys_map bp_map; + SGMutex lock; +}; diff --git a/src/Lib/terragear/tg_contour.cxx b/src/Lib/terragear/tg_contour.cxx new file mode 100644 index 00000000..d422edde --- /dev/null +++ b/src/Lib/terragear/tg_contour.cxx @@ -0,0 +1,855 @@ +#include +#include +#include + +#include "tg_misc.hxx" +#include "tg_accumulator.hxx" +#include "tg_contour.hxx" +#include "tg_polygon.hxx" + +tgContour tgContour::Snap( const tgContour& subject, double snap ) +{ + tgContour result; + SGGeod pt; + + for (unsigned int i = 0; i < subject.GetSize(); i++) { + pt = SGGeod_snap( subject.GetNode(i), snap ); + result.AddNode(pt); + } + result.SetHole( subject.GetHole() ); + + return result; +} + +double tgContour::GetMinimumAngle( void ) const +{ + unsigned int p1_index, p2_index, p3_index; + double angle, min_angle = 2.0 * SGD_PI; + unsigned int size = node_list.size(); + + SG_LOG(SG_GENERAL, SG_DEBUG, " tgContour::GetMinimumAngle() : contour size is " << size ); + + for ( unsigned int i = 0; i < size; i++ ) { + if ( i == 0) { + p1_index = size -1; + } else { + p1_index = i - 1; + } + + p2_index = i; + + if ( i == size ) { + p3_index = 0; + } else { + p3_index = i + 1; + } + + angle = SGGeod_CalculateTheta( node_list[p1_index], node_list[p2_index], node_list[p3_index] ); + if ( angle < min_angle ) { + min_angle = angle; + } + } + + return min_angle; +} + +double tgContour::GetArea( void ) const +{ + double area = 0.0; + SGVec2d a, b; + unsigned int i, j; + + if ( node_list.size() ) { + j = node_list.size() - 1; + for (i=0; i 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; +} + +tgContour tgContour::RemoveDups( const tgContour& subject ) +{ + tgContour result; + + int iters = 0; + bool found; + + for ( unsigned int i = 0; i < subject.GetSize(); i++ ) + { + result.AddNode( subject.GetNode( i ) ); + } + result.SetHole( subject.GetHole() ); + + SG_LOG( SG_GENERAL, SG_DEBUG, "remove contour dups : original contour has " << result.GetSize() << " points" ); + + do + { + SG_LOG( SG_GENERAL, SG_DEBUG, "remove_contour_dups: start new iteration" ); + found = false; + + // Step 1 - find a neighboring duplicate points + SGGeod cur, next; + unsigned int cur_loc, next_loc; + + for ( unsigned int i = 0; i < result.GetSize() && !found; i++ ) { + if (i == result.GetSize() - 1 ) { + cur_loc = i; + cur = result.GetNode(i); + + next_loc = 0; + next = result.GetNode(0); + + SG_LOG( SG_GENERAL, SG_DEBUG, " cur is last point: " << cur << " next is first point: " << next ); + } else { + cur_loc = i; + cur = result.GetNode(i); + + next_loc = i+1; + next = result.GetNode(i+1); + SG_LOG( SG_GENERAL, SG_DEBUG, " cur is: " << cur << " next is : " << next ); + } + + if ( SGGeod_isEqual2D( cur, next ) ) { + // keep the point with higher Z + if ( cur.getElevationM() < next.getElevationM() ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_dups: erasing " << cur ); + result.RemoveNodeAt( cur_loc ); + // result.RemoveNodeAt( result.begin()+cur ); + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_dups: erasing " << next ); + result.RemoveNodeAt( next_loc ); + } + found = true; + } + } + + iters++; + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_dups : after " << iters << " iterations, contour has " << result.GetSize() << " points" ); + } while( found ); + + return result; + +} + +tgContour tgContour::SplitLongEdges( const tgContour& subject, double max_len ) +{ + SGGeod p0, p1; + double dist; + tgContour result; + + for ( unsigned i = 0; i < subject.GetSize() - 1; i++ ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "point = " << i); + + p0 = subject.GetNode( i ); + p1 = subject.GetNode( i + 1 ); + + SG_LOG(SG_GENERAL, SG_DEBUG, " " << p0 << " - " << p1); + + if ( fabs(p0.getLatitudeDeg()) < (90.0 - SG_EPSILON) || + fabs(p1.getLatitudeDeg()) < (90.0 - SG_EPSILON) ) + { + dist = SGGeodesy::distanceM( p0, p1 ); + SG_LOG(SG_GENERAL, SG_DEBUG, "distance = " << dist); + + if ( dist > max_len ) { + unsigned int segments = (int)(dist / max_len) + 1; + SG_LOG(SG_GENERAL, SG_DEBUG, "segments = " << segments); + + double dx = (p1.getLongitudeDeg() - p0.getLongitudeDeg()) / segments; + double dy = (p1.getLatitudeDeg() - p0.getLatitudeDeg()) / segments; + + for ( unsigned int j = 0; j < segments; j++ ) { + SGGeod tmp = SGGeod::fromDeg( p0.getLongitudeDeg() + dx * j, p0.getLatitudeDeg() + dy * j ); + SG_LOG(SG_GENERAL, SG_DEBUG, tmp); + result.AddNode( tmp ); + } + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, p0); + result.AddNode( p0 ); + } + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, p0); + result.AddNode( p0 ); + } + + // end of segment is beginning of next segment + } + + p0 = subject.GetNode( subject.GetSize() - 1 ); + p1 = subject.GetNode( 0 ); + + dist = SGGeodesy::distanceM( p0, p1 ); + SG_LOG(SG_GENERAL, SG_DEBUG, "distance = " << dist); + + if ( dist > max_len ) { + unsigned int segments = (int)(dist / max_len) + 1; + SG_LOG(SG_GENERAL, SG_DEBUG, "segments = " << segments); + + double dx = (p1.getLongitudeDeg() - p0.getLongitudeDeg()) / segments; + double dy = (p1.getLatitudeDeg() - p0.getLatitudeDeg()) / segments; + + for ( unsigned int i = 0; i < segments; i++ ) { + SGGeod tmp = SGGeod::fromDeg( p0.getLongitudeDeg() + dx * i, p0.getLatitudeDeg() + dy * i ); + SG_LOG(SG_GENERAL, SG_DEBUG, tmp); + result.AddNode( tmp ); + } + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, p0); + result.AddNode( p0 ); + } + + // maintain original hole flag setting + result.SetHole( subject.GetHole() ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "split_long_edges() complete"); + + return result; +} + +tgContour tgContour::RemoveSpikes( const tgContour& subject ) +{ + tgContour result; + int iters = 0; + double theta; + bool found; + SGGeod cur, prev, next; + + for ( unsigned int i = 0; i < subject.GetSize(); i++ ) + { + result.AddNode( subject.GetNode( i ) ); + } + + SG_LOG(SG_GENERAL, SG_DEBUG, "remove contour spikes : original contour has " << result.GetSize() << " points" ); + + do + { + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_spikes: start new iteration"); + found = false; + + // Step 1 - find a duplicate point + for ( unsigned int i = 0; i < result.GetSize() && !found; i++ ) { + if (i == 0) { + SG_LOG(SG_GENERAL, SG_DEBUG, " cur is first point: " << i << ": " << result.GetNode(i) ); + cur = result.GetNode( 0 ); + prev = result.GetNode( result.GetSize()-1 ); + next = result.GetNode( 1 ); + } else if ( i == result.GetSize()-1 ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " cur is last point: " << i << ": " << result.GetNode(i) ); + cur = result.GetNode( i ); + prev = result.GetNode( i-1 ); + next = result.GetNode( 0 ); + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, " cur is: " << i << ": " << result.GetNode(i) ); + cur = result.GetNode( i ); + prev = result.GetNode( i-1 ); + next = result.GetNode( i+1 ); + } + + theta = SGMiscd::rad2deg( SGGeod_CalculateTheta(prev, cur, next) ); + + if ( abs(theta) < 0.1 ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_spikes: (theta is " << theta << ") erasing " << i << " prev is " << prev << " cur is " << cur << " next is " << next ); + result.RemoveNodeAt( i ); + found = true; + } + } + + iters++; + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_contour_spikes : after " << iters << " iterations, contour has " << result.GetSize() << " points" ); + } while( found ); + + return result; +} + +ClipperLib::Polygon tgContour::ToClipper( const tgContour& subject ) +{ + ClipperLib::Polygon contour; + + for ( unsigned int i=0; i::infinity(); + double miny = std::numeric_limits::infinity(); + double maxx = -std::numeric_limits::infinity(); + double maxy = -std::numeric_limits::infinity(); + + for (unsigned int i = 0; i < node_list.size(); i++) { + SGGeod pt = GetNode(i); + if ( pt.getLongitudeDeg() < minx ) { minx = pt.getLongitudeDeg(); } + if ( pt.getLongitudeDeg() > maxx ) { maxx = pt.getLongitudeDeg(); } + if ( pt.getLatitudeDeg() < miny ) { miny = pt.getLatitudeDeg(); } + if ( pt.getLatitudeDeg() > maxy ) { maxy = pt.getLatitudeDeg(); } + } + + min = SGGeod::fromDeg( minx, miny ); + max = SGGeod::fromDeg( maxx, maxy ); + + return tgRectangle( min, max ); +} + +tgPolygon tgContour::Union( const tgContour& subject, tgPolygon& clip ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { + all_nodes.add( subject.GetNode(i) ); + } + + ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); + ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); + ClipperLib::Polygons clipper_result; + + ClipperLib::Clipper c; + c.Clear(); + c.AddPolygon(clipper_subject, ClipperLib::ptSubject); + c.AddPolygons(clipper_clip, ClipperLib::ptClip); + c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + return result; +} + +tgPolygon tgContour::Diff( const tgContour& subject, tgPolygon& clip ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { + all_nodes.add( subject.GetNode(i) ); + } + + ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); + ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); + ClipperLib::Polygons clipper_result; + + ClipperLib::Clipper c; + c.Clear(); + c.AddPolygon(clipper_subject, ClipperLib::ptSubject); + c.AddPolygons(clipper_clip, ClipperLib::ptClip); + c.Execute(ClipperLib::ctDifference, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + return result; +} + +static bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, + const std::vector& nodes, SGGeod& 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]; + + 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 = current; + 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]; + + 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 = current; + 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; + + SG_LOG(SG_GENERAL, SG_BULK, " " << p0 << " <==> " << p1 ); + + bool found_extra = FindIntermediateNode( p0, p1, nodes, new_pt, bbEpsilon, errEpsilon ); + + if ( found_extra ) { + AddIntermediateNodes( p0, new_pt, nodes, result, bbEpsilon, errEpsilon ); + + result.AddNode( new_pt ); + SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt); + + AddIntermediateNodes( new_pt, p1, nodes, result, bbEpsilon, errEpsilon ); + } +} + +tgContour tgContour::AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ) +{ + SGGeod p0, p1; + tgContour result; + std::vector& tmp_nodes = nodes.get_list(); + + 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, tmp_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, tmp_nodes, result, SG_EPSILON*10, SG_EPSILON*4 ); + + // maintain original hole flag setting + result.SetHole( subject.GetHole() ); + + return result; +} + +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, 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, 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 ) +{ + SGGeod p0, p1; + SGGeod new_pt; + std::vector tmp_nodes; + + tmp_nodes.push_back( node ); + for ( unsigned int n = 0; n < subject.GetSize()-1; n++ ) { + p0 = subject.GetNode( n ); + p1 = subject.GetNode( n+1 ); + + // add intermediate points + bool found_extra = FindIntermediateNode( p0, p1, tmp_nodes, new_pt, SG_EPSILON*10, SG_EPSILON*4 ); + if ( found_extra ) { + start = p0; + end = p1; + return true; + } + } + + // check last segment + p0 = subject.GetNode( subject.GetSize() - 1 ); + p1 = subject.GetNode( 0 ); + + // add intermediate points + bool found_extra = FindIntermediateNode( p0, p1, tmp_nodes, new_pt, SG_EPSILON*10, SG_EPSILON*4 ); + if ( found_extra ) { + start = p0; + end = p1; + return true; + } + + return false; +} + +tgContour tgContour::Expand( const tgContour& subject, double offset ) +{ + tgPolygon poly; + tgContour result; + + poly.AddContour( subject ); + ClipperLib::Polygons clipper_src, clipper_dst; + + clipper_src = tgPolygon::ToClipper( poly ); + + // convert delta from meters to clipper units + OffsetPolygons( clipper_src, clipper_dst, Dist_ToClipper(offset) ); + + poly = tgPolygon::FromClipper( clipper_dst ); + + if ( poly.Contours() == 1 ) { + result = poly.GetContour( 0 ); + } else { + SG_LOG(SG_GENERAL, SG_INFO, "Expanding contour resulted in more than 1 contour ! "); + exit(0); + } + + return result; +} + +tgpolygon_list tgContour::ExpandToPolygons( const tgContour& subject, double width ) +{ + int turn_dir; + + SGGeod cur_inner; + SGGeod cur_outer; + SGGeod prev_inner; + SGGeod prev_outer; + SGGeod calc_inner; + SGGeod calc_outer; + + double last_end_v = 0.0f; + + tgContour expanded; + tgPolygon segment; + tgAccumulator accum; + tgpolygon_list result; + + // generate poly and texparam lists for each line segment + for (unsigned int i = 0; i < subject.GetSize(); i++) + { + last_end_v = 0.0f; + turn_dir = 0; + + sglog().setLogLevels( SG_ALL, SG_INFO ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "makePolygonsTP: calculating offsets for segment " << i); + + // for each point on the PointsList, generate a quad from + // start to next, offset by 1/2 width from the edge + if (i == 0) + { + // first point on the list - offset heading is 90deg + cur_outer = OffsetPointFirst( subject.GetNode(i), subject.GetNode(i+1), -width/2.0f ); + cur_inner = OffsetPointFirst( subject.GetNode(i), subject.GetNode(i+1), width/2.0f ); + } + else if (i == subject.GetSize()-1) + { + // last point on the list - offset heading is 90deg + cur_outer = OffsetPointLast( subject.GetNode(i-1), subject.GetNode(i), -width/2.0f ); + cur_inner = OffsetPointLast( subject.GetNode(i-1), subject.GetNode(i), width/2.0f ); + } + else + { + // middle section + cur_outer = OffsetPointMiddle( subject.GetNode(i-1), subject.GetNode(i), subject.GetNode(i+1), -width/2.0f, turn_dir ); + cur_inner = OffsetPointMiddle( subject.GetNode(i-1), subject.GetNode(i), subject.GetNode(i+1), width/2.0f, turn_dir ); + } + + if ( i > 0 ) + { + SGGeod prev_mp = midpoint( prev_outer, prev_inner ); + SGGeod cur_mp = midpoint( cur_outer, cur_inner ); + SGGeod intersect; + double heading; + double dist; + double az2; + + SGGeodesy::inverse( prev_mp, cur_mp, heading, az2, dist ); + + expanded.Erase(); + segment.Erase(); + + expanded.AddNode( prev_inner ); + expanded.AddNode( prev_outer ); + + // we need to extend one of the points so we're sure we don't create adjacent edges + if (turn_dir == 0) + { + // turned right - offset outer + if ( intersection( prev_inner, prev_outer, cur_inner, cur_outer, intersect ) ) + { + // yes - make a triangle with inner edge = 0 + expanded.AddNode( cur_outer ); + cur_inner = prev_inner; + } + else + { + expanded.AddNode( cur_outer ); + expanded.AddNode( cur_inner ); + } + } + else + { + // turned left - offset inner + if ( intersection( prev_inner, prev_outer, cur_inner, cur_outer, intersect ) ) + { + // yes - make a triangle with outer edge = 0 + expanded.AddNode( cur_inner ); + cur_outer = prev_outer; + } + else + { + expanded.AddNode( cur_outer ); + expanded.AddNode( cur_inner ); + } + } + + expanded.SetHole(false); + segment.AddContour(expanded); + segment.SetTexParams( prev_inner, width, 20.0f, heading ); + segment.SetTexLimits( 0, last_end_v, 1, 1 ); + segment.SetTexMethod( TG_TEX_BY_TPS_CLIPU, -1.0, 0.0, 1.0, 0.0 ); + result.push_back( segment ); + + last_end_v = 1.0f - (fmod( (double)(dist - last_end_v), (double)1.0f )); + } + + prev_outer = cur_outer; + prev_inner = cur_inner; + } + + sglog().setLogLevels( SG_ALL, SG_INFO ); + + return result; +} + +void tgContour::SaveToGzFile( gzFile& fp ) const +{ + // Save the nodelist + sgWriteUInt( fp, node_list.size() ); + for (unsigned int i = 0; i < node_list.size(); i++) { + sgWriteGeod( fp, node_list[i] ); + } + + // and the hole flag + sgWriteInt( fp, (int)hole ); +} + +void tgContour::LoadFromGzFile( gzFile& fp ) +{ + unsigned int count; + SGGeod node; + + // Start Clean + Erase(); + + // Load the nodelist + sgReadUInt( fp, &count ); + for (unsigned int i = 0; i < count; i++) { + sgReadGeod( fp, node ); + node_list.push_back( node ); + } + + sgReadInt( fp, (int *)&hole ); +} + +std::ostream& operator<< ( std::ostream& output, const tgContour& subject ) +{ + // Save the data + output << "NumNodes: " << subject.node_list.size() << "\n"; + + for( unsigned int n=0; n +#include +#include + +#include "tg_unique_geod.hxx" +#include "tg_rectangle.hxx" +#include "clipper.hpp" + +/* forward declarations */ +class tgPolygon; +typedef std::vector tgpolygon_list; + +class tgContour +{ +public: + tgContour() { + hole = false; + } + + void Erase() { + node_list.clear(); + } + + void SetHole( bool h ) { + hole = h; + } + bool GetHole( void ) const { + return hole; + } + + unsigned int GetSize( void ) const { + return node_list.size(); + } + + void Resize( int size ) { + node_list.resize( size ); + } + + void AddNode( SGGeod n ) { + node_list.push_back( n ); + } + void SetNode( unsigned int i, SGGeod n ) { + node_list[i] = n; + } + SGGeod GetNode( unsigned int i ) const { + return node_list[i]; + } + SGGeod const& operator[]( int index ) const { + return node_list[index]; + } + + void RemoveNodeAt( unsigned int idx ) { + if ( idx < node_list.size() ) { + node_list.erase( node_list.begin() + idx ); + } + } + void RemoveNodeRange( unsigned int from, unsigned int to ) { + if ( ( from < to ) && ( to < node_list.size() ) ) { + node_list.erase( node_list.begin()+from,node_list.begin()+to ); + } + } + + tgRectangle GetBoundingBox( void ) const; + + double GetMinimumAngle( void ) const; + double GetArea( void ) const; + + 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 tgPolygon Union( const tgContour& subject, tgPolygon& clip ); + static tgPolygon Diff( const tgContour& subject, tgPolygon& clip ); + + static tgContour AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ); + static tgContour AddColinearNodes( const tgContour& subject, std::vector& nodes ); + static bool FindColinearLine( const tgContour& subject, const SGGeod& node, SGGeod& start, SGGeod& end ); + + // conversions + static ClipperLib::Polygon ToClipper( const tgContour& subject ); + static tgContour FromClipper( const ClipperLib::Polygon& subject ); + + static tgContour Expand( const tgContour& subject, double offset ); + static tgpolygon_list ExpandToPolygons( const tgContour& subject, double width ); + + static void ToShapefile( const tgContour& subject, const std::string& datasource, const std::string& layer, const std::string& feature ); + + void SaveToGzFile( gzFile& fp ) const; + void LoadFromGzFile( gzFile& fp ); + + // Friend for output + friend std::ostream& operator<< ( std::ostream&, const tgContour& ); + +private: + std::vector node_list; + bool hole; +}; + +typedef std::vector tgcontour_list; +typedef tgcontour_list::iterator tgcontour_list_iterator; +typedef tgcontour_list::const_iterator const_tgcontour_list_iterator; + +#endif // _TGCONTOUR_HXX \ No newline at end of file diff --git a/src/Lib/terragear/tg_light.hxx b/src/Lib/terragear/tg_light.hxx new file mode 100644 index 00000000..09f7a5fd --- /dev/null +++ b/src/Lib/terragear/tg_light.hxx @@ -0,0 +1,95 @@ +#ifndef _TG_LIGHT_HXX +#define _TG_LIGHT_HXX + +#include + +#include + +class tgLight +{ +public: + SGGeod pos; + SGVec3f norm; +}; + +typedef std::vector tglight_list; +typedef tglight_list::iterator tglight_list_iterator; +typedef tglight_list::const_iterator const_tglight_list_iterator; + +class tgLightContour +{ +public: + unsigned int ContourSize( void ) const { + return lights.size(); + } + + void AddLight( SGGeod p, SGVec3f n ) { + tgLight light; + light.pos = p; + light.norm = n; + lights.push_back(light); + } + + void SetElevation( unsigned int i, double elev ) { + lights[i].pos.setElevationM( elev ); + } + + SGGeod GetNode( unsigned int i ) const { + return lights[i].pos; + } + + SGGeod GetPosition( unsigned int i ) const { + return lights[i].pos; + } + + SGVec3f GetNormal( unsigned int i ) const { + return lights[i].norm; + } + + std::string GetFlag( void ) const { + return flag; + } + void SetFlag( const std::string f ) { + flag = f; + } + + std::string GetType( void ) const { + return type; + } + void SetType( const std::string t ) { + type = t; + } + + std::vector GetPositionList( void ) { + std::vector positions; + + for (unsigned int i=0; i GetNormalList( void ) { + std::vector normals; + + for (unsigned int i=0; i tglightcontour_list; +typedef tglightcontour_list::iterator tglightcontour_list_iterator; +typedef tglightcontour_list::const_iterator const_tglightcontour_list_iterator; + +#endif // _TG_LIGHT_HXX \ No newline at end of file diff --git a/src/Lib/terragear/tg_misc.cxx b/src/Lib/terragear/tg_misc.cxx new file mode 100644 index 00000000..566d8b30 --- /dev/null +++ b/src/Lib/terragear/tg_misc.cxx @@ -0,0 +1,269 @@ +#include + +#include +#include + +#include + +#include "tg_polygon.hxx" +#include "tg_misc.hxx" + +const double isEqual2D_Epsilon = 0.000001; + +#define CLIPPER_FIXEDPT (10000000000000000) +#define CLIPPER_FIXED1M ( 90090) + +SGGeod SGGeod_snap( const SGGeod& in, double grid ) +{ + return SGGeod::fromDegM( grid * SGMisc::round( in.getLongitudeDeg()/grid ), + grid * SGMisc::round( in.getLatitudeDeg() /grid ), + grid * SGMisc::round( in.getElevationM() /grid ) ); +} + +bool SGGeod_isEqual2D( const SGGeod& g0, const SGGeod& g1 ) +{ + return ( (fabs( g0.getLongitudeDeg() - g1.getLongitudeDeg() ) < isEqual2D_Epsilon) && + (fabs( g0.getLatitudeDeg() - g1.getLatitudeDeg() ) < isEqual2D_Epsilon ) ); +} + +SGVec2d SGGeod_ToSGVec2d( const SGGeod& p ) +{ + return SGVec2d( p.getLongitudeDeg(), p.getLatitudeDeg() ); +} + +// Calculate theta of angle (a, b, c) +double SGGeod_CalculateTheta( const SGGeod& p0, const SGGeod& p1, const SGGeod& p2 ) +{ + SGVec2d v0 = SGGeod_ToSGVec2d( p0 ); + SGVec2d v1 = SGGeod_ToSGVec2d( p1 ); + SGVec2d v2 = SGGeod_ToSGVec2d( p2 ); + + return SGVec2d_CalculateTheta( v0, v1, v2 ); +} + +// Calculate theta of angle of v0, v1, and v2 - vectors represent cartestian positions +double SGVec2d_CalculateTheta( const SGVec2d& v0, const SGVec2d& v1, const SGVec2d& v2 ) +{ + SGVec2d u, v; + double udist, vdist, uv_dot; + + // u . v = ||u|| * ||v|| * cos(theta) + + u = v1 - v0; + udist = dist(v1, v2); + + v = v1 - v2; + vdist = dist(v1, v2); + + uv_dot = dot(u, v); + + return acos(uv_dot / (udist * vdist) ); +} + +// calculate theta from two directions +double CalculateTheta( const SGVec3d& dirCur, const SGVec3d& dirNext ) +{ + double dp = dot( dirCur, dirNext ); + + return acos( dp ); +} + +ClipperLib::IntPoint SGGeod_ToClipper( const SGGeod& p ) +{ + ClipperLib::long64 x, y; + + x = (ClipperLib::long64)( p.getLongitudeDeg() * CLIPPER_FIXEDPT ); + y = (ClipperLib::long64)( p.getLatitudeDeg() * CLIPPER_FIXEDPT ); + + return ClipperLib::IntPoint( x, y ); +} + +SGGeod SGGeod_FromClipper( const ClipperLib::IntPoint& p ) +{ + double lon, lat; + + lon = (double)( ((double)p.X) / (double)CLIPPER_FIXEDPT ); + lat = (double)( ((double)p.Y) / (double)CLIPPER_FIXEDPT ); + + return SGGeod::fromDeg( lon, lat ); +} + +double Dist_ToClipper( double dist ) +{ + return ( dist * ( CLIPPER_FIXEDPT / CLIPPER_FIXED1M ) ); +} + +#ifdef _MSC_VER +# define LONG_LONG_MAX LLONG_MAX +# define LONG_LONG_MIN LLONG_MIN +#endif + +tgRectangle BoundingBox_FromClipper( const ClipperLib::Polygons& subject ) +{ + ClipperLib::IntPoint min_pt, max_pt; + SGGeod min, max; + + min_pt.X = min_pt.Y = LONG_LONG_MAX; + max_pt.X = max_pt.Y = LONG_LONG_MIN; + + // for each polygon, we need to check the orientation, to set the hole flag... + for (unsigned int i=0; i max_pt.X ) { + max_pt.X = subject[i][j].X; + } + if ( subject[i][j].Y > max_pt.Y ) { + max_pt.Y = subject[i][j].Y; + } + } + } + + min = SGGeod_FromClipper( min_pt ); + max = SGGeod_FromClipper( max_pt ); + + return tgRectangle( min, max ); +} + +SGGeod OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by, int& turn_dir ) +{ + double courseCur, courseNext, courseAvg, theta; + SGVec3d dirCur, dirNext, dirAvg, cp; + double courseOffset, distOffset; + SGGeod pt; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Find average angle for contour: prev (" << gPrev << "), " + "cur (" << gCur << "), " + "next (" << gNext << ")" ); + + // first, find if the line turns left or right ar src + // for this, take the cross product of the vectors from prev to src, and src to next. + // if the cross product is negetive, we've turned to the left + // if the cross product is positive, we've turned to the right + courseCur = SGGeodesy::courseDeg( gCur, gPrev ); + dirCur = SGVec3d( sin( courseCur*SGD_DEGREES_TO_RADIANS ), cos( courseCur*SGD_DEGREES_TO_RADIANS ), 0.0f ); + + courseNext = SGGeodesy::courseDeg( gCur, gNext ); + dirNext = SGVec3d( sin( courseNext*SGD_DEGREES_TO_RADIANS ), cos( courseNext*SGD_DEGREES_TO_RADIANS ), 0.0f ); + + // Now find the average + dirAvg = normalize( dirCur + dirNext ); + courseAvg = SGMiscd::rad2deg( atan( dirAvg.x()/dirAvg.y() ) ); + if (courseAvg < 0) { + courseAvg += 180.0f; + } + + // check the turn direction + cp = cross( dirCur, dirNext ); + theta = SGMiscd::rad2deg(CalculateTheta( dirCur, dirNext ) ); + + if ( (abs(theta - 180.0) < 0.1) || (abs(theta) < 0.1) || (isnan(theta)) ) { + // straight line blows up math - offset 90 degree and dist is as given + courseOffset = SGMiscd::normalizePeriodic(0, 360, courseNext-90.0); + distOffset = offset_by; + } else { + // calculate correct distance for the offset point + if (cp.z() < 0.0f) { + courseOffset = SGMiscd::normalizePeriodic(0, 360, courseAvg+180); + turn_dir = 0; + } else { + courseOffset = SGMiscd::normalizePeriodic(0, 360, courseAvg); + turn_dir = 1; + } + distOffset = (offset_by)/sin(SGMiscd::deg2rad(courseNext-courseOffset)); + } + + // calculate the point from cur + pt = SGGeodesy::direct(gCur, courseOffset, distOffset); + SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << courseOffset << " distance is " << distOffset << " point is (" << pt.getLatitudeDeg() << "," << pt.getLongitudeDeg() << ")" ); + + return pt; +} + +SGGeod OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by ) +{ + int unused; + return OffsetPointMiddle( gPrev, gCur, gNext, offset_by, unused ); +} + +SGGeod OffsetPointFirst( const SGGeod& cur, const SGGeod& next, double offset_by ) +{ + double courseOffset; + SGGeod pt; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at Start : cur (" << cur << "), " + "next (" << next << ")" ); + + // find the offset angle + courseOffset = SGGeodesy::courseDeg( cur, next ) - 90; + courseOffset = SGMiscd::normalizePeriodic(0, 360, courseOffset); + + // calculate the point from cur + pt = SGGeodesy::direct( cur, courseOffset, offset_by ); + SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << courseOffset << " distance is " << offset_by << " point is (" << pt.getLatitudeDeg() << "," << pt.getLongitudeDeg() << ")" ); + + return pt; +} + +SGGeod OffsetPointLast( const SGGeod& prev, const SGGeod& cur, double offset_by ) +{ + double courseOffset; + SGGeod pt; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at End : prev (" << prev << "), " + "cur (" << cur << ")" ); + + // find the offset angle + courseOffset = SGGeodesy::courseDeg( prev, cur ) - 90; + courseOffset = SGMiscd::normalizePeriodic(0, 360, courseOffset); + + // calculate the point from cur + pt = SGGeodesy::direct( cur, courseOffset, offset_by ); + SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << courseOffset << " distance is " << offset_by << " point is (" << pt.getLatitudeDeg() << "," << pt.getLongitudeDeg() << ")" ); + + return pt; +} + +SGGeod midpoint( const SGGeod& p0, const SGGeod& p1 ) +{ + return SGGeod::fromDegM( (p0.getLongitudeDeg() + p1.getLongitudeDeg()) / 2, + (p0.getLatitudeDeg() + p1.getLatitudeDeg()) / 2, + (p0.getElevationM() + p1.getElevationM()) / 2 ); +} + +bool intersection(const SGGeod &p0, const SGGeod &p1, const SGGeod& p2, const SGGeod& p3, SGGeod& intersection) +{ + typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; + typedef Kernel::Point_2 Point_2; + typedef CGAL::Segment_2 Segment_2; + + Point_2 a1( p0.getLongitudeDeg(), p0.getLatitudeDeg() ); + Point_2 b1( p1.getLongitudeDeg(), p1.getLatitudeDeg() ); + Point_2 a2( p2.getLongitudeDeg(), p2.getLatitudeDeg() ); + Point_2 b2( p3.getLongitudeDeg(), p3.getLatitudeDeg() ); + + Segment_2 seg1( a1, b1 ); + Segment_2 seg2( a2, b2 ); + + CGAL::Object result = CGAL::intersection(seg1, seg2); + if (const CGAL::Point_2 *ipoint = CGAL::object_cast >(&result)) { + // handle the point intersection case with *ipoint. + return true; + } else { + if (const CGAL::Segment_2 *iseg = CGAL::object_cast >(&result)) { + // handle the segment intersection case with *iseg. + return false; + } else { + // handle the no intersection case. + return false; + } + } +} \ No newline at end of file diff --git a/src/Lib/terragear/tg_misc.hxx b/src/Lib/terragear/tg_misc.hxx new file mode 100644 index 00000000..7af9349e --- /dev/null +++ b/src/Lib/terragear/tg_misc.hxx @@ -0,0 +1,33 @@ +#include + +#include "clipper.hpp" +#include "tg_rectangle.hxx" + +// SGGeod Cleanup +SGGeod SGGeod_snap( const SGGeod& in, double grid ); + +// SGGeod Conversion +SGVec2d SGGeod_ToSGVec2d( const SGGeod& p ); +ClipperLib::IntPoint SGGeod_ToClipper( const SGGeod& p ); +SGGeod SGGeod_FromClipper( const ClipperLib::IntPoint& p ); + +// SGGeod Equivelence test +bool SGGeod_isEqual2D( const SGGeod& g0, const SGGeod& g1 ); + +// SGGeod Angle calculation +double SGGeod_CalculateTheta( const SGGeod& p0, const SGGeod& p1, const SGGeod& p2 ); + +// SGVec2d Angle calculation +double SGVec2d_CalculateTheta( const SGVec2d& v0, const SGVec2d& v1, const SGVec2d& v2); + +// Angle calculation from 2 directions +double CalculateTheta( const SGVec3d& dirCur, const SGVec3d& dirNext ); + +// Clipper Conversion +double Dist_ToClipper( double dist ); + +// should be in rectangle +tgRectangle BoundingBox_FromClipper( const ClipperLib::Polygons& subject ); + +bool intersection(const SGGeod &p0, const SGGeod &p1, const SGGeod& p2, const SGGeod& p3, SGGeod& intersection); + diff --git a/src/Lib/Polygon/tg_nodes.cxx b/src/Lib/terragear/tg_nodes.cxx similarity index 100% rename from src/Lib/Polygon/tg_nodes.cxx rename to src/Lib/terragear/tg_nodes.cxx diff --git a/src/Lib/Polygon/tg_nodes.hxx b/src/Lib/terragear/tg_nodes.hxx similarity index 99% rename from src/Lib/Polygon/tg_nodes.hxx rename to src/Lib/terragear/tg_nodes.hxx index 3578dce8..fc1e9089 100644 --- a/src/Lib/Polygon/tg_nodes.hxx +++ b/src/Lib/terragear/tg_nodes.hxx @@ -39,7 +39,7 @@ 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 diff --git a/src/Lib/terragear/tg_polygon.cxx b/src/Lib/terragear/tg_polygon.cxx new file mode 100644 index 00000000..644cddd6 --- /dev/null +++ b/src/Lib/terragear/tg_polygon.cxx @@ -0,0 +1,519 @@ +// polygon.cxx -- polygon (with holes) management class +// +// Written by Curtis Olson, started March 1999. +// +// Copyright (C) 1999 Curtis L. Olson - http://www.flightgear.org/~curt +// +// This program is free software; you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. +// + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "tg_misc.hxx" +#include "tg_polygon.hxx" + +// tgPolygon static functions +unsigned int tgPolygon::TotalNodes( void ) const +{ + unsigned int total_nodes = 0; + + for (unsigned int c = 0; c < contours.size(); c++) { + total_nodes += contours[c].GetSize(); + } + + return total_nodes; +} + +tgPolygon tgPolygon::Expand( const tgPolygon& subject, double offset ) +{ + ClipperLib::Polygons clipper_src, clipper_dst; + clipper_src = tgPolygon::ToClipper( subject ); + tgPolygon result; + + // convert delta from meters to clipper units + OffsetPolygons( clipper_src, clipper_dst, Dist_ToClipper(offset) ); + + result = tgPolygon::FromClipper( clipper_dst ); + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + return result; +} + +tgPolygon tgPolygon::Expand( const SGGeod& subject, double offset ) +{ + tgPolygon result; + tgContour contour; + SGGeod pt; + + pt = SGGeodesy::direct( subject, 90, offset/2.0 ); + double dlon = pt.getLongitudeDeg() - subject.getLongitudeDeg(); + + pt = SGGeodesy::direct( subject, 0, offset/2.0 ); + double dlat = pt.getLatitudeDeg() - subject.getLatitudeDeg(); + + contour.AddNode( SGGeod::fromDeg( subject.getLongitudeDeg() - dlon, subject.getLatitudeDeg() - dlat ) ); + contour.AddNode( SGGeod::fromDeg( subject.getLongitudeDeg() + dlon, subject.getLatitudeDeg() - dlat ) ); + 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; +} + +tgRectangle tgPolygon::GetBoundingBox( void ) const +{ + SGGeod min, max; + + double minx = std::numeric_limits::infinity(); + double miny = std::numeric_limits::infinity(); + double maxx = -std::numeric_limits::infinity(); + double maxy = -std::numeric_limits::infinity(); + + for ( unsigned int i = 0; i < Contours(); i++ ) { + for (unsigned int j = 0; j < ContourSize(i); j++) { + SGGeod pt = GetNode(i,j); + if ( pt.getLongitudeDeg() < minx ) { minx = pt.getLongitudeDeg(); } + if ( pt.getLongitudeDeg() > maxx ) { maxx = pt.getLongitudeDeg(); } + if ( pt.getLatitudeDeg() < miny ) { miny = pt.getLatitudeDeg(); } + if ( pt.getLatitudeDeg() > maxy ) { maxy = pt.getLatitudeDeg(); } + } + } + + min = SGGeod::fromDeg( minx, miny ); + max = SGGeod::fromDeg( maxx, maxy ); + + return tgRectangle( min, max ); +} + +tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, std::vector& nodes ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + result.AddContour( tgContour::AddColinearNodes( subject.GetContour(c), nodes ) ); + } + + return result; +} + +tgPolygon tgPolygon::AddColinearNodes( const tgPolygon& subject, UniqueSGGeodSet& nodes ) +{ + return AddColinearNodes( subject, nodes.get_list() ); +} + +// 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 ) +{ + bool found = false; + + for ( unsigned int c = 0; c < subject.Contours() && !found; c++ ) { + found = tgContour::FindColinearLine( subject.GetContour(c), node, start, end ); + } + + return found; +} + +SGGeod InterpolateElevation( const SGGeod& dst_node, const SGGeod& start, const SGGeod& end ) +{ + double total_dist = SGGeodesy::distanceM( start, end ); + double inter_dist = SGGeodesy::distanceM( start, dst_node ); + double delta = inter_dist/total_dist; + + double dest_elevation = start.getElevationM() + (delta * ( end.getElevationM() - start.getElevationM() )); + + return SGGeod::fromDegM( dst_node.getLongitudeDeg(), dst_node.getLatitudeDeg(), dest_elevation ); +} + + +void tgPolygon::InheritElevations( const tgPolygon& source ) +{ + UniqueSGGeodSet src_nodes; + + // build a list of points from the source polygon + for ( unsigned int i = 0; i < source.Contours(); ++i ) { + for ( unsigned int j = 0; j < source.ContourSize(i); ++j ) { + src_nodes.add( source.GetNode( i, j ) ); + } + } + + // traverse the dest polygon and build a mirror image but with + // elevations from the source polygon + for ( unsigned int i = 0; i < contours.size(); ++i ) { + for ( unsigned int j = 0; j < contours[i].GetSize(); ++j ) { + SGGeod dst_node = GetNode(i,j); + int index = src_nodes.find( dst_node ); + if ( index >= 0 ) { + SetNode( i, j, src_nodes.get_list()[index] ); + } else { + /* node not is source - we need to find the two points to interpolate from */ + SGGeod start, end, result; + if ( FindColinearLine( source, dst_node, start, end ) ) { + dst_node = InterpolateElevation( dst_node, start, end ); + SetNode( i, j, dst_node ); + } + } + } + } +} + +void tgPolygon::Texture( void ) +{ + SGGeod p; + SGVec2f t; + double x, y; + float tx, ty; + + SG_LOG(SG_GENERAL, SG_DEBUG, "Texture Poly with material " << material << " method " << tp.method << " tpref " << tp.ref << " heading " << tp.heading ); + + switch( tp.method ) { + case TG_TEX_BY_GEODE: + { + // The Simgear General texture coordinate routine takes a fan. + // Simgear could probably use a new function that just takes a Geod vector + // For now, just create an identity fan... + std::vector< int > node_idxs; + for (int i = 0; i < 3; i++) { + node_idxs.push_back(i); + } + + for ( unsigned int i = 0; i < triangles.size(); i++ ) { + std::vector< SGVec2f > tc_list; + std::vector< SGGeod > nodes; + + nodes = triangles[i].GetNodeList(); + tc_list = sgCalcTexCoords( tp.center_lat, nodes, node_idxs ); + triangles[i].SetTexCoordList( tc_list ); + } + } + break; + + case TG_TEX_BY_TPS_NOCLIP: + case TG_TEX_BY_TPS_CLIPU: + case TG_TEX_BY_TPS_CLIPV: + case TG_TEX_BY_TPS_CLIPUV: + { + for ( unsigned int i = 0; i < triangles.size(); i++ ) { + for ( unsigned int j = 0; j < 3; j++ ) { + p = triangles[i].GetNode( j ); + SG_LOG(SG_GENERAL, SG_DEBUG, "point = " << p); + + // + // 1. Calculate distance and bearing from the center of + // the poly + // + + // given alt, lat1, lon1, lat2, lon2, calculate starting + // and ending az1, az2 and distance (s). Lat, lon, and + // azimuth are in degrees. distance in meters + double az1, az2, dist; + SGGeodesy::inverse( tp.ref, p, az1, az2, dist ); + SG_LOG(SG_GENERAL, SG_DEBUG, "basic course = " << az2); + + // + // 2. Rotate this back into a coordinate system where Y + // runs the length of the poly and X runs crossways. + // + + double course = SGMiscd::normalizePeriodic(0, 360, az2 - tp.heading); + SG_LOG( SG_GENERAL, SG_DEBUG," course = " << course << " dist = " << dist ); + + // + // 3. Convert from polar to cartesian coordinates + // + + x = sin( course * SGD_DEGREES_TO_RADIANS ) * dist; + y = cos( course * SGD_DEGREES_TO_RADIANS ) * dist; + SG_LOG(SG_GENERAL, SG_DEBUG, " x = " << x << " y = " << y); + + // + // 4. Map x, y point into texture coordinates + // + float tmp; + + tmp = (float)x / (float)tp.width; + tx = tmp * (float)(tp.maxu - tp.minu) + (float)tp.minu; + SG_LOG(SG_GENERAL, SG_DEBUG, " (" << tx << ")"); + + // clip u? + if ( (tp.method == TG_TEX_BY_TPS_CLIPU) || (tp.method == TG_TEX_BY_TPS_CLIPUV) ) { + if ( tx < (float)tp.min_clipu ) { tx = (float)tp.min_clipu; } + if ( tx > (float)tp.max_clipu ) { tx = (float)tp.max_clipu; } + } + + tmp = (float)y / (float)tp.length; + ty = tmp * (float)(tp.maxv - tp.minv) + (float)tp.minv; + SG_LOG(SG_GENERAL, SG_DEBUG, " (" << ty << ")"); + + // clip v? + if ( (tp.method == TG_TEX_BY_TPS_CLIPV) || (tp.method == TG_TEX_BY_TPS_CLIPUV) ) { + if ( ty < (float)tp.min_clipv ) { ty = (float)tp.min_clipv; } + if ( ty > (float)tp.max_clipv ) { ty = (float)tp.max_clipv; } + } + + t = SGVec2f( tx, ty ); + SG_LOG(SG_GENERAL, SG_DEBUG, " (" << tx << ", " << ty << ")"); + + triangles[i].SetTexCoord( j, t ); + } + } + } + break; + } +} + +void tgPolygon::SaveToGzFile( gzFile& fp ) const +{ + // Save the contours + sgWriteUInt( fp, contours.size() ); + for (unsigned int i = 0; i < contours.size(); i++) { + contours[i].SaveToGzFile( fp ); + } + + // Save the triangles + sgWriteUInt( fp, triangles.size() ); + for (unsigned int i = 0; i < triangles.size(); i++) { + triangles[i].SaveToGzFile( fp ); + } + + // Save the tex params + tp.SaveToGzFile( fp ); + + // and the rest + sgWriteString( fp, material.c_str() ); + sgWriteString( fp, flag.c_str() ); + sgWriteInt( fp, (int)preserve3d ); +} + +void tgPolygon::LoadFromGzFile( gzFile& fp ) +{ + unsigned int count; + tgContour contour; + tgTriangle triangle; + char *strbuff; + + // Start clean + Erase(); + + // Load the contours + sgReadUInt( fp, &count ); + for (unsigned int i = 0; i < count; i++) { + contour.LoadFromGzFile( fp ); + AddContour(contour); + } + + // load the triangles + sgReadUInt( fp, &count ); + for (unsigned int i = 0; i < count; i++) { + triangle.LoadFromGzFile( fp ); + AddTriangle(triangle); + } + + // Load the tex params + tp.LoadFromGzFile( fp ); + + // and the rest + sgReadString( fp, &strbuff ); + if ( strbuff ) { + material = strbuff; + delete strbuff; + } + + sgReadString( fp, &strbuff ); + if ( strbuff ) { + flag = strbuff; + delete strbuff; + } + + sgReadInt( fp, (int *)&preserve3d ); +} + +// Friends for serialization +std::ostream& operator<< ( std::ostream& output, const tgPolygon& subject ) +{ + // Save the data + output << "NumContours: " << subject.contours.size() << "\n"; + + for( unsigned int c=0; c #include #include #include -#include - -#include "tg_unique_geod.hxx" +#include "tg_rectangle.hxx" +#include "tg_contour.hxx" // utilities - belong is simgear? double CalculateTheta( const SGVec3d& dirCur, const SGVec3d& dirNext, const SGVec3d& cp ); @@ -89,98 +86,6 @@ typedef std::vector tgpolygon_list; typedef tgpolygon_list::iterator tgpolygon_list_iterator; typedef tgpolygon_list::const_iterator const_tgpolygon_list_iterator; -class tgContour -{ -public: - tgContour() { - hole = false; - } - - void Erase() { - node_list.clear(); - } - - void SetHole( bool h ) { - hole = h; - } - bool GetHole( void ) const { - return hole; - } - - unsigned int GetSize( void ) const { - return node_list.size(); - } - - void Resize( int size ) { - node_list.resize( size ); - } - - void AddNode( SGGeod n ) { - node_list.push_back( n ); - } - void SetNode( unsigned int i, SGGeod n ) { - node_list[i] = n; - } - SGGeod GetNode( unsigned int i ) const { - return node_list[i]; - } - SGGeod const& operator[]( int index ) const { - return node_list[index]; - } - - void RemoveNodeAt( unsigned int idx ) { - if ( idx < node_list.size() ) { - node_list.erase( node_list.begin() + idx ); - } - } - void RemoveNodeRange( unsigned int from, unsigned int to ) { - if ( ( from < to ) && ( to < node_list.size() ) ) { - node_list.erase( node_list.begin()+from,node_list.begin()+to ); - } - } - - tgRectangle GetBoundingBox( void ) const; - - double GetMinimumAngle( void ) const; - double GetArea( void ) const; - - 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 tgPolygon Union( const tgContour& subject, tgPolygon& clip ); - static tgPolygon Diff( const tgContour& subject, tgPolygon& clip ); - - static tgContour AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ); - static tgContour AddColinearNodes( const tgContour& subject, std::vector& nodes ); - static bool FindColinearLine( const tgContour& subject, const SGGeod& node, SGGeod& start, SGGeod& end ); - - // conversions - static ClipperLib::Polygon ToClipper( const tgContour& subject ); - static tgContour FromClipper( const ClipperLib::Polygon& subject ); - - static tgContour Expand( const tgContour& subject, double offset ); - static tgpolygon_list ExpandToPolygons( const tgContour& subject, double width ); - - static void ToShapefile( const tgContour& subject, const std::string& datasource, const std::string& layer, const std::string& feature ); - - void SaveToGzFile( gzFile& fp ) const; - void LoadFromGzFile( gzFile& fp ); - - // Friend for output - friend std::ostream& operator<< ( std::ostream&, const tgContour& ); - -private: - std::vector node_list; - bool hole; -}; - -typedef std::vector tgcontour_list; -typedef tgcontour_list::iterator tgcontour_list_iterator; -typedef tgcontour_list::const_iterator const_tgcontour_list_iterator; - class tgTriangle { public: @@ -410,7 +315,7 @@ public: void SetPreserve3D( bool p ) { preserve3d = p; } - + unsigned int GetId( void ) const { return id; } @@ -418,6 +323,8 @@ public: id = i; } + + // Texturing void SetTexParams( const SGGeod& ref, double width, double length, double heading ) { tp.ref = ref; tp.width = width; @@ -455,62 +362,50 @@ public: tgTexMethod GetTexMethod( void ) const { return tp.method; } + void Texture( void ); + // Tesselation void Tesselate( void ); void Tesselate( const std::vector& extra ); - void Texture( void ); - // Boolean operations - - // TODO : Both should be constant - // what we really need is multiple accumulators - // init_accumulator should return a handle... - static bool ChopIdxInit( const std::string& path ); - - static void SetDump( bool dmp ); - static tgPolygon Union( const tgContour& subject, tgPolygon& clip ); + static void SetClipperDump( bool dmp ); static tgPolygon Union( const tgPolygon& subject, tgPolygon& clip ); static tgPolygon Union( const tgpolygon_list& polys ); static tgPolygon Diff( const tgPolygon& subject, tgPolygon& clip ); static tgPolygon Intersect( const tgPolygon& subject, const tgPolygon& clip ); - // Conversions - static ClipperLib::Polygons ToClipper( const tgPolygon& subject ); - static tgPolygon FromClipper( const ClipperLib::Polygons& subject ); - - static tgPolygon FromOGR( const OGRPolygon* subject ); - - // other operations + // cleanup operations static tgPolygon Snap( const tgPolygon& subject, double snap ); static tgPolygon StripHoles( const tgPolygon& subject ); static tgPolygon SplitLongEdges( const tgPolygon& subject, double dist ); - static tgPolygon RemoveCycles( const tgPolygon& subject ); static tgPolygon RemoveDups( const tgPolygon& subject ); static tgPolygon RemoveBadContours( const tgPolygon& subject ); static tgPolygon Simplify( const tgPolygon& subject ); static tgPolygon RemoveTinyContours( const tgPolygon& subject ); static tgPolygon RemoveSpikes( const tgPolygon& subject ); + static void RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers ); + static tgcontour_list MergeSlivers( tgpolygon_list& subjects, tgcontour_list& slivers ); static tgPolygon Expand( const tgPolygon& subject, double offset ); static tgPolygon Expand( const SGGeod& subject, double offset ); - + + // Conversions + static ClipperLib::Polygons ToClipper( const tgPolygon& subject ); + static tgPolygon FromClipper( const ClipperLib::Polygons& subject ); + static void ToShapefile( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& feature ); - static void Tesselate( const tgPolygon& subject ); - + // T-Junctions and segment search static tgPolygon AddColinearNodes( const tgPolygon& subject, UniqueSGGeodSet& nodes ); static tgPolygon AddColinearNodes( const tgPolygon& subject, std::vector& nodes ); static bool FindColinearLine( const tgPolygon& subject, SGGeod& node, SGGeod& start, SGGeod& end ); - static void RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers ); - static tgcontour_list MergeSlivers( tgpolygon_list& subjects, tgcontour_list& slivers ); - + // IO void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); - // Friend for output friend std::ostream& operator<< ( std::ostream&, const tgPolygon& ); private: @@ -524,143 +419,4 @@ private: tgTexParams tp; }; -// for ogr-decode : generate a bunch of polygons, mapped by bucket id -typedef std::map bucket_polys_map; -typedef bucket_polys_map::iterator bucket_polys_map_interator; - -class tgChopper -{ -public: - tgChopper( const std::string& path ) { - root_path = path; - } - - void Add( const tgPolygon& poly, const std::string& type ); - void Save( void ); - -private: - long int GenerateIndex( std::string path ); - void Clip( const tgPolygon& subject, const std::string& type, SGBucket& b ); - void Chop( const tgPolygon& subject, const std::string& type ); - - std::string root_path; - bucket_polys_map bp_map; - SGMutex lock; -}; - -class tgLight -{ -public: - SGGeod pos; - SGVec3f norm; -}; - -typedef std::vector tglight_list; -typedef tglight_list::iterator tglight_list_iterator; -typedef tglight_list::const_iterator const_tglight_list_iterator; - -class tgLightContour -{ -public: - unsigned int ContourSize( void ) const { - return lights.size(); - } - - void AddLight( SGGeod p, SGVec3f n ) { - tgLight light; - light.pos = p; - light.norm = n; - lights.push_back(light); - } - - void SetElevation( unsigned int i, double elev ) { - lights[i].pos.setElevationM( elev ); - } - - SGGeod GetNode( unsigned int i ) const { - return lights[i].pos; - } - - SGGeod GetPosition( unsigned int i ) const { - return lights[i].pos; - } - - SGVec3f GetNormal( unsigned int i ) const { - return lights[i].norm; - } - - std::string GetFlag( void ) const { - return flag; - } - void SetFlag( const std::string f ) { - flag = f; - } - - std::string GetType( void ) const { - return type; - } - void SetType( const std::string t ) { - type = t; - } - - std::vector GetPositionList( void ) { - std::vector positions; - - for (unsigned int i=0; i GetNormalList( void ) { - std::vector normals; - - for (unsigned int i=0; i tglightcontour_list; -typedef tglightcontour_list::iterator tglightcontour_list_iterator; -typedef tglightcontour_list::const_iterator const_tglightcontour_list_iterator; - -typedef std::vector < ClipperLib::Polygons > clipper_polygons_list; - -class tgAccumulator -{ -public: - tgPolygon Diff( const tgContour& subject ); - tgPolygon Diff( const tgPolygon& subject ); - - void Add( const tgContour& subject ); - void Add( const tgPolygon& subject ); - - void ToShapefiles( const std::string& path, const std::string& layer ); - -private: - clipper_polygons_list accum; -}; - -class tgShapefile -{ -public: - static void Init( void ); - static void* OpenDatasource( const char* datasource_name ); - static void* OpenLayer( void* ds_id, const char* layer_name ); -// static void CreateFeature( void* ds_id, void* l_id, const TGPolygon &poly, const char* description ); - static void CloseLayer( void* l_id ); - static void* CloseDatasource( void* ds_id ); -}; - #endif // _POLYGON_HXX diff --git a/src/Lib/terragear/tg_polygon_clean.cxx b/src/Lib/terragear/tg_polygon_clean.cxx new file mode 100644 index 00000000..11785628 --- /dev/null +++ b/src/Lib/terragear/tg_polygon_clean.cxx @@ -0,0 +1,272 @@ +#include + +#include "tg_polygon.hxx" + +tgPolygon tgPolygon::Snap( const tgPolygon& subject, double snap ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for (unsigned int c = 0; c < subject.Contours(); c++) { + result.AddContour( tgContour::Snap( subject.GetContour( c ), snap ) ); + } + + return result; +} + +tgPolygon tgPolygon::RemoveDups( const tgPolygon& subject ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + result.AddContour( tgContour::RemoveDups( subject.GetContour( c ) ) ); + } + + return result; +} + +tgPolygon tgPolygon::RemoveBadContours( const tgPolygon& subject ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + tgContour contour = subject.GetContour(c); + if ( contour.GetSize() >= 3 ) { + /* keeping the contour */ + result.AddContour( contour ); + } + } + + return result; +} + +tgPolygon tgPolygon::RemoveCycles( const tgPolygon& subject ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + result.AddContour( tgContour::RemoveCycles( subject.GetContour( c ) ) ); + } + + return result; +} + +tgPolygon tgPolygon::SplitLongEdges( const tgPolygon& subject, double dist ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned c = 0; c < subject.Contours(); c++ ) + { + result.AddContour( tgContour::SplitLongEdges( subject.GetContour(c), dist ) ); + } + + return result; +} + +tgPolygon tgPolygon::StripHoles( const tgPolygon& subject ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.Contours(); ++i ) { + for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { + all_nodes.add( subject.GetNode(i, j) ); + } + } + + ClipperLib::Polygons clipper_result; + ClipperLib::Clipper c; + c.Clear(); + + for ( unsigned int i = 0; i < subject.Contours(); i++ ) { + tgContour contour = subject.GetContour( i ); + if ( !contour.GetHole() ) { + c.AddPolygon( tgContour::ToClipper( contour ), ClipperLib::ptClip ); + } + } + c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + return result; +} + +tgPolygon tgPolygon::Simplify( const tgPolygon& subject ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.Contours(); ++i ) { + for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { + all_nodes.add( subject.GetNode(i, j) ); + } + } + + ClipperLib::Polygons clipper_poly = tgPolygon::ToClipper( subject ); + SimplifyPolygons( clipper_poly ); + + result = tgPolygon::FromClipper( clipper_poly ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + return result; +} + +tgPolygon tgPolygon::RemoveTinyContours( const tgPolygon& subject ) +{ + double min_area = SG_EPSILON*SG_EPSILON; + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + tgContour contour = subject.GetContour( c ); + double area = contour.GetArea(); + + if ( area >= min_area) { + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_tiny_contours NO - " << c << " area is " << area << " requirement is " << min_area); + result.AddContour( contour ); + } else { + SG_LOG(SG_GENERAL, SG_DEBUG, "remove_tiny_contours " << c << " area is " << area << ": removing"); + } + } + + return result; +} + +tgPolygon tgPolygon::RemoveSpikes( const tgPolygon& subject ) +{ + tgPolygon result; + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + for ( unsigned int c = 0; c < subject.Contours(); c++ ) { + result.AddContour( tgContour::RemoveSpikes( subject.GetContour(c) ) ); + } + + return result; +} + +// Move slivers from in polygon to out polygon. +void tgPolygon::RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers ) +{ +#if 0 + // traverse each contour of the polygon and attempt to identify + // likely slivers + SG_LOG(SG_GENERAL, SG_DEBUG, "tgPolygon::RemoveSlivers()"); + + tgPolygon result; + tgContour contour; + int i; + + double angle_cutoff = 10.0 * SGD_DEGREES_TO_RADIANS; + double area_cutoff = 0.000000001; + double min_angle; + double area; + + // process contours in reverse order so deleting a contour doesn't + // foul up our sequence + for ( i = subject.Contours() - 1; i >= 0; --i ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "contour " << i ); + + contour = subject.GetContour(i); + + SG_LOG(SG_GENERAL, SG_DEBUG, " calc min angle for contour " << i); + min_angle = contour.GetMinimumAngle(); + SG_LOG(SG_GENERAL, SG_DEBUG, " min_angle (rad) = " << min_angle ); + + area = contour.GetArea(); + SG_LOG(SG_GENERAL, SG_DEBUG, " area = " << area ); + + if ( ((min_angle < angle_cutoff) && (area < area_cutoff)) || + ( area < area_cutoff / 10.0) ) + { + if ((min_angle < angle_cutoff) && (area < area_cutoff)) + { + SG_LOG(SG_GENERAL, SG_DEBUG, " WE THINK IT'S A SLIVER! - min angle < 10 deg, and area < 10 sq meters"); + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, " WE THINK IT'S A SLIVER! - min angle > 10 deg, but area < 1 sq meters"); + } + + // Remove the sliver from source + subject.DeleteContourAt( i ); + + // And add it to the slive list if it isn't a hole + if ( !contour.GetHole() ) { + // move sliver contour to sliver list + SG_LOG(SG_GENERAL, SG_DEBUG, " Found SLIVER!"); + + slivers.push_back( contour ); + } + } + } +#endif +} + +tgcontour_list tgPolygon::MergeSlivers( tgpolygon_list& polys, tgcontour_list& sliver_list ) { + tgPolygon poly, result; + tgContour sliver; + tgContour contour; + tgcontour_list unmerged; + unsigned int original_contours, result_contours; + bool done; + + for ( unsigned int i = 0; i < sliver_list.size(); i++ ) { + sliver = sliver_list[i]; + SG_LOG(SG_GENERAL, SG_DEBUG, "Merging sliver = " << i ); + + sliver.SetHole( false ); + + done = false; + + // try to merge the slivers with the list of clipped polys + for ( unsigned int j = 0; j < polys.size() && !done; j++ ) { + poly = polys[j]; + original_contours = poly.Contours(); + result = tgContour::Union( sliver, poly ); + result_contours = result.Contours(); + + if ( original_contours == result_contours ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " FOUND a poly to merge the sliver with"); + result.SetMaterial( polys[j].GetMaterial() ); + result.SetTexParams( polys[j].GetTexParams() ); + polys[j] = result; + done = true; + } + } + + if ( !done ) { + SG_LOG(SG_GENERAL, SG_DEBUG, "couldn't merge sliver " << i ); + unmerged.push_back( sliver ); + } + } + + return unmerged; +} + diff --git a/src/Lib/terragear/tg_polygon_clip.cxx b/src/Lib/terragear/tg_polygon_clip.cxx new file mode 100644 index 00000000..04b996be --- /dev/null +++ b/src/Lib/terragear/tg_polygon_clip.cxx @@ -0,0 +1,175 @@ +#include +#include + +#include + +#include "tg_polygon.hxx" + +static bool clipper_dump = false; +void tgPolygon::SetClipperDump( bool dmp ) +{ + clipper_dump = dmp; +} + +tgPolygon tgPolygon::Union( const tgPolygon& subject, tgPolygon& clip ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + std::ofstream dmpfile; + + /* before union - gather all nodes */ + for ( unsigned int i = 0; i < subject.Contours(); ++i ) { + for ( unsigned int j = 0; j < subject.ContourSize( i ); ++j ) { + all_nodes.add( subject.GetNode(i, j) ); + } + } + + ClipperLib::Polygons clipper_subject = tgPolygon::ToClipper( subject ); + ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); + ClipperLib::Polygons clipper_result; + + if ( clipper_dump ) { + dmpfile.open ("subject.txt"); + dmpfile << clipper_subject; + dmpfile.close(); + + dmpfile.open ("clip.txt"); + dmpfile << clipper_clip; + dmpfile.close(); + } + + ClipperLib::Clipper c; + c.Clear(); + c.AddPolygons(clipper_subject, ClipperLib::ptSubject); + c.AddPolygons(clipper_clip, ClipperLib::ptClip); + c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + + if ( clipper_dump ) { + dmpfile.open ("result.txt"); + dmpfile << clipper_result; + dmpfile.close(); + } + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + result.SetMaterial( subject.GetMaterial() ); + result.SetTexParams( subject.GetTexParams() ); + + return result; +} + +tgPolygon tgPolygon::Union( const tgpolygon_list& polys ) +{ + ClipperLib::Polygons clipper_result; + ClipperLib::Clipper c; + UniqueSGGeodSet all_nodes; + tgPolygon result; + + /* before union - gather all nodes */ + for ( unsigned int i=0; i + +#include +#include +#include +#include +#include + +#include + +#include "tg_polygon.hxx" + +/* determining if a face is within the reulting poly */ +struct FaceInfo2 +{ + FaceInfo2() {} + int nesting_level; + + bool in_domain(){ + return nesting_level%2 == 1; + } +}; + +typedef CGAL::Exact_predicates_inexact_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::Constrained_Delaunay_triangulation_2 CDT; +typedef CDT::Point Point; +typedef CGAL::Polygon_2 Polygon_2; +typedef CGAL::Triangle_2 Triangle_2; + +static void tg_mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list& border ) +{ + if(start->info().nesting_level != -1) { + return; + } + + std::list queue; + queue.push_back(start); + + while( !queue.empty() ){ + CDT::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); + if(n->info().nesting_level == -1) { + if(ct.is_constrained(e)) border.push_back(e); + else queue.push_back(n); + } + } + } + } +} + +//explore set of facets connected with non constrained edges, +//and attribute to each such set a nesting level. +//We start from facets incident to the infinite vertex, with a nesting +//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) +{ + for(CDT::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; + tg_mark_domains(cdt, cdt.infinite_face(), index++, border); + while(! border.empty()) { + CDT::Edge e = border.front(); + border.pop_front(); + CDT::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) +{ + if ( polygon.is_empty() ) return; + + CDT::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); + cdt.insert_constraint(vh,v_prev); + v_prev=vh; + } +} + +void tgPolygon::Tesselate( const std::vector& extra ) +{ + CDT cdt; + + SG_LOG( SG_GENERAL, SG_DEBUG, "Tess with extra" ); + + // Bail right away if polygon is empty + if ( contours.size() != 0 ) { + // First, convert the extra points to cgal Points + std::vector points; + points.reserve(extra.size()); + for (unsigned int n = 0; n < extra.size(); n++) { + points.push_back( Point(extra[n].getLongitudeDeg(), extra[n].getLatitudeDeg() ) ); + } + + cdt.insert(points.begin(), points.end()); + + // then insert each polygon as a constraint into the triangulation + for ( unsigned int c = 0; c < contours.size(); c++ ) { + tgContour contour = contours[c]; + Polygon_2 poly; + + for (unsigned int n = 0; n < contour.GetSize(); n++ ) { + SGGeod node = contour.GetNode(n); + poly.push_back( Point( node.getLongitudeDeg(), node.getLatitudeDeg() ) ); + } + + 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) { + 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() ); + + AddTriangle( p0, p1, p2 ); + + ++count; + } + } + } +} + +void tgPolygon::Tesselate() +{ + CDT cdt; + + SG_LOG( SG_GENERAL, SG_DEBUG, "Tess" ); + + // Bail right away if polygon is empty + if ( contours.size() != 0 ) { + // insert each polygon as a constraint into the triangulation + for ( unsigned int c = 0; c < contours.size(); c++ ) { + tgContour contour = contours[c]; + Polygon_2 poly; + + for (unsigned int n = 0; n < contour.GetSize(); n++ ) { + SGGeod node = contour.GetNode(n); + SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : Adding GEOD " << node); + poly.push_back( Point( node.getLongitudeDeg(), node.getLatitudeDeg() ) ); + } + + 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) { + 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() ); + + AddTriangle( p0, p1, p2 ); + + ++count; + } else { + SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : face not in domain"); + } + } + } + else + { + SG_LOG( SG_GENERAL, SG_DEBUG, "Tess : no contours" ); + } +} \ No newline at end of file diff --git a/src/Lib/Polygon/rectangle.cxx b/src/Lib/terragear/tg_rectangle.cxx similarity index 98% rename from src/Lib/Polygon/rectangle.cxx rename to src/Lib/terragear/tg_rectangle.cxx index b3d135e1..8e03279f 100644 --- a/src/Lib/Polygon/rectangle.cxx +++ b/src/Lib/terragear/tg_rectangle.cxx @@ -4,7 +4,7 @@ // // This file is in the Public Domain and comes with NO WARRANTY OF ANY KIND. -#include "rectangle.hxx" +#include "tg_rectangle.hxx" tgRectangle::tgRectangle() { diff --git a/src/Lib/Polygon/rectangle.hxx b/src/Lib/terragear/tg_rectangle.hxx similarity index 100% rename from src/Lib/Polygon/rectangle.hxx rename to src/Lib/terragear/tg_rectangle.hxx diff --git a/src/Lib/terragear/tg_shapefile.cxx b/src/Lib/terragear/tg_shapefile.cxx new file mode 100644 index 00000000..2122fc74 --- /dev/null +++ b/src/Lib/terragear/tg_shapefile.cxx @@ -0,0 +1,208 @@ +#include + +#include + +#include "tg_misc.hxx" +#include "tg_shapefile.hxx" + +bool tgShapefile::initialized = false; + +void* tgShapefile::OpenDatasource( const char* datasource_name ) +{ + OGRDataSource* datasource; + OGRSFDriver* ogrdriver; + const char* format_name = "ESRI Shapefile"; + + if (!tgShapefile::initialized) { + OGRRegisterAll(); + tgShapefile::initialized = true; + } + + ogrdriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName( format_name ); + if ( !ogrdriver ) { + SG_LOG( SG_GENERAL, SG_ALERT, "Unknown datasource format driver: " << format_name ); + exit(1); + } + + datasource = ogrdriver->Open( datasource_name, TRUE ); + if ( !datasource ) { + datasource = ogrdriver->CreateDataSource( datasource_name, NULL ); + } + + if ( !datasource ) { + SG_LOG( SG_GENERAL, SG_ALERT, "Unable to open or create datasource: " << datasource_name ); + exit(1); + } + + return (void*)datasource; +} + +void* tgShapefile::OpenLayer( void* ds_id, const char* layer_name ) { + OGRDataSource* datasource = ( OGRDataSource * )ds_id; + OGRLayer* layer; + + OGRSpatialReference srs; + srs.SetWellKnownGeogCS("WGS84"); + + layer = datasource->GetLayerByName( layer_name ); + + if ( !layer ) { + layer = datasource->CreateLayer( layer_name, &srs, wkbPolygon25D, NULL ); + + OGRFieldDefn descriptionField( "ID", OFTString ); + descriptionField.SetWidth( 128 ); + + if( layer->CreateField( &descriptionField ) != OGRERR_NONE ) { + SG_LOG( SG_GENERAL, SG_ALERT, "Creation of field 'Description' failed" ); + } + } + + if ( !layer ) { + SG_LOG( SG_GENERAL, SG_ALERT, "Creation of layer '" << layer_name << "' failed" ); + return NULL; + } + + return (void*)layer; +} + +void* tgShapefile::CloseDatasource( void* ds_id ) +{ + OGRDataSource* datasource = ( OGRDataSource * )ds_id; + OGRDataSource::DestroyDataSource( datasource ); + + return (void *)-1; +} + +void tgShapefile::FromClipper( const ClipperLib::Polygons& subject, const std::string& datasource, const std::string& layer, const std::string& description ) +{ + void* ds_id = tgShapefile::OpenDatasource( datasource.c_str() ); + SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenDatasource returned " << (unsigned long)ds_id); + + OGRLayer* l_id = (OGRLayer *)tgShapefile::OpenLayer( ds_id, layer.c_str() ); + SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenLayer returned " << (unsigned long)l_id); + + OGRPolygon* polygon = new OGRPolygon(); + SG_LOG(SG_GENERAL, SG_DEBUG, "subject has " << subject.size() << " contours "); + + for ( unsigned int i = 0; i < subject.size(); i++ ) { + ClipperLib::Polygon const& contour = subject[i]; + + if (contour.size() < 3) { + SG_LOG(SG_GENERAL, SG_DEBUG, "Polygon with less than 3 points"); + } else { + OGRLinearRing *ring = new OGRLinearRing(); + + // FIXME: Current we ignore the hole-flag and instead assume + // that the first ring is not a hole and the rest + // are holes + for (unsigned int pt = 0; pt < contour.size(); pt++) { + OGRPoint *point = new OGRPoint(); + SGGeod geod = SGGeod_FromClipper( contour[pt] ); + + point->setX( geod.getLongitudeDeg() ); + point->setY( geod.getLatitudeDeg() ); + point->setZ( 0.0 ); + + ring->addPoint( point ); + } + + ring->closeRings(); + polygon->addRingDirectly( ring ); + } + + OGRFeature* feature = new OGRFeature( l_id->GetLayerDefn() ); + + feature->SetField("ID", description.c_str()); + feature->SetGeometry(polygon); + if( l_id->CreateFeature( feature ) != OGRERR_NONE ) + { + SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile"); + } + + OGRFeature::DestroyFeature(feature); + } +} + +void tgShapefile::FromPolygon( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description ) +{ + void* ds_id = tgShapefile::OpenDatasource( datasource.c_str() ); + SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenDatasource returned " << (unsigned long)ds_id); + + OGRLayer* l_id = (OGRLayer *)tgShapefile::OpenLayer( ds_id, layer.c_str() ); + SG_LOG(SG_GENERAL, SG_DEBUG, "tgShapefile::OpenLayer returned " << (unsigned long)l_id); + + OGRPolygon* polygon = new OGRPolygon(); + + SG_LOG(SG_GENERAL, SG_DEBUG, "subject has " << subject.Contours() << " contours "); + + for ( unsigned int i = 0; i < subject.Contours(); i++ ) { + bool skip_ring=false; + tgContour contour = subject.GetContour( i ); + + if (contour.GetSize() < 3) { + SG_LOG(SG_GENERAL, SG_DEBUG, "Polygon with less than 3 points"); + skip_ring=true; + } + + // FIXME: Current we ignore the hole-flag and instead assume + // that the first ring is not a hole and the rest + // are holes + OGRLinearRing *ring=new OGRLinearRing(); + for (unsigned int pt = 0; pt < contour.GetSize(); pt++) { + OGRPoint *point=new OGRPoint(); + + point->setX( contour.GetNode(pt).getLongitudeDeg() ); + point->setY( contour.GetNode(pt).getLatitudeDeg() ); + point->setZ( 0.0 ); + ring->addPoint(point); + } + ring->closeRings(); + + if (!skip_ring) { + polygon->addRingDirectly(ring); + } + + OGRFeature* feature = NULL; + feature = new OGRFeature( l_id->GetLayerDefn() ); + feature->SetField("ID", description.c_str()); + feature->SetGeometry(polygon); + if( l_id->CreateFeature( feature ) != OGRERR_NONE ) + { + SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile"); + } + OGRFeature::DestroyFeature(feature); + } + + // close after each write + ds_id = tgShapefile::CloseDatasource( ds_id ); +} + +tgPolygon tgShapefile::ToPolygon( const void* subject ) +{ + const OGRPolygon* ogr_poly = (const OGRPolygon*)subject; + OGRLinearRing const *ring = ogr_poly->getExteriorRing(); + + tgContour contour; + tgPolygon result; + + for (int i = 0; i < ring->getNumPoints(); i++) { + contour.AddNode( SGGeod::fromDegM( ring->getX(i), ring->getY(i), ring->getZ(i)) ); + } + contour.SetHole( false ); + result.AddContour( contour ); + + // then add the inner rings + for ( int j = 0 ; j < ogr_poly->getNumInteriorRings(); j++ ) { + ring = ogr_poly->getInteriorRing( j ); + contour.Erase(); + + for (int i = 0; i < ring->getNumPoints(); i++) { + contour.AddNode( SGGeod::fromDegM( ring->getX(i), ring->getY(i), ring->getZ(i)) ); + } + contour.SetHole( true ); + result.AddContour( contour ); + } + result.SetTexMethod( TG_TEX_BY_GEODE ); + + return result; +} \ No newline at end of file diff --git a/src/Lib/terragear/tg_shapefile.hxx b/src/Lib/terragear/tg_shapefile.hxx new file mode 100644 index 00000000..f3917c50 --- /dev/null +++ b/src/Lib/terragear/tg_shapefile.hxx @@ -0,0 +1,20 @@ + +#include "tg_polygon.hxx" +#include "tg_contour.hxx" +#include "clipper.hpp" + +class tgShapefile +{ +public: + static void FromPolygon( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description ); + static tgPolygon ToPolygon( const void* subject ); + + static void FromClipper( const ClipperLib::Polygons& subject, const std::string& datasource, const std::string& layer, const std::string& description ); + +private: + static bool initialized; + + static void* OpenDatasource( const char* datasource_name ); + static void* OpenLayer( void* ds_id, const char* layer_name ); + static void* CloseDatasource( void* ds_id ); +}; \ No newline at end of file diff --git a/src/Lib/Polygon/surface.cxx b/src/Lib/terragear/tg_surface.cxx similarity index 99% rename from src/Lib/Polygon/surface.cxx rename to src/Lib/terragear/tg_surface.cxx index 15b6bbad..71ab44c4 100644 --- a/src/Lib/Polygon/surface.cxx +++ b/src/Lib/terragear/tg_surface.cxx @@ -32,7 +32,7 @@ #include #include "TNT/jama_qr.h" -#include "surface.hxx" +#include "tg_surface.hxx" // Final grid size for surface (in meters) const double coarse_grid = 300.0; diff --git a/src/Lib/Polygon/surface.hxx b/src/Lib/terragear/tg_surface.hxx similarity index 99% rename from src/Lib/Polygon/surface.hxx rename to src/Lib/terragear/tg_surface.hxx index f3a7cead..096d6105 100644 --- a/src/Lib/Polygon/surface.hxx +++ b/src/Lib/terragear/tg_surface.hxx @@ -28,7 +28,7 @@ #include #include "TNT/tnt_array2d.h" -#include "polygon.hxx" +#include "tg_polygon.hxx" /*** * A dirt simple matrix class for our convenience based on top of SGGeod diff --git a/src/Lib/Polygon/tg_unique_geod.hxx b/src/Lib/terragear/tg_unique_geod.hxx similarity index 100% rename from src/Lib/Polygon/tg_unique_geod.hxx rename to src/Lib/terragear/tg_unique_geod.hxx diff --git a/src/Lib/Polygon/tg_unique_tgnode.hxx b/src/Lib/terragear/tg_unique_tgnode.hxx similarity index 100% rename from src/Lib/Polygon/tg_unique_tgnode.hxx rename to src/Lib/terragear/tg_unique_tgnode.hxx diff --git a/src/Lib/Polygon/tg_unique_vec2f.hxx b/src/Lib/terragear/tg_unique_vec2f.hxx similarity index 100% rename from src/Lib/Polygon/tg_unique_vec2f.hxx rename to src/Lib/terragear/tg_unique_vec2f.hxx diff --git a/src/Lib/Polygon/tg_unique_vec3d.hxx b/src/Lib/terragear/tg_unique_vec3d.hxx similarity index 100% rename from src/Lib/Polygon/tg_unique_vec3d.hxx rename to src/Lib/terragear/tg_unique_vec3d.hxx diff --git a/src/Lib/Polygon/tg_unique_vec3f.hxx b/src/Lib/terragear/tg_unique_vec3f.hxx similarity index 100% rename from src/Lib/Polygon/tg_unique_vec3f.hxx rename to src/Lib/terragear/tg_unique_vec3f.hxx diff --git a/src/Prep/OGRDecode/CMakeLists.txt b/src/Prep/OGRDecode/CMakeLists.txt index d9f58b70..e2a41fc4 100644 --- a/src/Prep/OGRDecode/CMakeLists.txt +++ b/src/Prep/OGRDecode/CMakeLists.txt @@ -3,7 +3,7 @@ add_executable(ogr-decode ogr-decode.cxx) target_link_libraries(ogr-decode ${GDAL_LIBRARY} - Polygon + terragear ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} ) diff --git a/src/Prep/OGRDecode/ogr-decode.cxx b/src/Prep/OGRDecode/ogr-decode.cxx index 0ae9840d..acfe9645 100644 --- a/src/Prep/OGRDecode/ogr-decode.cxx +++ b/src/Prep/OGRDecode/ogr-decode.cxx @@ -32,13 +32,15 @@ #include #include #include - #include #include #include #include -#include + +#include +#include +#include /* stretch endpoints to reduce slivers in linear data ~.1 meters */ #define EP_STRETCH (0.1) @@ -195,7 +197,7 @@ void Decoder::processPolygon(OGRPolygon* poGeometry, const string& area_type ) // bool preserve3D = ((poGeometry->getGeometryType()&wkb25DBit)==wkb25DBit); // first add the outer ring - tgPolygon shape = tgPolygon::FromOGR( poGeometry ); + tgPolygon shape = tgShapefile::ToPolygon( poGeometry ); if ( max_segment_length > 0 ) { shape = tgPolygon::SplitLongEdges( shape, max_segment_length ); diff --git a/src/Utils/poly2ogr/CMakeLists.txt b/src/Utils/poly2ogr/CMakeLists.txt index 4101fbd1..96d5ca4a 100644 --- a/src/Utils/poly2ogr/CMakeLists.txt +++ b/src/Utils/poly2ogr/CMakeLists.txt @@ -7,7 +7,8 @@ endif(MSVC) target_link_libraries(poly2ogr ${GDAL_LIBRARY} - Polygon ${GETOPT_LIB} + terragear + ${GETOPT_LIB} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} ) diff --git a/src/Utils/poly2ogr/poly2ogr.cxx b/src/Utils/poly2ogr/poly2ogr.cxx index fec6fb98..73d1c4b8 100644 --- a/src/Utils/poly2ogr/poly2ogr.cxx +++ b/src/Utils/poly2ogr/poly2ogr.cxx @@ -45,7 +45,7 @@ #include #include -#include +#include using std::string;