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;