From cd3a1d647eff439e95e85f2c684ac262c932cd1c Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Fri, 14 Dec 2012 10:14:54 -0500 Subject: [PATCH] - Use the tgPolygon SaveToGzFile / LoadFromGzFile functions to store the results of OGRDecode, instead of ascii, completely different format - move Utility functions to generate polys from line data into tgPolygon ( use this for both GenApts850 linear features, and OGRDecode ) - kill off old superpoly, and texparams - add tgChopper which queues chopped polygons per bucket. When saved, result is less files (each tile has 1 file per shapefile decoded) Both genapt and ogrdecode now use tgChopper instead of Geometry/Util.cxx functions - tgChooper acquires a global file name lock, so we can safely run genapts in paralell with ogrdecode. (this was dangerous before, as both could try to open the index file as write, the loser would segfault, later ) --- src/Airports/GenAirports850/CMakeLists.txt | 1 + src/Airports/GenAirports850/airport.cxx | 10 +- src/Airports/GenAirports850/linearfeature.cxx | 4 +- src/Airports/GenAirports850/linearfeature.hxx | 6 +- src/Airports/GenAirports850/main.cxx | 9 +- src/BuildTiles/Main/CMakeLists.txt | 2 + src/BuildTiles/Main/tgconstruct.hxx | 3 - src/BuildTiles/Main/tgconstruct_elevation.cxx | 2 + src/BuildTiles/Main/tgconstruct_poly.cxx | 228 +---- src/Lib/Array/array.hxx | 4 - src/Lib/Geometry/CMakeLists.txt | 2 - src/Lib/Geometry/tg_nodes.hxx | 4 - src/Lib/Geometry/util.cxx | 853 ---------------- src/Lib/Geometry/util.hxx | 129 --- src/Lib/Polygon/CMakeLists.txt | 4 - src/Lib/Polygon/chop-bin.cxx | 5 +- src/Lib/Polygon/chop.hxx | 6 - src/Lib/Polygon/polygon.cxx | 952 +++++++++++++----- src/Lib/Polygon/polygon.hxx | 93 +- src/Lib/Polygon/superpoly.cxx | 267 ----- src/Lib/Polygon/superpoly.hxx | 264 ----- src/Lib/Polygon/texparams.cxx | 46 - src/Lib/Polygon/texparams.hxx | 160 --- src/Prep/OGRDecode/ogr-decode.cxx | 505 ++++------ 24 files changed, 1013 insertions(+), 2546 deletions(-) delete mode 100644 src/Lib/Geometry/util.cxx delete mode 100644 src/Lib/Geometry/util.hxx delete mode 100644 src/Lib/Polygon/superpoly.cxx delete mode 100644 src/Lib/Polygon/superpoly.hxx delete mode 100644 src/Lib/Polygon/texparams.cxx delete mode 100644 src/Lib/Polygon/texparams.hxx diff --git a/src/Airports/GenAirports850/CMakeLists.txt b/src/Airports/GenAirports850/CMakeLists.txt index 01178f2e..6c2dca2a 100644 --- a/src/Airports/GenAirports850/CMakeLists.txt +++ b/src/Airports/GenAirports850/CMakeLists.txt @@ -1,3 +1,4 @@ +include_directories(${GDAL_INCLUDE_DIR}) include_directories(${PROJECT_SOURCE_DIR}/src/Lib) add_executable(genapts850 diff --git a/src/Airports/GenAirports850/airport.cxx b/src/Airports/GenAirports850/airport.cxx index f8640859..75b5e5fc 100644 --- a/src/Airports/GenAirports850/airport.cxx +++ b/src/Airports/GenAirports850/airport.cxx @@ -1360,7 +1360,11 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src ) } std::string holepath = root + "/AirportArea"; + tgChopper chopper( holepath ); + divided_base.SetPreserve3D( true ); - tgPolygon::Chop( divided_base, holepath, "Hole", false, true ); - tgPolygon::Chop( apt_clearing, holepath, "Airport", false, false ); -} + chopper.Add( divided_base, "Hole" ); + chopper.Add( apt_clearing, "Airport" ); + + chopper.Save(); +} \ No newline at end of file diff --git a/src/Airports/GenAirports850/linearfeature.cxx b/src/Airports/GenAirports850/linearfeature.cxx index 2e71f77c..69a20d57 100644 --- a/src/Airports/GenAirports850/linearfeature.cxx +++ b/src/Airports/GenAirports850/linearfeature.cxx @@ -1,7 +1,6 @@ #include #include -#include #include "global.hxx" #include "beznode.hxx" @@ -365,7 +364,7 @@ LinearFeature::~LinearFeature() } } - +#if 0 SGGeod LinearFeature::OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by ) { double courseCur, courseNext, courseAvg, theta; @@ -467,6 +466,7 @@ SGGeod midpoint( const SGGeod& p0, const SGGeod& p1 ) (p0.getLatitudeDeg() + p1.getLatitudeDeg()) / 2, (p0.getElevationM() + p1.getElevationM()) / 2 ); } +#endif int LinearFeature::Finish( bool closed, unsigned int idx ) { diff --git a/src/Airports/GenAirports850/linearfeature.hxx b/src/Airports/GenAirports850/linearfeature.hxx index 9923246c..61131597 100644 --- a/src/Airports/GenAirports850/linearfeature.hxx +++ b/src/Airports/GenAirports850/linearfeature.hxx @@ -95,9 +95,9 @@ public: int BuildBtg( tgpolygon_list& line_polys, tglightcontour_list& lights, tgAccumulator& accum, bool debug ); private: - SGGeod OffsetPointFirst( const SGGeod& cur, const SGGeod& next, double offset_by ); - SGGeod OffsetPointMiddle( const SGGeod& prev, const SGGeod& cur, const SGGeod& next, double offset_by ); - SGGeod OffsetPointLast( const SGGeod& prev, const SGGeod& cur, double offset_by ); +// SGGeod OffsetPointFirst( const SGGeod& cur, const SGGeod& next, double offset_by ); +// SGGeod OffsetPointMiddle( const SGGeod& prev, const SGGeod& cur, const SGGeod& next, double offset_by ); +// SGGeod OffsetPointLast( const SGGeod& prev, const SGGeod& cur, double offset_by ); double offset; double width; diff --git a/src/Airports/GenAirports850/main.cxx b/src/Airports/GenAirports850/main.cxx index 989ee451..0e035b8a 100644 --- a/src/Airports/GenAirports850/main.cxx +++ b/src/Airports/GenAirports850/main.cxx @@ -22,8 +22,6 @@ #include #include -#include -#include #include #include @@ -104,9 +102,6 @@ int nudge = 10; double slope_max = 0.02; double gSnap = 0.00000001; // approx 1 mm -//TODO : new polygon chop API -extern bool tgPolygon_index_init( const std::string& path ); - int main(int argc, char **argv) { SGGeod min = SGGeod::fromDeg( -180, -90 ); @@ -185,6 +180,7 @@ int main(int argc, char **argv) { max.setLatitudeDeg(atof( arg.substr(10).c_str() )); } +#if 0 // relly? - do we need this? else if ( arg.find("--chunk=") == 0 ) { tg::Rectangle rectangle = tg::parseChunk(arg.substr(8).c_str(), 10.0); @@ -197,6 +193,7 @@ int main(int argc, char **argv) min = rectangle.getMin(); max = rectangle.getMax(); } +#endif else if ( arg.find("--airport=") == 0 ) { airport_id = arg.substr(10).c_str(); @@ -290,7 +287,7 @@ int main(int argc, char **argv) // initialize persistant polygon counter std::string counter_file = airportareadir+"/poly_counter"; - tgPolygon_index_init( counter_file ); + tgPolygon::ChopIdxInit( counter_file ); tg::Rectangle boundingBox(min, max); boundingBox.sanify(); diff --git a/src/BuildTiles/Main/CMakeLists.txt b/src/BuildTiles/Main/CMakeLists.txt index f2ef5661..33eb62d2 100644 --- a/src/BuildTiles/Main/CMakeLists.txt +++ b/src/BuildTiles/Main/CMakeLists.txt @@ -1,3 +1,5 @@ +include_directories(${GDAL_INCLUDE_DIR}) + add_executable(tg-construct tgconstruct.hxx tgconstruct.cxx diff --git a/src/BuildTiles/Main/tgconstruct.hxx b/src/BuildTiles/Main/tgconstruct.hxx index db2d2665..018e2f24 100644 --- a/src/BuildTiles/Main/tgconstruct.hxx +++ b/src/BuildTiles/Main/tgconstruct.hxx @@ -122,9 +122,6 @@ private: // Load Data void LoadElevationArray( bool add_nodes ); int LoadLandclassPolys( void ); - // Load Data Helpers - bool load_poly(const std::string& path); - void add_poly(int area, tgPolygon& poly, std::string material); // Clip Data bool ClipLandclassPolys( void ); diff --git a/src/BuildTiles/Main/tgconstruct_elevation.cxx b/src/BuildTiles/Main/tgconstruct_elevation.cxx index cb1df472..35dc5081 100644 --- a/src/BuildTiles/Main/tgconstruct_elevation.cxx +++ b/src/BuildTiles/Main/tgconstruct_elevation.cxx @@ -81,6 +81,8 @@ void TGConstruct::CalcElevations( void ) } } + return; + nodes.get_geod_nodes(raw_nodes); // now flatten some stuff diff --git a/src/BuildTiles/Main/tgconstruct_poly.cxx b/src/BuildTiles/Main/tgconstruct_poly.cxx index e3edb37e..fb367ad3 100644 --- a/src/BuildTiles/Main/tgconstruct_poly.cxx +++ b/src/BuildTiles/Main/tgconstruct_poly.cxx @@ -37,187 +37,6 @@ using std::string; static unsigned int cur_poly_id = 0; -bool TGConstruct::load_poly(const string& path) { - bool poly3d = false; - bool with_tp = false; - string first_line; - string poly_name; - AreaType poly_type; - int contours, count, i, j, k; - int hole_flag; - int num_polys; - double startx, starty, startz, x, y, z, lastx, lasty, lastz; - - sg_gzifstream in( path ); - - if ( !in ) { - SG_LOG( SG_CLIPPER, SG_ALERT, "Cannot open file: " << path ); - exit(-1); - } - - tgPolygon poly; - tgTexParams tp; - - // (this could break things, why is it here) in >> skipcomment; - while ( !in.eof() ) { - in >> first_line; - if ( first_line == "#2D" ) { - poly3d = false; - with_tp = false; - - in >> poly_name; - num_polys = 1; - } else if ( first_line == "#2D_WITH_MASK" ) { - poly3d = false; - with_tp = false; - - in >> poly_name; - in >> num_polys; - } else if ( first_line == "#2D_WITH_TPS" ) { - poly3d = false; - with_tp = true; - - in >> poly_name; - in >> num_polys; - } else if ( first_line == "#3D" ) { - poly3d = true; - with_tp = false; - - in >> poly_name; - num_polys = 1; - } else { - // support old format (default to 2d) - poly3d = false; - with_tp = false; - - poly_name = first_line; - num_polys = 1; - } - poly_type = get_area_type( poly_name ); - - int area = (int)poly_type; - string material = get_area_name( area ); - - - // Generate a new Shape for the poly - tgPolygon poly; - SGGeod p; - - for (k=0; k> x; - in >> y; - tp.ref = SGGeod::fromDeg(x,y); - - in >> tp.width; - in >> tp.length; - in >> tp.heading; - in >> tp.minu; - in >> tp.maxu; - in >> tp.minv; - in >> tp.maxv; - poly.SetTexParams( tp ); - } - - in >> contours; - - poly.Erase(); - for ( i = 0; i < contours; ++i ) { - in >> count; - - if ( count < 3 ) { - SG_LOG( SG_CLIPPER, SG_ALERT, "Polygon with less than 3 data points." ); - exit(-1); - } - - in >> hole_flag; - - in >> startx; - in >> starty; - if ( poly3d ) { - in >> startz; - } else { - startz = -9999.0; - } - - p = SGGeod::fromDegM(startx+nudge, starty+nudge, startz ); - poly.AddNode( i, p ); - - if ( poly3d ) { - nodes.unique_add_fixed_elevation( p ); - } else { - nodes.unique_add( p ); - } - - for ( j = 1; j < count - 1; ++j ) { - in >> x; - in >> y; - if ( poly3d ) { - in >> z; - } else { - z = -9999.0; - } - - p = SGGeod::fromDegM( x+nudge, y+nudge, z ); - poly.AddNode( i, p ); - - if ( poly3d ) { - nodes.unique_add_fixed_elevation( p ); - } else { - nodes.unique_add( p ); - } - } - - in >> lastx; - in >> lasty; - if ( poly3d ) { - in >> lastz; - } else { - lastz = -9999.0; - } - - if ( (fabs(startx - lastx) < SG_EPSILON) && - (fabs(starty - lasty) < SG_EPSILON) && - (fabs(startz - lastz) < SG_EPSILON) ) { - // last point same as first, discard - } else { - p = SGGeod::fromDegM( lastx+nudge, lasty+nudge, lastz ); - poly.AddNode( i, p ); - - if ( poly3d ) { - nodes.unique_add_fixed_elevation( p ); - } else { - nodes.unique_add( p ); - } - } - } - - poly = tgPolygon::Snap( poly, gSnap ); - poly = tgPolygon::RemoveDups( poly ); - poly.SetMaterial( material ); - - if ( with_tp ) { - poly.SetTexMethod( TG_TEX_BY_TPS_CLIPU, -1, 0, 1, 0 ); - } else { - poly.SetTexMethod( TG_TEX_BY_GEODE, bucket.get_center_lat() ); - } - - in >> skipcomment; - - poly.SetId( cur_poly_id++ ); - polys_in.add_poly( area, poly ); - - if ( IsDebugShape( poly.GetId() ) ) { - tgPolygon::ToShapefile( poly, ds_name, "loaded", "" ); - } - } - } - - return true; -} - // load all 2d polygons from the specified load disk directories and // clip against each other to resolve any overlaps int TGConstruct::LoadLandclassPolys( void ) { @@ -225,7 +44,7 @@ int TGConstruct::LoadLandclassPolys( void ) { string base = bucket.gen_base_path(); string poly_path; - int count = 0; + int total_polys_read = 0; polys_in.clear(); @@ -254,12 +73,49 @@ int TGConstruct::LoadLandclassPolys( void ) { { // skipped! } else { - load_poly( p.str() ); + int area; + std::string material; + gzFile fp =gzopen( p.c_str(), "rb" ); + unsigned int count; + + sgReadUInt( fp, &count ); + SG_LOG( SG_GENERAL, SG_DEBUG, " Load " << count << " polys from " << p.realpath() ); + + for ( unsigned int i=0; i -#endif - #include #include #include diff --git a/src/Lib/Geometry/CMakeLists.txt b/src/Lib/Geometry/CMakeLists.txt index 47a4a50e..b16e164a 100644 --- a/src/Lib/Geometry/CMakeLists.txt +++ b/src/Lib/Geometry/CMakeLists.txt @@ -10,7 +10,6 @@ add_library(Geometry STATIC tg_nodes.cxx trinodes.cxx trisegs.cxx - util.cxx contour_tree.hxx line.hxx point3d.hxx @@ -20,5 +19,4 @@ add_library(Geometry STATIC tg_nodes.hxx trinodes.hxx trisegs.hxx - util.hxx ) diff --git a/src/Lib/Geometry/tg_nodes.hxx b/src/Lib/Geometry/tg_nodes.hxx index 840eabaf..39f9df3a 100644 --- a/src/Lib/Geometry/tg_nodes.hxx +++ b/src/Lib/Geometry/tg_nodes.hxx @@ -1,10 +1,6 @@ #ifndef _TG_NODES_HXX #define _TG_NODES_HXX -#ifdef HAVE_CONFIG_H -# include -#endif - #ifndef __cplusplus # error This library requires C++ #endif diff --git a/src/Lib/Geometry/util.cxx b/src/Lib/Geometry/util.cxx deleted file mode 100644 index c2050943..00000000 --- a/src/Lib/Geometry/util.cxx +++ /dev/null @@ -1,853 +0,0 @@ -// util.cxx - a collection of simple geometry utility functions. -// -// Started by David Megginson, July 2002 -// -// This file is in the Public Domain and comes with NO WARRANTY OF ANY KIND. - -#include "util.hxx" - -#include -#include -#include - -#include -#include -#include - -#include - -// use cgal for intersection, old implementation returns true, even when line SEGMENTS don't intersect -// CGAL intersection can do that with lines, but we want to use segments (end at the first and last point) - -#include -#include - -typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel; -typedef Kernel::Point_2 Point_2; -typedef CGAL::Segment_2 Segment_2; - -#define MP_STRETCH (0.000001) - -using std::string; -using std::vector; - -namespace tg { - - -/* - * Computer the intersection point of two lines. Reference: - * http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/ - */ -bool -getIntersection (const Point3D &p0, const Point3D &p1, - const Point3D &p2, const Point3D &p3, - Point3D &intersection) -{ - const double my_eps = 1E-12; - - double u_num = - ((p3.x()-p2.x())*(p0.y()-p2.y()))-((p3.y()-p2.y())*(p0.x()-p2.x())); - double u_den = - ((p3.y()-p2.y())*(p1.x()-p0.x()))-((p3.x()-p2.x())*(p1.y()-p0.y())); - - if ( fabs(u_den) < my_eps ) { - if ( fabs(u_num) < my_eps ) { - SG_LOG(SG_GENERAL, SG_ALERT, "Intersection: coincident lines"); - } else { - SG_LOG(SG_GENERAL, SG_ALERT, "Intersection: parallel lines"); - } - return false; - } else { - double u = u_num/u_den; - // cout << "u = " << u << " u_num = " << u_num << " u_den = " << u_den << endl; - intersection = Point3D((p0.x()+u*(p1.x()-p0.x())), - (p0.y()+u*(p1.y()-p0.y())), - 0); - return true; - } -} - -// use CGAL - -bool getIntersection_cgal(const Point3D &p0, const Point3D &p1, - const Point3D &p2, const Point3D &p3, - Point3D &intersection) -{ - Point_2 a1( p0.x(), p0.y() ); - Point_2 b1( p1.x(), p1.y() ); - Point_2 a2( p2.x(), p2.y() ); - Point_2 b2( p3.x(), p3.y() ); - - 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; - } - } -} - -/** - * Create a polygon out of a point. - * - * Note that simple geometry doesn't work here, because the scale is - * not even -- the points on the x-axis (longitude) become closer and - * closer as the y-axis (latitude) approaches the poles, meeting in - * a single point at y=90 and y=-90. As a result, this function - * uses the WGS80 functions, rather than simple Pythagorean stuff. - */ -void -makePolygon (const Point3D &p, double width, TGPolygon &polygon) -{ - double x = 0.0f, y = 0.0f, az = 0.0f; - double lon = p.x(); - double lat = p.y(); - - polygon.erase(); - - geo_direct_wgs_84(0, lat, lon, 90, width/2.0, &y, &x, &az); - double dlon = x - lon; - - geo_direct_wgs_84(0, lat, lon, 0, width/2.0, &y, &x, &az); - double dlat = y - lat; - - polygon.add_node(0, Point3D(lon - dlon, lat - dlat, 0)); - polygon.add_node(0, Point3D(lon + dlon, lat - dlat, 0)); - polygon.add_node(0, Point3D(lon + dlon, lat + dlat, 0)); - polygon.add_node(0, Point3D(lon - dlon, lat + dlat, 0)); - - polygon.set_hole_flag(0, 0); // mark as solid -} - - -void -makePolygon (const Line &line, double width, TGPolygon &polygon) -{ - vector segment_list; - - int nPoints = line.getPointCount(); - int i; - for (i = 0; i < nPoints - 1; i++) { - const Point3D p1 = line.getPoint(i); - const Point3D p2 = line.getPoint(i+1); - - double angle1 = 0.0f; - double angle2 = 0.0f; - double dist = 0.0f; - double x = 0.0f; - double y = 0.0f; - double az = 0.0f; - - geo_inverse_wgs_84(0, p1.y(), p1.x(), p2.y(), p2.x(), &angle1, &angle2, &dist); - polygon.erase(); - - // Wind each rectangle counterclockwise - - // Corner 1 - // cout << "point = " << i << endl; - geo_direct_wgs_84(0, p1.y(), p1.x(), CLAMP_ANGLE(angle1+90), width/2.0, &y, &x, &az); - // Sometimes creating the new point can wrap the date line, even if the - // original line does not. This can cause problems later, so lets just - // clip the result to date line and hope that the typical error is very - // small. FIXME: we should handle this more robustly in case larger - // structures cross the date line, but that implies a much broader swath - // of changes. - if ( x - p1.x() > 180.0 ) x = -180.0; - else if ( x - p1.x() < -180.0 ) x = 180.0; - polygon.add_node(0, Point3D(x, y, 0)); - // cout << "corner 1 = " << Point3D(x, y, 0) << endl; - - // Corner 2 - geo_direct_wgs_84(0, p2.y(), p2.x(), CLAMP_ANGLE(angle1+90), width/2.0, &y, &x, &az); - // Sometimes creating the new point can wrap the date line, even if the - // original line does not. This can cause problems later, so lets just - // clip the result to date line and hope that the typical error is very - // small. FIXME: we should handle this more robustly in case larger - // structures cross the date line, but that implies a much broader swath - // of changes. - if ( x - p2.x() > 180.0 ) x = -180.0; - else if ( x - p2.x() < -180.0 ) x = 180.0; - polygon.add_node(0, Point3D(x, y, 0)); - // cout << "corner 2 = " << Point3D(x, y, 0) << endl; - - // Corner 3 - geo_direct_wgs_84(0, p2.y(), p2.x(), CLAMP_ANGLE(angle1-90), width/2.0, &y, &x, &az); - // Sometimes creating the new point can wrap the date line, even if the - // original line does not. This can cause problems later, so lets just - // clip the result to date line and hope that the typical error is very - // small. FIXME: we should handle this more robustly in case larger - // structures cross the date line, but that implies a much broader swath - // of changes. - if ( x - p2.x() > 180.0 ) x = -180.0; - else if ( x - p2.x() < -180.0 ) x = 180.0; - polygon.add_node(0, Point3D(x, y, 0)); - // cout << "corner 3 = " << Point3D(x, y, 0) << endl; - - // Corner 4 - geo_direct_wgs_84(0, p1.y(), p1.x(), CLAMP_ANGLE(angle1-90), width/2.0, &y, &x, &az); - // Sometimes creating the new point can wrap the date line, even if the - // original line does not. This can cause problems later, so lets just - // clip the result to date line and hope that the typical error is very - // small. FIXME: we should handle this more robustly in case larger - // structures cross the date line, but that implies a much broader swath - // of changes. - if ( x - p1.x() > 180.0 ) x = -180.0; - else if ( x - p1.x() < -180.0 ) x = 180.0; - polygon.add_node(0, Point3D(x, y, 0)); - // cout << "corner 4 = " << Point3D(x, y, 0) << endl; - - // Save this rectangle - segment_list.push_back(polygon); - } - - // Build one big polygon out of all the rectangles by intersecting - // the lines running through the bottom and top sides - - polygon.erase(); - - // Connect the bottom part. - int nSegments = segment_list.size(); - Point3D intersection; - polygon.add_node(0, segment_list[0].get_pt(0, 0)); - // cout << "bnode = " << segment_list[0].get_pt(0, 0) << endl; - for (i = 0; i < nSegments - 1; i++) { - // cout << "point = " << i << endl; - if (getIntersection(segment_list[i].get_pt(0, 0), - segment_list[i].get_pt(0, 1), - segment_list[i+1].get_pt(0, 0), - segment_list[i+1].get_pt(0, 1), - intersection)) { - polygon.add_node(0, intersection); - // cout << "bnode (int) = " << intersection << endl; - } else { - polygon.add_node(0, segment_list[i].get_pt(0, 1)); - // cout << "bnode = " << segment_list[i].get_pt(0, 1) << endl; - } - } - polygon.add_node(0, segment_list[nSegments-1].get_pt(0, 1)); - // cout << "bnode = " << segment_list[nSegments-1].get_pt(0, 1) << endl; - - // Connect the top part - polygon.add_node(0, segment_list[nSegments-1].get_pt(0, 2)); - // cout << "tnode = " << segment_list[nSegments-1].get_pt(0, 2) << endl; - for (i = nSegments - 1; i > 0; i--) { - // cout << "point = " << i << endl; - if (getIntersection(segment_list[i].get_pt(0, 2), - segment_list[i].get_pt(0, 3), - segment_list[i-1].get_pt(0, 2), - segment_list[i-1].get_pt(0, 3), - intersection)) { - polygon.add_node(0, intersection); - // cout << "tnode (int) = " << intersection << endl; - } else { - polygon.add_node(0, segment_list[i].get_pt(0, 3)); - // cout << "tnode = " << segment_list[i].get_pt(0, 3) << endl; - } - } - polygon.add_node(0, segment_list[0].get_pt(0, 3)); - // cout << "tnode = " << segment_list[0].get_pt(0, 3) << endl; - - polygon.set_hole_flag(0, 0); // mark as solid -} - -inline double CalculateTheta( Point3D p0, Point3D p1, Point3D p2 ) -{ - Point3D u, v; - double udist, vdist, uv_dot, tmp; - - // u . v = ||u|| * ||v|| * cos(theta) - - u = p1 - p0; - udist = sqrt( u.x() * u.x() + u.y() * u.y() ); - // printf("udist = %.6f\n", udist); - - v = p1 - p2; - vdist = sqrt( v.x() * v.x() + v.y() * v.y() ); - // printf("vdist = %.6f\n", vdist); - - uv_dot = u.x() * v.x() + u.y() * v.y(); - // printf("uv_dot = %.6f\n", uv_dot); - - tmp = uv_dot / (udist * vdist); - // printf("tmp = %.6f\n", tmp); - - return acos(tmp); -} - -Point3D OffsetPointMiddle( Point3D prev, Point3D cur, Point3D next, double offset_by, int& turn_dir ) -{ - double offset_dir; - double next_dir; - double az2; - double dist; - double theta; - double pt_x = 0, pt_y = 0; - - SG_LOG(SG_GENERAL, SG_DEBUG, "Find average angle for contour: prev (" << prev << "), " - "cur (" << cur << "), " - "next (" << next << ")" ); - - - // 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 - // if the cross product is 0, then we need to use the direction passed in - SGVec3d dir1 = prev.toSGVec3d() - cur.toSGVec3d(); - dir1 = normalize(dir1); - - SGVec3d dir2 = next.toSGVec3d() - cur.toSGVec3d(); - dir2 = normalize(dir2); - - // Now find the average - SGVec3d avg = dir1 + dir2; - avg = normalize(avg); - - // check the turn direction - SGVec3d cp = cross( dir1, dir2 ); - SG_LOG(SG_GENERAL, SG_DEBUG, "\tcross product of dir1: " << dir1 << " and dir2: " << dir2 << " is " << cp ); - - // calculate the angle between cur->prev and cur->next - theta = SGMiscd::rad2deg(CalculateTheta(prev, cur, next)); - - if ( abs(theta - 180.0) < 0.1 ) - { - SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta close to 180) theta is " << theta ); - - // find the direction to the next point - geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &next_dir, &az2, &dist); - - offset_dir = next_dir - 90.0; - while (offset_dir < 0.0) - { - offset_dir += 360.0; - } - - // straight line blows up math - dist should be exactly as given - dist = offset_by; - } - else if ( abs(theta) < 0.1 ) - { - SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta close to 0) : theta is " << theta ); - - // find the direction to the next point - geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &next_dir, &az2, &dist); - - offset_dir = next_dir - 90; - while (offset_dir < 0.0) - { - offset_dir += 360.0; - } - - // straight line blows up math - dist should be exactly as given - dist = offset_by; - } - else - { - SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta NOT close to 180) : theta is " << theta ); - - // find the offset angle - geo_inverse_wgs_84( avg.y(), avg.x(), 0.0f, 0.0f, &offset_dir, &az2, &dist); - - // if we turned right, reverse the heading - if (cp.z() < 0.0f) - { - turn_dir = 0; - offset_dir += 180.0; - } - else - { - turn_dir = 1; - } - - while (offset_dir >= 360.0) - { - offset_dir -= 360.0; - } - - // find the direction to the next point - geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &next_dir, &az2, &dist); - - // calculate correct distance for the offset point - dist = (offset_by)/sin(SGMiscd::deg2rad(next_dir-offset_dir)); - } - - SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << offset_dir << " distance is " << dist ); - - // calculate the point from cur - geo_direct_wgs_84( cur.y(), cur.x(), offset_dir, dist, &pt_y, &pt_x, &az2 ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "\tpoint is (" << pt_x << "," << pt_y << ")" ); - - return Point3D(pt_x, pt_y, 0.0f); -} - -Point3D OffsetPointFirst( Point3D cur, Point3D next, double offset_by ) -{ - double offset_dir; - double az2; - double dist; - double pt_x = 0, pt_y = 0; - - SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at Start : cur (" << cur << "), " - "next (" << next << ")" ); - - // find the offset angle - geo_inverse_wgs_84( cur.y(), cur.x(), next.y(), next.x(), &offset_dir, &az2, &dist); - offset_dir -= 90; - if (offset_dir < 0) - { - offset_dir += 360; - } - - SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << offset_dir << " distance is " << offset_by ); - - // calculate the point from cur - geo_direct_wgs_84( cur.y(), cur.x(), offset_dir, offset_by, &pt_y, &pt_x, &az2 ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "\tpoint is (" << pt_x << "," << pt_y << ")" ); - - return Point3D(pt_x, pt_y, 0.0f); -} - -Point3D OffsetPointLast( Point3D prev, Point3D cur, double offset_by ) -{ - double offset_dir; - double az2; - double dist; - double pt_x = 0, pt_y = 0; - - SG_LOG(SG_GENERAL, SG_DEBUG, "Find OffsetPoint at End : prev (" << prev << "), " - "cur (" << cur << ")" ); - - // find the offset angle - geo_inverse_wgs_84( prev.y(), prev.x(), cur.y(), cur.x(), &offset_dir, &az2, &dist); - offset_dir -= 90; - if (offset_dir < 0) - { - offset_dir += 360; - } - - SG_LOG(SG_GENERAL, SG_DEBUG, "\theading is " << offset_dir << " distance is " << offset_by ); - - // calculate the point from cur - geo_direct_wgs_84( cur.y(), cur.x(), offset_dir, offset_by, &pt_y, &pt_x, &az2 ); - - SG_LOG(SG_GENERAL, SG_DEBUG, "\tpoint is (" << pt_x << "," << pt_y << ")" ); - - return Point3D(pt_x, pt_y, 0.0f); -} - -Point3D midpoint( Point3D p0, Point3D p1 ) -{ - return Point3D( (p0.x() + p1.x()) / 2, (p0.y() + p1.y()) / 2, (p0.z() + p1.z()) / 2 ); -} - -void -makePolygons (const Line &line, double width, poly_list& polys) -{ - int nPoints = line.getPointCount(); - int i; - int turn_dir; - - Point3D cur_inner; - Point3D cur_outer; - Point3D prev_inner = Point3D(0.0f, 0.0f, 0.0f); - Point3D prev_outer = Point3D(0.0f, 0.0f, 0.0f); - - double heading = 0.0f; - double az2 = 0.0f; - double dist = 0.0f; - double pt_x = 0.0f; - double pt_y = 0.0f; - - TGPolygon poly; - - // generate poly and texparam lists for each line segment - for (i=0; i p.x()) - min_x = p.x(); - if (max_x < p.x()) - max_x = p.x(); - if (min_y > p.y()) - min_y = p.y(); - if (max_y < p.y()) - max_y = p.y(); - } - } - } - return Rectangle(SGGeod::fromDeg(min_x, min_y), SGGeod::fromDeg(max_x, max_y)); -} - - -Rectangle -parseChunk (const string &s, double delta) -{ - Rectangle bounds; - int x_factor; - int y_factor; - - if (s.size() != 7) - throw sg_exception(string("Bad length for chunk specifier: ") + s); - - if (s[0] == 'w') - x_factor = -1; - else if (s[0] == 'e') - x_factor = 1; - else - throw sg_exception(string("Chunk specifier must begin with 'e' or 'w': " - + s)); - - if (s[4] == 's') - y_factor = -1; - else if (s[4] == 'n') - y_factor = 1; - else - throw sg_exception("Second part of chunk specifier must begin with 's' or 'n': " + s); - - - double x = atoi(s.substr(1,3).c_str()) * x_factor; - double y = atoi(s.substr(5).c_str()) * y_factor; - bounds.setMin(SGGeod::fromDeg(x, y)); - bounds.setMax(SGGeod::fromDeg(x + delta, y + delta)); - - return bounds; -} - -Rectangle -parseTile (const string &s) -{ - return( parseChunk(s, 1.0) ); -} - -}; - -// end of util.cxx diff --git a/src/Lib/Geometry/util.hxx b/src/Lib/Geometry/util.hxx deleted file mode 100644 index 22ba6df0..00000000 --- a/src/Lib/Geometry/util.hxx +++ /dev/null @@ -1,129 +0,0 @@ -// util.hxx - a collection of simple WGS84 utility functions. -// -// Started by David Megginson, July 2002 -// -// This file is in the Public Domain and comes with NO WARRANTY OF ANY KIND. - -#ifndef __UTIL_HXX -#define __UTIL_HXX 1 - -#ifndef __cplusplus -# error This library requires C++ -#endif - -#include -#include - -#include - -#include -#include -#include - -#include "line.hxx" - -namespace tg { - - -/** - * Inline function to clamp an angle between 0 and 360 degrees. - * - * @param a The angle to clamp, in degrees. - * @return The clamped angle. - */ -inline double -CLAMP_ANGLE (double a) -{ - while (a < 0.0) - a += 360.0; - while (a >= 360.0) - a -= 360.0; - return a; -} - - -/** - * Calculate the intersection of two line segments. - * - * @param p0 First point on the first segment. - * @param p1 A second point on the second segment. - * @param p2 First point on the second segment. - * @param p3 A second point on the second segment. - * @param intersection A variable to hold the calculated intersection. - * @return true if there was an intersection, false if the segments. - * are parallel or coincident. - */ -bool getIntersection (const Point3D &p0, const Point3D &p1, - const Point3D &p2, const Point3D &p3, - Point3D &intersection); - - -/** - * Create a polygon out of a point. - * - * Note that simple geometry doesn't work here, because the scale is - * not even -- the points on the x-axis (longitude) become closer and - * closer as the y-axis (latitude) approaches the poles, meeting in - * a single point at y=90 and y=-90. As a result, this function - * uses the WGS80 functions, rather than simple Pythagorean stuff. - * - * @param p The point at the centre of the new polygon. - * @param width The width in standard units (meters for FlightGear). - * @param polygon The object that will hold the new polygon. - */ -void makePolygon (const Point3D &p, double width, TGPolygon &polygon); - - -/** - * Create a polygon out of a line. - * - * Note that simple geometry doesn't work here, because the scale is - * not even -- the points on the x-axis (longitude) become closer and - * closer as the y-axis (latitude) approaches the poles, meeting in - * a single point at y=90 and y=-90. As a result, this function - * uses the WGS80 functions, rather than simple Pythagorean stuff. - * - * Update CLO 7 Mar 2003: Function rewritten to carefully preserve - * specifed polygon width through any kind of point to point angle - * change. - * - * @param line The multi-segment line inside the new polygon. - * @param width The width in standard units (meters for FlightGear). - * @param polygon The object that will hold the new polygon. - */ -void makePolygon (const Line &line, double width, TGPolygon &polygon); - -void makePolygons (const Line &line, double width, poly_list &segments); -void makePolygonsTP (const Line &line, double width, poly_list &segments, texparams_list &tps); - -/** - * Make a 2D bounding rectangle for a polygon. - * - * @param polygon The polygon to make the rectangle for. - * @return The bounding rectangle. - */ -Rectangle makeBounds (const TGPolygon &polygon); - - -/** - * Parse a chunk string (like "w080n40") into a rectangle. - * Defaults to 10 x 10 degrees. - * - * @param s The string. - * @return A rectangle containing the bounds. - */ -Rectangle parseChunk (const std::string &s, double delta); - - -/** - * Parse a tile string (like "w080n41") into a 1x1 degree rectangle. - * - * @param s The string. - * @return A rectangle containing the bounds. - */ -Rectangle parseTile (const std::string &s); - - -}; // namespace tg - -#endif // __UTIL_HXX diff --git a/src/Lib/Polygon/CMakeLists.txt b/src/Lib/Polygon/CMakeLists.txt index 221b0e01..6e67da07 100644 --- a/src/Lib/Polygon/CMakeLists.txt +++ b/src/Lib/Polygon/CMakeLists.txt @@ -7,15 +7,11 @@ add_library(Polygon STATIC index.cxx polygon.cxx simple_clip.cxx - superpoly.cxx chop.hxx index.hxx names.hxx polygon.hxx simple_clip.hxx - superpoly.hxx - texparams.cxx - texparams.hxx clipper.cpp clipper.hpp ) diff --git a/src/Lib/Polygon/chop-bin.cxx b/src/Lib/Polygon/chop-bin.cxx index ffe13dda..e762d5fc 100644 --- a/src/Lib/Polygon/chop-bin.cxx +++ b/src/Lib/Polygon/chop-bin.cxx @@ -261,6 +261,7 @@ static void clip_and_write_polys_with_mask( string root, long int p_index, fclose( rfp ); } +#if 0 static void clip_and_write_polys_with_tps( string root, long int p_index, const string &poly_type, SGBucket b, @@ -493,7 +494,7 @@ static void clip_and_write_poly_tp( string root, long int p_index, fclose( rfp ); } } - +#endif // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) @@ -799,6 +800,7 @@ void tgChopNormalPolygonsWithMask(const std::string& path, const std::string& po } } +#if 0 // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) void tgChopNormalPolygonTP( const string& path, const string& poly_type, @@ -1106,6 +1108,7 @@ void tgChopNormalPolygonsWithTP(const std::string& path, const std::string& poly tgChopNormalPolygonsWithTP( path, poly_type, top_clip_list, top_tp_list, preserve3d ); } } +#endif // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) This has a front end to a crude clipper diff --git a/src/Lib/Polygon/chop.hxx b/src/Lib/Polygon/chop.hxx index 7f97b9c0..9ffa0612 100644 --- a/src/Lib/Polygon/chop.hxx +++ b/src/Lib/Polygon/chop.hxx @@ -29,7 +29,6 @@ #include #include "polygon.hxx" -#include "texparams.hxx" // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) @@ -39,11 +38,6 @@ void tgChopNormalPolygon( const std::string& path, const std::string& poly_type, void tgChopNormalPolygonsWithMask(const std::string& path, const std::string& poly_type, const poly_list& segments, bool preserve3d ); -// process polygon shape (chop up along tile boundaries and write each -// polygon piece to a file) -void tgChopNormalPolygonsWithTP( const std::string& path, const std::string& poly_type, - const poly_list& segments, const texparams_list& tps, bool preserve3d ); - // process polygon shape (chop up along tile boundaries and write each // polygon piece to a file) This has a front end to a crude clipper // that doesn't handle holes so beware. This routine is appropriate diff --git a/src/Lib/Polygon/polygon.cxx b/src/Lib/Polygon/polygon.cxx index 1fada60c..103fe914 100644 --- a/src/Lib/Polygon/polygon.cxx +++ b/src/Lib/Polygon/polygon.cxx @@ -970,7 +970,7 @@ std::ostream& operator << (std::ostream &output, const TGPolygon &poly) return output; } -void TGPolygon::SaveToGzFile(gzFile& fp) +void TGPolygon::SaveToGzFile(gzFile& fp) const { int nContours = poly.size(); @@ -1055,6 +1055,51 @@ void TGPolygon::LoadFromGzFile(gzFile& fp) // NEW IMPLEMENTATIONS ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include +#include +#include +#include +#include +#include + + + + +// 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; @@ -1962,7 +2007,7 @@ bool tgContour::FindColinearLine( const tgContour& subject, const SGGeod& node, return false; } -void tgContour::SaveToGzFile( gzFile& fp ) +void tgContour::SaveToGzFile( gzFile& fp ) const { // Save the nodelist sgWriteUInt( fp, node_list.size() ); @@ -2240,6 +2285,314 @@ tgPolygon tgPolygon::Expand( const tgPolygon& subject, double offset ) return result; } +#if 0 +inline double CalculateTheta( const SGGeod& p0, const SGGeod& p1, const SGGeod& p2 ) +{ + SGVec2d v0, v1, v2; + SGVec2d u, v; + double udist, vdist, uv_dot; + + v0 = SGVec2d( p0.getLongitudeDeg(), p0.getLatitudeDeg() ); + v1 = SGVec2d( p1.getLongitudeDeg(), p1.getLatitudeDeg() ); + v2 = SGVec2d( p2.getLongitudeDeg(), p2.getLatitudeDeg() ); + + u = v1 - v0; + udist = norm(u); + + v = v1 - v2; + vdist = norm(v); + + uv_dot = dot(u, v); + + return acos( uv_dot / (udist * vdist) ); +} +#endif + +inline double CalculateTheta( const SGVec3d& dirCur, const SGVec3d& dirNext, const SGVec3d& cp ) +{ + double dp = dot( dirCur, dirNext ); + + return acos( dp ); +} + +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, cp ) ); + + 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 getIntersection_cgal(const SGGeod &p0, const SGGeod &p1, const SGGeod& p2, const SGGeod& p3, SGGeod& intersection) +{ + 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; + } + } +} + +#include + +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; +// double heading = 0.0f; +// double az2 = 0.0f; +// double dist = 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; @@ -2698,7 +3051,6 @@ void tgPolygon::Texture( void ) } } -#include void tgPolygon::ToShapefile( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& description ) { void* ds_id = tgShapefileOpenDatasource( datasource.c_str() ); @@ -2753,15 +3105,35 @@ void tgPolygon::ToShapefile( const tgPolygon& subject, const std::string& dataso ds_id = tgShapefileCloseDatasource( ds_id ); } - tgcontour_list contours; - tgtriangle_list triangles; +tgPolygon tgPolygon::FromOGR( const OGRPolygon* subject ) +{ + OGRLinearRing const *ring = subject->getExteriorRing(); + tgContour contour; + tgPolygon result; - std::string material; - std::string flag; // let's get rid of this.... - tgTexParams tp; + 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(); -void tgPolygon::SaveToGzFile( gzFile& fp ) + 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() ); @@ -2781,6 +3153,7 @@ void tgPolygon::SaveToGzFile( gzFile& fp ) // and the rest sgWriteString( fp, material.c_str() ); sgWriteString( fp, flag.c_str() ); + sgWriteInt( fp, (int)preserve3d ); } void tgPolygon::LoadFromGzFile( gzFile& fp ) @@ -2822,6 +3195,8 @@ void tgPolygon::LoadFromGzFile( gzFile& fp ) flag = strbuff; delete strbuff; } + + sgReadInt( fp, (int *)&preserve3d ); } ////////////////////////////// CHOP //////////////////////////////// @@ -2832,7 +3207,7 @@ void tgPolygon::LoadFromGzFile( gzFile& fp ) static long int poly_index = 0; static std::string poly_path; -bool tgPolygon_index_init( const std::string& path ) +bool tgPolygon::ChopIdxInit( const std::string& path ) { poly_path = path; FILE *fp = fopen( poly_path.c_str(), "r" ); @@ -2868,264 +3243,7 @@ static long int tgPolygon_index_next() return poly_index; } -static void ClipToFile( const tgPolygon& subject, std::string root, - long int p_index, - const std::string& type, - SGBucket& b, - bool withTexparams, - bool preserve3d ) -{ - Point3D p; - - SGGeod min, max; - SGGeod c = b.get_center(); - double span = b.get_width(); - - tgPolygon base, result; - char tile_name[256]; - char poly_index[256]; - - // 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 ( preserve3d ) { - result.InheritElevations( subject ); - } - - tgTexParams tp = subject.GetTexParams(); - - if ( result.Contours() > 0 ) { - long int t_index = b.gen_index(); - std::string path = root + "/" + b.gen_base_path(); - - SGPath sgp( path ); - sgp.append( "dummy" ); - sgp.create_dir( 0755 ); - - sprintf( tile_name, "%ld", t_index ); - std::string polyfile = path + "/" + tile_name; - - sprintf( poly_index, "%ld", p_index ); - polyfile += "."; - polyfile += poly_index; - - FILE *rfp = fopen( polyfile.c_str(), "w" ); - if ( preserve3d && withTexparams ) { - fprintf( rfp, "#3D_TP\n" ); - } else if ( withTexparams ) { - fprintf( rfp, "#2D_TP\n" ); - } else if ( preserve3d ) { - fprintf( rfp, "#3D\n" ); - } else { - fprintf( rfp, "#2D\n" ); - } - - fprintf( rfp, "%s\n", type.c_str() ); - - if ( withTexparams ) { - fprintf( rfp, "%.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", - tp.ref.getLongitudeDeg(), tp.ref.getLatitudeDeg(), - tp.width, tp.length, - tp.heading, - tp.minu, tp.maxu, tp.minu, tp.maxu); - } - - fprintf( rfp, "%d\n", result.Contours() ); - for ( unsigned int i = 0; i < result.Contours(); ++i ) { - fprintf( rfp, "%d\n", result.ContourSize(i) ); - fprintf( rfp, "%d\n", result.GetContour(i).GetHole() ); - for ( unsigned int j = 0; j < result.ContourSize(i); ++j ) { - p = Point3D::fromSGGeod( result.GetNode( i, j ) ); - if ( preserve3d ) - fprintf( rfp, "%.15f %.15f %.15f\n", p.x(), p.y(), p.z() ); - else - fprintf( rfp, "%.15f %.15f\n", p.x(), p.y() ); - } - } - fclose( rfp ); - } -} - -void tgPolygon::Chop( const tgPolygon& subject, const std::string& path, const std::string& type, bool withTexparams, bool preserve3d ) -{ - tg::Rectangle bb; - SGGeod p; - long int index; - - // bail out immediately if polygon is empty - if ( subject.Contours() == 0 ) - return; - - bb = subject.GetBoundingBox(); - - // get next polygon index - index = tgPolygon_index_next(); - - 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 - ClipToFile( subject, path, index, type, b_min, withTexparams, preserve3d ); - 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); - ClipToFile( subject, path, index, type, b_cur, withTexparams, preserve3d ); - } - } - 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().getLongitudeDeg(), 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( bottom, subject ); - - // the texparam should be constant over each clipped poly. - // when they are reassembled, we want the texture map to - // be seamless - tgPolygon::Chop( bottom_clip, path, type, withTexparams, preserve3d ); - } - - { - // - // 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( top, subject ); - - tgPolygon::Chop( top_clip, path, type, withTexparams, preserve3d ); - } -} - /************************ TESSELATION ***********************************/ -#include -#include -#include -#include -#include -#include - -/* 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; void tg_mark_domains(CDT& ct, CDT::Face_handle start, int index, std::list& border ) { @@ -3488,7 +3606,7 @@ void tgAccumulator::Add( const tgPolygon& subject ) accum.push_back( clipper_subject ); } -void tgTriangle::SaveToGzFile( gzFile& fp ) +void tgTriangle::SaveToGzFile( gzFile& fp ) const { // Save the three nodes, and their attributes for (unsigned int i = 0; i < 3; i++) { @@ -3518,7 +3636,7 @@ void tgTriangle::LoadFromGzFile( gzFile& fp ) // sgReadDouble( fp, &face_area ); } -void tgTexParams::SaveToGzFile( gzFile& fp ) +void tgTexParams::SaveToGzFile( gzFile& fp ) const { // Save the parameters sgWriteInt( fp, (int)method ); @@ -3580,4 +3698,298 @@ void tgTexParams::LoadFromGzFile( gzFile& fp ) sgReadDouble( fp, &max_clipv ); } } -} \ No newline at end of file +} + +// CHOPPER +void tgChopper::Clip( const tgPolygon& subject, + const std::string& type, + SGBucket& b ) +{ + Point3D p; + + SGGeod min, max; + SGGeod c = b.get_center(); + double span = b.get_width(); + + tgPolygon base, result; +// char tile_name[256]; + + // 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); + +// fprintf( rfp, "%s\n", type.c_str() ); +// +// if ( withTexparams ) { +// fprintf( rfp, "%.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f %.15f\n", +// tp.ref.getLongitudeDeg(), tp.ref.getLatitudeDeg(), +// tp.width, tp.length, +// tp.heading, +// tp.minu, tp.maxu, tp.minv, tp.maxv); +// } + +// fprintf( rfp, "%d\n", result.Contours() ); +// for ( unsigned int i = 0; i < result.Contours(); ++i ) { +// fprintf( rfp, "%d\n", result.ContourSize(i) ); +// fprintf( rfp, "%d\n", result.GetContour(i).GetHole() ); +// for ( unsigned int j = 0; j < result.ContourSize(i); ++j ) { +// p = Point3D::fromSGGeod( result.GetNode( i, j ) ); +// if ( preserve3d ) +// fprintf( rfp, "%.15f %.15f %.15f\n", p.x(), p.y(), p.z() ); +// else +// fprintf( rfp, "%.15f %.15f\n", p.x(), p.y() ); +// } +// } +// fclose( rfp ); + + lock.lock(); + bp_map[b.gen_index()].push_back( result ); + lock.unlock(); + } +} + +void tgChopper::Add( const tgPolygon& subject, const std::string& type ) +{ + tg::Rectangle 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 ); + } +} + +#include +#include + +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 -#include #include @@ -206,7 +205,7 @@ public: // output void write_contour( const int contour, const std::string& file ) const; - void SaveToGzFile( gzFile& fp ); + void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); // Friends for serialization @@ -311,14 +310,33 @@ std::ostream &operator<<(std::ostream &output, const TGPolygon &poly); // For the first attempt, I will just have TGPolygon with a list of contours. // the extra data are stored in paralell vectors. // This should also make TGSuperPoly obsolete +#include +#include + +#include +#include + #include #include #include "tg_unique_geod.hxx" +// utilities - belong is simgear? +double CalculateTheta( const SGVec3d& dirCur, const SGVec3d& dirNext, const SGVec3d& cp ); +SGGeod midpoint( const SGGeod& p0, const SGGeod& p1 ); +SGGeod OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by, int& turn_dir ); +SGGeod OffsetPointMiddle( const SGGeod& gPrev, const SGGeod& gCur, const SGGeod& gNext, double offset_by ); +SGGeod OffsetPointFirst( const SGGeod& cur, const SGGeod& next, double offset_by ); +SGGeod OffsetPointLast( const SGGeod& prev, const SGGeod& cur, double offset_by ); // forward declaration +// Forward Declaration: class tgPolygon; +class tgChoppedPolygons; + +typedef std::vector tgpolygon_list; +typedef tgpolygon_list::iterator tgpolygon_list_iterator; +typedef tgpolygon_list::const_iterator const_tgpolygon_list_iterator; class tgContour { @@ -392,10 +410,11 @@ public: 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 ); + void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); // Friend for output @@ -467,7 +486,7 @@ public: return face_area; } - void SaveToGzFile( gzFile& fp ); + void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); // Friend for output @@ -517,23 +536,20 @@ public: double center_lat; - void SaveToGzFile( gzFile& fp ); + void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); // Friend for output friend std::ostream& operator<< ( std::ostream&, const tgTexParams& ); }; -// Forward Declaration: -class tgPolygon; -typedef std::vector tgpolygon_list; -typedef tgpolygon_list::iterator tgpolygon_list_iterator; -typedef tgpolygon_list::const_iterator const_tgpolygon_list_iterator; -typedef void (*tgpolygon_texfunc)(void); - class tgPolygon { public: + tgPolygon() { + preserve3d = false; + } + void Erase( void ) { contours.clear(); triangles.clear(); @@ -623,6 +639,20 @@ public: material = m; } + std::string const& GetFlag( void ) const { + return flag; + } + void SetFlag( const std::string& f ) { + flag = f; + } + + bool GetPreserve3D( void ) const { + return preserve3d; + } + void SetPreserve3D( bool p ) { + preserve3d = p; + } + unsigned int GetId( void ) const { return id; } @@ -649,6 +679,7 @@ public: tp.maxu = maxu; tp.maxv = maxv; } + void SetTexMethod( tgTexMethod m ) { tp.method = m; } @@ -663,6 +694,9 @@ public: tp.method = m; tp.center_lat = cl; } + tgTexMethod GetTexMethod( void ) const { + return tp.method; + } void Tesselate( void ); void Tesselate( const std::vector& extra ); @@ -674,6 +708,8 @@ public: // 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 tgPolygon Union( const tgPolygon& subject, tgPolygon& clip ); @@ -685,6 +721,8 @@ public: static ClipperLib::Polygons ToClipper( const tgPolygon& subject ); static tgPolygon FromClipper( const ClipperLib::Polygons& subject ); + static tgPolygon FromOGR( const OGRPolygon* subject ); + // other operations static tgPolygon Snap( const tgPolygon& subject, double snap ); static tgPolygon StripHoles( const tgPolygon& subject ); @@ -698,8 +736,8 @@ public: static tgPolygon RemoveSpikes( const tgPolygon& subject ); static tgPolygon Expand( const tgPolygon& subject, double offset ); - - static void Chop( const tgPolygon& subject, const std::string& path, const std::string& type, bool withTexparams, bool preserve3d ); + static tgPolygon Expand( const SGGeod& subject, double offset ); + static void ToShapefile( const tgPolygon& subject, const std::string& datasource, const std::string& layer, const std::string& feature ); static void Tesselate( const tgPolygon& subject ); @@ -711,7 +749,7 @@ public: static void RemoveSlivers( tgPolygon& subject, tgcontour_list& slivers ); static tgcontour_list MergeSlivers( tgpolygon_list& subjects, tgcontour_list& slivers ); - void SaveToGzFile( gzFile& fp ); + void SaveToGzFile( gzFile& fp ) const; void LoadFromGzFile( gzFile& fp ); // Friend for output @@ -723,10 +761,35 @@ private: std::string material; std::string flag; // let's get rid of this.... + bool preserve3d; unsigned int id; // unique polygon id for debug 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: diff --git a/src/Lib/Polygon/superpoly.cxx b/src/Lib/Polygon/superpoly.cxx deleted file mode 100644 index a380a35d..00000000 --- a/src/Lib/Polygon/superpoly.cxx +++ /dev/null @@ -1,267 +0,0 @@ -// superpoly.cxx -- Manage all aspects of a rendered polygon -// -// Written by Curtis Olson, started June 2000. -// -// Copyright (C) 2000 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. -// -// $Id: superpoly.cxx,v 1.6 2004-11-19 22:25:50 curt Exp $ - - -#include "superpoly.hxx" -#include - - -// Constructor -TGSuperPoly::TGSuperPoly() : - flag( "" ) -{ -} - - -// Destructor -TGSuperPoly::~TGSuperPoly() -{ -} - -// erase the "super" polygon -void TGSuperPoly::erase() -{ - material = ""; - poly.erase(); - normals.erase(); - texcoords.erase(); - tris.erase(); - face_normals.clear(); -} - -// Friends for serialization -std::ostream& operator<< ( std::ostream& output, const TGSuperPoly& sp ) -{ - int nFaceNormals; - int nFaceAreas; - - // Save the data - output << sp.material << "\n"; - output << sp.poly; - output << sp.normals; - output << sp.texcoords; - - output << sp.tris; - output << sp.tri_idxs; - - nFaceNormals = sp.face_normals.size(); - output << nFaceNormals << "\n"; - for ( int i = 0; i < nFaceNormals; i++ ) { - output << sp.face_normals[i]; - } - - nFaceAreas = sp.face_areas.size(); - output << nFaceAreas; - for ( int i = 0; i < nFaceAreas; i++ ) { - output << sp.face_areas[i] << " "; - } - output << "\n"; - - if ( sp.flag.empty() ) { - output << "none\n"; - } else { - output << sp.flag << "\n"; - } - - return output; -} - -void TGSuperPoly::SaveToGzFile(gzFile& fp) -{ - int nFaceNormals; - int nFaceAreas; - - // Save the data - sgWriteString( fp, material.c_str() ); - - poly.SaveToGzFile( fp ); - normals.SaveToGzFile( fp ); - texcoords.SaveToGzFile( fp ); - - tris.SaveToGzFile( fp ); - tri_idxs.SaveToGzFile( fp ); - - nFaceNormals = face_normals.size(); - sgWriteInt( fp, nFaceNormals ); - for ( int i = 0; i < nFaceNormals; i++ ) { - sgWritePoint3D( fp, face_normals[i] ); - } - - nFaceAreas = face_areas.size(); - sgWriteInt( fp, nFaceAreas ); - for ( int i = 0; i < nFaceAreas; i++ ) { - sgWriteDouble( fp, face_areas[i] ); - } - - sgWriteString( fp, flag.c_str() ); -} -std::istream& operator>> ( std::istream& input, TGSuperPoly& sp ) -{ - int nFaceNormals; - int nFaceAreas; - Point3D normal; - double area; - - // Load the data - input >> sp.material; - input >> sp.poly; - input >> sp.normals; - input >> sp.texcoords; - input >> sp.tris; - input >> sp.tri_idxs; - - input >> nFaceNormals; - for ( int i = 0; i < nFaceNormals; i++ ) { - input >> normal; - sp.face_normals.push_back(normal); - } - - input >> nFaceAreas; - for ( int i = 0; i < nFaceAreas; i++ ) { - input >> area; - sp.face_areas.push_back(area); - } - - input >> sp.flag; - - return input; -} - -void TGSuperPoly::LoadFromGzFile(gzFile& fp) -{ - int nFaceNormals; - int nFaceAreas; - Point3D normal; - double area; - char* strbuf; - - // Load the data - sgReadString( fp, &strbuf ); - if ( strbuf ) { - material = strbuf; - delete strbuf; - } - - poly.LoadFromGzFile( fp ); - normals.LoadFromGzFile( fp ); - texcoords.LoadFromGzFile( fp ); - - tris.LoadFromGzFile( fp ); - tri_idxs.LoadFromGzFile( fp ); - - sgReadInt( fp, &nFaceNormals ); - for ( int i = 0; i < nFaceNormals; i++ ) { - sgReadPoint3D( fp, normal ); - face_normals.push_back( normal ); - } - - sgReadInt( fp, &nFaceAreas ); - for ( int i = 0; i < nFaceAreas; i++ ) { - sgReadDouble( fp, &area ); - face_areas.push_back( area ); - } - - sgReadString( fp, &strbuf ); - if ( strbuf ) { - flag = strbuf; - delete strbuf; - } -} - -std::ostream& operator<< ( std::ostream& output, const TGPolyNodes& pn ) -{ - int nContours; - int nPoints; - - // Save the data - nContours = pn.poly.size(); - output << nContours << "\n"; - for(int i=0; i> ( std::istream& input, TGPolyNodes& pn ) -{ - int nContours; - int nPoints; - int point; - - // Load the data - input >> nContours; - for(int i=0; i> nPoints; - for (int j=0; j> point; - points.push_back( point ); - } - pn.poly.push_back( points ); - } - - return input; -} - -void TGPolyNodes::LoadFromGzFile(gzFile& fp) -{ - int nContours; - int nPoints; - int point; - - // Load the data - sgReadInt( fp, &nContours ); - for( int i=0; i -#include - -#include -#include - -#include - -// TODO : Needs to be its own class -typedef std::vector < int > int_list; -typedef std::vector < int_list > idx_list; -typedef std::vector < double > double_list; -typedef idx_list::iterator idx_list_iterator; -typedef idx_list::const_iterator const_idx_list_iterator; - -class TGPolyNodes { - -private: - - idx_list poly; // polygon node indexes - -public: - - // Constructor and destructor - TGPolyNodes( void ) {} - ~TGPolyNodes( void ) {} - - // Add a contour - inline void add_contour( const int_list contour ) - { - poly.push_back( contour ); - } - - // Get a contour - inline int_list get_contour( const int i ) const - { - return poly[i]; - } - - // Delete a contour - inline void delete_contour( const int i ) - { - idx_list_iterator start_poly = poly.begin(); - poly.erase( start_poly + i ); - } - - // Add the specified node (index) to the polygon - inline void add_node( int contour, int n ) - { - if ( contour >= (int)poly.size() ) { - // extend polygon - int_list empty_contour; - empty_contour.clear(); - for ( int i = 0; i < contour - (int)poly.size() + 1; ++i ) { - poly.push_back( empty_contour ); - } - } - poly[contour].push_back( n ); - } - - // return size - inline int contours() const - { - return poly.size(); - } - inline int contour_size( int contour ) const - { - return poly[contour].size(); - } - inline int total_size() const - { - int size = 0; - - for ( int i = 0; i < contours(); ++i ) - size += poly[i].size(); - return size; - } - - // return the ith point from the specified contour - inline int get_pt( int contour, int i ) const - { - return poly[contour][i]; - } - - // update the value of a point - inline void set_pt( int contour, int i, const int n ) - { - poly[contour][i] = n; - } - - void SaveToGzFile( gzFile& fp ); - void LoadFromGzFile( gzFile& fp ); - - // Friends for serialization - friend std::istream& operator>> ( std::istream&, TGPolyNodes& ); - friend std::ostream& operator<< ( std::ostream&, const TGPolyNodes& ); -}; -// END TODO - - - -class TGSuperPoly { - -private: - - std::string material; // material/texture name - TGPolygon poly; // master polygon - TGPolygon normals; // corresponding normals - TGPolygon texcoords; // corresponding texture coordinates - - TGPolygon tris; // triangulation - TGPolyNodes tri_idxs; // triangle node indexes - - point_list face_normals; // triangle normals - double_list face_areas; // triangle areas - std::string flag; // For various potential record keeping needs - -public: - - // Constructor and destructor - TGSuperPoly( void ); - ~TGSuperPoly( void ); - - inline std::string get_material() const - { - return material; - } - inline void set_material( const std::string &m ) - { - material = m; - } - - inline TGPolygon get_poly() const - { - return poly; - } - inline void set_poly( const TGPolygon &p ) - { - poly = p; - } - - inline TGPolygon get_normals() const - { - return normals; - } - inline void set_normals( const TGPolygon &p ) - { - normals = p; - } - - inline TGPolygon get_texcoords() const - { - return texcoords; - } - inline void set_texcoords( const TGPolygon &p ) - { - texcoords = p; - } - - inline TGPolygon get_tris() const - { - return tris; - } - inline void set_tris( const TGPolygon &p ) - { - tris = p; - } - - inline TGPolyNodes get_tri_idxs() const - { - return tri_idxs; - } - inline void set_tri_idxs( const TGPolyNodes &p ) - { - tri_idxs = p; - } - - inline Point3D get_face_normal( int tri ) const - { - return face_normals[tri]; - } - - inline point_list get_face_normals() const - { - return face_normals; - } - inline void set_face_normals( const point_list &fns ) - { - face_normals = fns; - } - - inline double get_face_area( int tri ) const - { - return face_areas[tri]; - } - - inline double_list get_face_areas() const - { - return face_areas; - } - inline void set_face_areas( const double_list &fas ) - { - face_areas = fas; - } - - inline std::string get_flag() const - { - return flag; - } - inline void set_flag( const std::string f ) - { - flag = f; - } - - // erase the polygon - void erase(); - - void SaveToGzFile( gzFile& fp ); - void LoadFromGzFile( gzFile& fp ); - - // Friends for serialization - friend std::istream& operator>> ( std::istream&, TGSuperPoly& ); - friend std::ostream& operator<< ( std::ostream&, const TGSuperPoly& ); -}; - - -typedef std::vector < TGSuperPoly > superpoly_list; -typedef superpoly_list::iterator superpoly_list_iterator; -typedef superpoly_list::const_iterator const_superpoly_list_iterator; - - -#endif // _SUPERPOLY_HXX diff --git a/src/Lib/Polygon/texparams.cxx b/src/Lib/Polygon/texparams.cxx deleted file mode 100644 index 433cee3d..00000000 --- a/src/Lib/Polygon/texparams.cxx +++ /dev/null @@ -1,46 +0,0 @@ -#include - -#include "texparams.hxx" - -// Send a TexParam to standard output. -std::ostream& operator << (std::ostream &output, const TGTexParams &tp) -{ - // Save the data - output << tp.ref; - output << tp.width << " "; - output << tp.length << " "; - output << tp.heading << "\n"; - - output << tp.minu << " "; - output << tp.maxu << " "; - output << tp.minv << " "; - output << tp.maxv << "\n"; - - return output; -} - -void TGTexParams::SaveToGzFile(gzFile& fp) -{ - sgWriteGeod( fp, ref ); - sgWriteDouble( fp, width ); - sgWriteDouble( fp, length ); - sgWriteDouble( fp, heading ); - - sgWriteDouble( fp, minu ); - sgWriteDouble( fp, maxu ); - sgWriteDouble( fp, minv ); - sgWriteDouble( fp, maxv ); -} - -void TGTexParams::LoadFromGzFile(gzFile& fp) -{ - sgReadGeod( fp, ref ); - sgReadDouble( fp, &width ); - sgReadDouble( fp, &length ); - sgReadDouble( fp, &heading ); - - sgReadDouble( fp, &minu ); - sgReadDouble( fp, &maxu ); - sgReadDouble( fp, &minv ); - sgReadDouble( fp, &maxv ); -} diff --git a/src/Lib/Polygon/texparams.hxx b/src/Lib/Polygon/texparams.hxx deleted file mode 100644 index 29aa2135..00000000 --- a/src/Lib/Polygon/texparams.hxx +++ /dev/null @@ -1,160 +0,0 @@ -// texparams.hxx -- A simple class to hold texture application -// parameters for sections of the runway -// -// Written by Curtis Olson, started August 2000. -// -// Copyright (C) 2000 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. -// -// $Id: texparams.hxx,v 1.7 2004-11-19 22:25:49 curt Exp $ - - -#ifndef _TEXPARAMS_HXX -#define _TEXPARAMS_HXX - - -#ifndef __cplusplus -# error This library requires C++ -#endif - -#include -#include -#include - -#include -#include - -#include /* for gzFile */ - - -class TGTexParams { - -private: - - SGGeod ref; - double width; - double length; - double heading; - - double minu; - double maxu; - double minv; - double maxv; - -public: - - // Constructor and destructor - inline TGTexParams( void ) - { - } - inline TGTexParams( const SGGeod &r, const double w, const double l, const double h ) - { - ref = r; - width = w; - length = l; - heading = h; - - minu = minv = 0.0; - maxu = maxv = 1.0; - } - inline ~TGTexParams( void ) - { - } - - inline SGGeod get_ref() const - { - return ref; - } - inline void set_ref( const SGGeod &r ) - { - ref = r; - } - - inline double get_width() const - { - return width; - } - inline void set_width( const double w ) - { - width = w; - } - - inline double get_length() const - { - return length; - } - inline void set_length( const double l ) - { - length = l; - } - - inline double get_heading() const - { - return heading; - } - inline void set_heading( const double h ) - { - heading = h; - } - - inline double get_minu() const - { - return minu; - } - inline void set_minu( const double x ) - { - minu = x; - } - - inline double get_maxu() const - { - return maxu; - } - inline void set_maxu( const double x ) - { - maxu = x; - } - - inline double get_minv() const - { - return minv; - } - inline void set_minv( const double x ) - { - minv = x; - } - - inline double get_maxv() const - { - return maxv; - } - inline void set_maxv( const double x ) - { - maxv = x; - } - - void SaveToGzFile( gzFile& fp ); - void LoadFromGzFile( gzFile& fp ); - - // Friend for output - friend std::ostream& operator<< ( std::ostream&, const TGTexParams& ); -}; - -typedef std::vector < TGTexParams > texparams_list; -typedef texparams_list::iterator texparams_list_iterator; -typedef texparams_list::const_iterator const_texparams_list_iterator; - -#endif // _TEXPARAMS_HXX diff --git a/src/Prep/OGRDecode/ogr-decode.cxx b/src/Prep/OGRDecode/ogr-decode.cxx index f860a384..1e983738 100644 --- a/src/Prep/OGRDecode/ogr-decode.cxx +++ b/src/Prep/OGRDecode/ogr-decode.cxx @@ -23,34 +23,27 @@ // // $Id$ - -#include - #include #include +#include + +#include #include #include #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.000001) #define EP_STRETCH (0.1) using std::string; -using std::map; int line_width=50; string line_width_col; @@ -59,7 +52,7 @@ string point_width_col; string area_type="Default"; string area_type_col; int continue_on_errors=0; -int texture_lines = 0; +bool texture_lines = false; int seperate_segments = 0; int max_segment_length=0; // ==0 => don't split int start_record=0; @@ -68,62 +61,27 @@ string attribute_query; bool use_spatial_query=false; double spat_min_x, spat_min_y, spat_max_x, spat_max_y; -void processPoint(OGRPoint* poGeometry, - const string& work_dir, - const string& area_type, - int width) +void processPoint(OGRPoint* poGeometry, const string& area_type, int width, tgChopper& results ) { - TGPolygon shape; - - shape.erase(); - tg::makePolygon(Point3D(poGeometry->getX(),poGeometry->getY(),0), - width, - shape); + SGGeod point = SGGeod::fromDeg( poGeometry->getX(),poGeometry->getY() ); + tgPolygon shape = tgPolygon::Expand( point, width ); if ( max_segment_length > 0 ) { - shape = tgPolygonSplitLongEdges( shape, max_segment_length ); + shape = tgPolygon::SplitLongEdges( shape, max_segment_length ); } - tgChopNormalPolygon(work_dir, area_type, shape, false); + shape.SetPreserve3D( false ); + shape.SetTexMethod( TG_TEX_BY_GEODE ); + + results.Add( shape, area_type ); } -void processLineString(OGRLineString* poGeometry, - const string& work_dir, - const string& area_type, - int width) +void processLineString(OGRLineString* poGeometry, const string& area_type, int width, int with_texture, tgChopper& results ) { - tg::Line line; - TGPolygon shape; + tgpolygon_list segments; + tgContour line; - if (poGeometry->getNumPoints()<2) { - SG_LOG( SG_GENERAL, SG_WARN, "Skipping line with less than two points" ); - return; - } - - shape.erase(); - for (int i=0;igetNumPoints();i++) { - line.addPoint(Point3D(poGeometry->getX(i),poGeometry->getY(i),0)); - } - - tg::makePolygon(line,width,shape); - - if ( max_segment_length > 0 ) { - shape = tgPolygonSplitLongEdges( shape, max_segment_length ); - } - - tgChopNormalPolygon(work_dir, area_type, shape, false); -} - -void processLineStringWithMask(OGRLineString* poGeometry, - const string& work_dir, - const string& area_type, - int width) -{ - poly_list segments; - TGPolygon segment; - tg::Line line; - Point3D p0, p1; + SGGeod p0, p1; double heading, dist, az2; - double pt_x = 0.0f, pt_y = 0.0f; int i, j, numPoints, numSegs; double max_dist; @@ -136,18 +94,17 @@ void processLineStringWithMask(OGRLineString* poGeometry, max_dist = (double)width * 10.0f; // because vector data can generate adjacent polys, lets stretch the two enpoints by a little bit - p0 = Point3D(poGeometry->getX(0),poGeometry->getY(0),0); - p1 = Point3D(poGeometry->getX(1),poGeometry->getY(1),0); + p0 = SGGeod::fromDeg( poGeometry->getX(0), poGeometry->getY(0) ); + p1 = SGGeod::fromDeg( poGeometry->getX(1), poGeometry->getY(1) ); - geo_inverse_wgs_84( p1.y(), p1.x(), p0.y(), p0.x(), &heading, &az2, &dist); - geo_direct_wgs_84( p0.y(), p0.x(), heading, EP_STRETCH, &pt_y, &pt_x, &az2 ); - line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); + heading = SGGeodesy::courseDeg( p1, p0 ); + line.AddNode( SGGeodesy::direct( p0, heading, EP_STRETCH ) ); // now add the middle points : if they are too far apart, add intermediate nodes for ( i=1;igetX(i-1),poGeometry->getY(i-1),0); - p1 = Point3D(poGeometry->getX(i ),poGeometry->getY(i ),0); - geo_inverse_wgs_84( p0.y(), p0.x(), p1.y(), p1.x(), &heading, &az2, &dist); + p0 = SGGeod::fromDeg( poGeometry->getX(i-1), poGeometry->getY(i-1) ); + p1 = SGGeod::fromDeg( poGeometry->getX(i ), poGeometry->getY(i ) ); + SGGeodesy::inverse( p0, p1, heading, az2, dist ); if (dist > max_dist) { @@ -156,151 +113,55 @@ void processLineStringWithMask(OGRLineString* poGeometry, for (j=0; jgetX(numPoints-2),poGeometry->getY(numPoints-2),0); - p1 = Point3D(poGeometry->getX(numPoints-1),poGeometry->getY(numPoints-1),0); - - geo_inverse_wgs_84( p0.y(), p0.x(), p1.y(), p1.x(), &heading, &az2, &dist); - geo_direct_wgs_84( p1.y(), p1.x(), heading, EP_STRETCH, &pt_y, &pt_x, &az2 ); - line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); - - // make a plygons from the line segments - tg::makePolygons(line,width,segments); - - tgChopNormalPolygonsWithMask(work_dir, area_type, segments, false); -} - -void processLineStringWithTextureInfo(OGRLineString* poGeometry, - const string& work_dir, - const string& area_type, - int width) -{ - poly_list segments; - texparams_list tps; - tg::Line line; - Point3D p0, p1; - double heading, dist, az2; - double pt_x = 0.0f, pt_y = 0.0f; - int i, j, numPoints, numSegs; - double max_dist; -// double min_dist; -// double cur_dist; - - numPoints = poGeometry->getNumPoints(); - if (numPoints < 2) { - SG_LOG( SG_GENERAL, SG_WARN, "Skipping line with less than two points" ); - return; - } - - max_dist = (double)width * 10.0f; -// min_dist = (double)width * 1.5f; -// cur_dist = 0.0f; - - // because vector data can generate adjacent polys, lets stretch the two enpoints by a little bit - p0 = Point3D(poGeometry->getX(0),poGeometry->getY(0),0); - p1 = Point3D(poGeometry->getX(1),poGeometry->getY(1),0); - - geo_inverse_wgs_84( p1.y(), p1.x(), p0.y(), p0.x(), &heading, &az2, &dist); - geo_direct_wgs_84( p0.y(), p0.x(), heading, EP_STRETCH, &pt_y, &pt_x, &az2 ); - line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); - - // now add the middle points : if they are too far apart, add intermediate nodes - for ( i=1;igetX(i-1),poGeometry->getY(i-1),0); - p1 = Point3D(poGeometry->getX(i ),poGeometry->getY(i ),0); - geo_inverse_wgs_84( p0.y(), p0.x(), p1.y(), p1.x(), &heading, &az2, &dist); - - if (dist > max_dist) - { - numSegs = (dist / max_dist) + 1; - dist = dist / (double)numSegs; - - for (j=0; jgetX(numPoints-2),poGeometry->getY(numPoints-2),0); - p1 = Point3D(poGeometry->getX(numPoints-1),poGeometry->getY(numPoints-1),0); + p0 = SGGeod::fromDeg( poGeometry->getX(numPoints-2), poGeometry->getY(numPoints-2) ); + p1 = SGGeod::fromDeg( poGeometry->getX(numPoints-1), poGeometry->getY(numPoints-1) ); - geo_inverse_wgs_84( p0.y(), p0.x(), p1.y(), p1.x(), &heading, &az2, &dist); - geo_direct_wgs_84( p1.y(), p1.x(), heading, EP_STRETCH, &pt_y, &pt_x, &az2 ); - line.addPoint( Point3D( pt_x, pt_y, 0.0f ) ); + heading = SGGeodesy::courseDeg( p0, p1 ); + line.AddNode( SGGeodesy::direct(p1, heading, EP_STRETCH) ); // make a plygons from the line segments - tg::makePolygonsTP(line,width,segments,tps); + segments = tgContour::ExpandToPolygons( line, width ); + for ( unsigned int i=0; igetGeometryType()&wkb25DBit)==wkb25DBit); - TGPolygon shape; - - shape.erase(); - - OGRLinearRing *ring; + // bool preserve3D = ((poGeometry->getGeometryType()&wkb25DBit)==wkb25DBit); // first add the outer ring - ring=poGeometry->getExteriorRing(); - for (int i=0;igetNumPoints();i++) { - shape.add_node(0, - Point3D(ring->getX(i), - ring->getY(i), - ring->getZ(i))); - } - shape.set_hole_flag( 0, 0 ); - - // then add the inner rings - for ( int j = 0 ; j < poGeometry->getNumInteriorRings() ; j ++ ) { - ring=poGeometry->getInteriorRing( j ); - for (int i=0;igetNumPoints();i++) { - shape.add_node(j+1, - Point3D(ring->getX(i), - ring->getY(i), - ring->getZ(i))); - } - shape.set_hole_flag( j+1 , 1 ); - } + tgPolygon shape = tgPolygon::FromOGR( poGeometry ); if ( max_segment_length > 0 ) { - shape = tgPolygonSplitLongEdges( shape, max_segment_length ); + shape = tgPolygon::SplitLongEdges( shape, max_segment_length ); } - tgChopNormalPolygon(work_dir, area_type, shape, preserve3D); + // shape.SetPreserve3D( preserve3D ); + shape.SetPreserve3D( false ); + + results.Add( shape, area_type ); } -void processLayer(OGRLayer* poLayer, - const string& work_dir) { - +void processLayer(OGRLayer* poLayer, tgChopper& results ) +{ int feature_count=poLayer->GetFeatureCount(); if (feature_count!=-1 && start_record>0 && start_record>=feature_count) { @@ -388,9 +249,12 @@ void processLayer(OGRLayer* poLayer, } } + // PSADRO TODO : Generate the work queue for this layer + // Then process the workqueue with threads OGRFeature *poFeature; - poLayer->SetNextByIndex(start_record); + int numFeatures = poLayer->GetFeatureCount(); + int cur_feature = 1; for ( ; (poFeature = poLayer->GetNextFeature()) != NULL; OGRFeature::DestroyFeature( poFeature ) ) { OGRGeometry *poGeometry; @@ -423,29 +287,29 @@ void processLayer(OGRLayer* poLayer, area_type_name=poFeature->GetFieldAsString(area_type_field); } - if ( is_ocean_area(area_type_name) ) { - // interior of polygon is ocean, holes are islands + if ( is_ocean_area(area_type_name) ) { + // interior of polygon is ocean, holes are islands - SG_LOG( SG_GENERAL, SG_ALERT, "Ocean area ... SKIPPING!" ); + SG_LOG( SG_GENERAL, SG_ALERT, "Ocean area ... SKIPPING!" ); - // Ocean data now comes from GSHHS so we want to ignore - // all other ocean data - continue; - } else if ( is_void_area(area_type_name) ) { - // interior is ???? + // Ocean data now comes from GSHHS so we want to ignore + // all other ocean data + continue; + } else if ( is_void_area(area_type_name) ) { + // interior is ???? - // skip for now - SG_LOG( SG_GENERAL, SG_ALERT, "Void area ... SKIPPING!" ); + // skip for now + SG_LOG( SG_GENERAL, SG_ALERT, "Void area ... SKIPPING!" ); - continue; - } else if ( is_null_area(area_type_name) ) { - // interior is ???? + continue; + } else if ( is_null_area(area_type_name) ) { + // interior is ???? - // skip for now - SG_LOG( SG_GENERAL, SG_ALERT, "Null area ... SKIPPING!" ); + // skip for now + SG_LOG( SG_GENERAL, SG_ALERT, "Null area ... SKIPPING!" ); - continue; - } + continue; + } poGeometry->transform( poCT ); @@ -459,7 +323,7 @@ void processLayer(OGRLayer* poLayer, width=point_width; } } - processPoint((OGRPoint*)poGeometry,work_dir,area_type_name,width); + processPoint((OGRPoint*)poGeometry, area_type_name, width, results); break; } case wkbMultiPoint: { @@ -473,7 +337,7 @@ void processLayer(OGRLayer* poLayer, } OGRMultiPoint* multipt=(OGRMultiPoint*)poGeometry; for (int i=0;igetNumGeometries();i++) { - processPoint((OGRPoint*)(multipt->getGeometryRef(i)),work_dir,area_type_name,width); + processPoint((OGRPoint*)(multipt->getGeometryRef(i)), area_type_name, width, results); } break; } @@ -487,13 +351,7 @@ void processLayer(OGRLayer* poLayer, } } - if (texture_lines) { - processLineStringWithTextureInfo((OGRLineString*)poGeometry,work_dir,area_type_name,width); - } else if (seperate_segments) { - processLineStringWithMask((OGRLineString*)poGeometry,work_dir,area_type_name,width); - } else { - processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); - } + processLineString((OGRLineString*)poGeometry, area_type_name, width, texture_lines, results); break; } case wkbMultiLineString: { @@ -508,26 +366,20 @@ void processLayer(OGRLayer* poLayer, OGRMultiLineString* multils=(OGRMultiLineString*)poGeometry; for (int i=0;igetNumGeometries();i++) { - if (texture_lines) { - processLineStringWithTextureInfo((OGRLineString*)poGeometry,work_dir,area_type_name,width); - } else if (seperate_segments) { - processLineStringWithMask((OGRLineString*)poGeometry,work_dir,area_type_name,width); - } else { - processLineString((OGRLineString*)poGeometry,work_dir,area_type_name,width); - } + processLineString((OGRLineString*)poGeometry, area_type_name, width, texture_lines, results); } break; } case wkbPolygon: { SG_LOG( SG_GENERAL, SG_DEBUG, "Polygon feature" ); - processPolygon((OGRPolygon*)poGeometry,work_dir,area_type_name); + processPolygon((OGRPolygon*)poGeometry, area_type_name, results); break; } case wkbMultiPolygon: { SG_LOG( SG_GENERAL, SG_DEBUG, "MultiPolygon feature" ); OGRMultiPolygon* multipoly=(OGRMultiPolygon*)poGeometry; for (int i=0;igetNumGeometries();i++) { - processPolygon((OGRPolygon*)(multipoly->getGeometryRef(i)),work_dir,area_type_name); + processPolygon((OGRPolygon*)(multipoly->getGeometryRef(i)), area_type_name, results); } break; } @@ -589,94 +441,105 @@ int main( int argc, char **argv ) { sglog().setLogLevels( SG_ALL, SG_INFO ); while (argc>1) { - if (!strcmp(argv[1],"--line-width")) { - if (argc<3) - usage(progname); - line_width=atoi(argv[2]); - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--line-width-column")) { - if (argc<3) - usage(progname); - line_width_col=argv[2]; - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--point-width")) { - if (argc<3) - usage(progname); - point_width=atoi(argv[2]); - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--point-width-column")) { - if (argc<3) - usage(progname); - point_width_col=argv[2]; - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--area-type")) { - if (argc<3) - usage(progname); - area_type=argv[2]; - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--area-type-column")) { - if (argc<3) - usage(progname); - area_type_col=argv[2]; - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--texture-lines")) { - argv++; - argc--; - texture_lines=1; - } else if (!strcmp(argv[1],"--seperate-segments")) { - argv++; - argc--; - seperate_segments=1; - } else if (!strcmp(argv[1],"--continue-on-errors")) { + if (!strcmp(argv[1],"--line-width")) { + if (argc<3) { + usage(progname); + } + line_width=atoi(argv[2]); + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--line-width-column")) { + if (argc<3) { + usage(progname); + } + line_width_col=argv[2]; + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--point-width")) { + if (argc<3) { + usage(progname); + } + point_width=atoi(argv[2]); + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--point-width-column")) { + if (argc<3) { + usage(progname); + } + point_width_col=argv[2]; + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--area-type")) { + if (argc<3) { + usage(progname); + } + area_type=argv[2]; + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--area-type-column")) { + if (argc<3) { + usage(progname); + } + area_type_col=argv[2]; + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--texture-lines")) { argv++; argc--; - continue_on_errors=1; - } else if (!strcmp(argv[1],"--max-segment")) { - if (argc<3) - usage(progname); - max_segment_length=atoi(argv[2]); - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--start-record")) { - if (argc<3) - usage(progname); - start_record=atoi(argv[2]); - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--where")) { - if (argc<3) - usage(progname); + texture_lines=true; + } else if (!strcmp(argv[1],"--seperate-segments")) { + argv++; + argc--; + seperate_segments=1; + } else if (!strcmp(argv[1],"--continue-on-errors")) { + argv++; + argc--; + continue_on_errors=1; + } else if (!strcmp(argv[1],"--max-segment")) { + if (argc<3) { + usage(progname); + } + max_segment_length=atoi(argv[2]); + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--start-record")) { + if (argc<3) { + usage(progname); + } + start_record=atoi(argv[2]); + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--where")) { + if (argc<3) { + usage(progname); + } use_attribute_query=true; - attribute_query=argv[2]; - argv+=2; - argc-=2; - } else if (!strcmp(argv[1],"--spat")) { - if (argc<6) - usage(progname); + attribute_query=argv[2]; + argv+=2; + argc-=2; + } else if (!strcmp(argv[1],"--spat")) { + if (argc<6) { + usage(progname); + } use_spatial_query=true; spat_min_x=atof(argv[2]); spat_min_y=atof(argv[3]); spat_max_x=atof(argv[4]); spat_max_y=atof(argv[5]); - argv+=5; - argc-=5; - } else if (!strcmp(argv[1],"--help")) { - usage(progname); - } else - break; + argv+=5; + argc-=5; + } else if (!strcmp(argv[1],"--help")) { + usage(progname); + } else { + break; + } } SG_LOG( SG_GENERAL, SG_ALERT, "ogr-decode version " << getTGVersion() << "\n" ); - if (argc<3) - usage(progname); - + if (argc<3) { + usage(progname); + } work_dir=argv[1]; datasource=argv[2]; @@ -684,14 +547,19 @@ int main( int argc, char **argv ) { sgp.append( "dummy" ); sgp.create_dir( 0755 ); + tgChopper results( work_dir ); + // initialize persistant polygon counter - string counter_file = work_dir + "/poly_counter"; - poly_index_init( counter_file ); + //string counter_file = work_dir + "/poly_counter"; + //poly_index_init( counter_file ); + + // new chop api + //std::string counter_file2 = work_dir + "/poly_counter2"; + //tgPolygon::ChopIdxInit( counter_file ); SG_LOG( SG_GENERAL, SG_DEBUG, "Opening datasource " << datasource << " for reading." ); OGRRegisterAll(); - OGRDataSource *poDS; poDS = OGRSFDriverRegistrar::Open( datasource.c_str(), FALSE ); @@ -704,27 +572,28 @@ int main( int argc, char **argv ) { OGRLayer *poLayer; if (argc>3) { - for (int i=3;iGetLayerByName( argv[i] ); + for (int i=3;iGetLayerByName( argv[i] ); - if (poLayer == NULL ) - { - SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource ); - exit( 1 ); - } - - processLayer(poLayer, work_dir); - } + if (poLayer == NULL ) + { + SG_LOG( SG_GENERAL, SG_ALERT, "Failed opening layer " << argv[i] << " from datasource " << datasource ); + exit( 1 ); + } + processLayer(poLayer, results ); + } } else { - for (int i=0;iGetLayerCount();i++) { - poLayer = poDS->GetLayer(i); + for (int i=0;iGetLayerCount();i++) { + poLayer = poDS->GetLayer(i); - assert(poLayer != NULL); + assert(poLayer != NULL); - processLayer(poLayer, work_dir); - } + processLayer(poLayer, results ); + } } + results.Save(); + OGRDataSource::DestroyDataSource( poDS ); return 0;