From 071d98b3bc25c808d031715c3d7b12cbe2272a6f Mon Sep 17 00:00:00 2001 From: Peter Sadrozinski Date: Tue, 14 Feb 2012 18:08:11 -0500 Subject: [PATCH] Lot's pf success on the quest to a full file parse. Just three airports crashed, and no hangs. One of the crashes is still NZSP I don't know if I'll fix this soon. It's the south pole problem with stuff crossing the IDL. - Added dome debug support for generating shapefiles from TGPolygon. This allows directly creating shapefiles instead of chopping a poly, then having to run poly2ogr on it. The API is flexible, allowing multiple datafiles, and layers open simutaneously, and the features in the shapefiles have a description field that can be set. - Added a restart-id= which allows a bash script to check the exit value of genapts, and automatically restart on the airport after the one that crashed. I've included an example script I have used. - Stability fixes: 1) actually, this is both stability, and graphical fix. Some airport linear features and polygon contours were using bezier nodes for straight lines. When the control point lies just outside the line, the parser would attempt to draw a 'hook' unfortunately, some of my calcs return NAN, and sometimes you just get crazy looking polys to attempt to clip - not good. Fix is to check cubic and quadratic nodes to see if they should really be linear. Substitute linear segments where apropriate. 2) Along the same line as 1, There are some really short bezier segments, that were getting subdivided into the hard coded 8 segments. Now, if the approximately linear distance of the curve is less than 4 meters, the curve is broken into 1/2 meter segments. Likewise, if the curve is greater than 800 meters, the curve is broken into 100 meter segments. 3) Sometimes, during clipping, degenerate ontours are created. Although these are cleaned before triangulation, there were still issues with generating the accumulation poly with the union operator, causing a crash. I now simplify, and reduce degeneracy on all polys before clipping against the accumulator. - Extra special bonus - successful completeion of an airport is printed in green :) --- src/Airports/GenAirports/CMakeLists.txt | 1 + src/Airports/GenAirports850/CMakeLists.txt | 1 + src/Airports/GenAirports850/airport.cxx | 51 ++- src/Airports/GenAirports850/airport.hxx | 5 + src/Airports/GenAirports850/beznode.hxx | 17 +- src/Airports/GenAirports850/closedpoly.cxx | 208 ++++++---- src/Airports/GenAirports850/linearfeature.cxx | 273 +++++++++---- src/Airports/GenAirports850/linearfeature.hxx | 4 +- src/Airports/GenAirports850/main.cxx | 43 +- src/Airports/GenAirports850/parser.cxx | 22 +- src/Airports/GenAirports850/parser.hxx | 3 +- src/Airports/GenAirports850/runway.cxx | 2 + src/Airports/GenAirports850/runway.hxx | 5 +- src/Airports/GenAirports850/rwy_gen.cxx | 383 +++++++++++------- src/BuildTiles/Main/CMakeLists.txt | 1 + src/Lib/Geometry/CMakeLists.txt | 4 +- src/Lib/Geometry/poly_support.cxx | 152 ++++++- src/Lib/Geometry/poly_support.hxx | 10 + 18 files changed, 845 insertions(+), 340 deletions(-) diff --git a/src/Airports/GenAirports/CMakeLists.txt b/src/Airports/GenAirports/CMakeLists.txt index f00c83df..045435be 100644 --- a/src/Airports/GenAirports/CMakeLists.txt +++ b/src/Airports/GenAirports/CMakeLists.txt @@ -24,6 +24,7 @@ target_link_libraries(genapts TriangleJRS ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} + ${GDAL_LIBRARY} ${GPC_LIBRARY} ${NEWMAT_LIBRARY}) diff --git a/src/Airports/GenAirports850/CMakeLists.txt b/src/Airports/GenAirports850/CMakeLists.txt index 9791bf19..4655c681 100644 --- a/src/Airports/GenAirports850/CMakeLists.txt +++ b/src/Airports/GenAirports850/CMakeLists.txt @@ -25,6 +25,7 @@ target_link_libraries(genapts850 Polygon Geometry Array Optimize Output poly2tri TriangleJRS + ${GDAL_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} ${GPC_LIBRARY} diff --git a/src/Airports/GenAirports850/airport.cxx b/src/Airports/GenAirports850/airport.cxx index cbfb2eff..c8872009 100644 --- a/src/Airports/GenAirports850/airport.cxx +++ b/src/Airports/GenAirports850/airport.cxx @@ -29,6 +29,12 @@ #include "elevations.hxx" + +#include + +string SGLOG_GREEN = "\033[0;32m"; +string SGLOG_NORMAL = "\033[0m"; + Airport::Airport( int c, char* def) { int numParams; @@ -426,6 +432,14 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) } } +#if 0 + if (icao == "SSOK") + { + SG_LOG(SG_GENERAL, SG_ALERT, "Injecting error at icao " << icao ); + exit(-1); + } +#endif + // Starting to clip the polys gettimeofday(&build_start, NULL); @@ -435,7 +449,17 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) for ( unsigned int i=0; iGetDescription() ); - features[i]->BuildBtg( altitude, &line_polys, &line_tps, &line_accum, &rwy_lights ); + +#if 0 + if ( i == 22 ) { + features[i]->BuildBtg( altitude, &line_polys, &line_tps, &line_accum, &rwy_lights, true ); + } else { + //features[i]->BuildBtg( altitude, &line_polys, &line_tps, &line_accum, &rwy_lights, false ); + } +#else + features[i]->BuildBtg( altitude, &line_polys, &line_tps, &line_accum, &rwy_lights, false ); +#endif + } } else @@ -526,6 +550,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) gettimeofday(&log_time, NULL); SG_LOG( SG_GENERAL, SG_ALERT, "Finished building pavements for " << icao << " at " << ctime(&log_time.tv_sec) ); +#if 0 // Build runway shoulders here for ( unsigned int i=0; iGetsShoulder() ) { - runways[i]->BuildShoulder( altitude, &rwy_polys, &rwy_tps, &accum ); + runways[i]->BuildShoulder( altitude, &rwy_polys, &rwy_tps, &accum, &apt_base, &apt_clearing ); } } +#endif // build the base and clearing if there's a boundary if (boundary) @@ -671,6 +697,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) poly = remove_spikes( poly ); poly = remove_dups( poly ); poly = remove_bad_contours( poly ); + poly = remove_small_contours( poly ); rwy_polys[k].set_poly( poly ); } @@ -688,6 +715,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) poly = remove_spikes( poly ); poly = remove_dups( poly ); poly = remove_bad_contours( poly ); + poly = remove_small_contours( poly ); pvmt_polys[k].set_poly( poly ); } @@ -701,10 +729,12 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) poly = remove_cycles( poly ); poly = remove_dups( poly ); poly = remove_bad_contours( poly ); + poly = tgPolygonSimplify( poly ); poly = remove_tiny_contours( poly ); poly = remove_spikes( poly ); poly = remove_dups( poly ); poly = remove_bad_contours( poly ); + poly = remove_tiny_contours( poly ); line_polys[k].set_poly( poly ); } @@ -728,6 +758,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) base_poly = remove_spikes( base_poly ); base_poly = remove_dups( base_poly ); base_poly = remove_bad_contours( base_poly ); + base_poly = remove_small_contours( base_poly ); gettimeofday(&cleanup_end, NULL); timersub(&cleanup_end, &cleanup_start, &cleanup_time); @@ -778,7 +809,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) TGPolygon poly = pvmt_polys[i].get_poly(); #if 0 - if ( i == 1 ) { + if ( i == 0 ) { SG_LOG(SG_GENERAL, SG_INFO, "Problem poly: " << poly ); tgChopNormalPolygon( "/home/pete", "Base", poly, false ); @@ -814,10 +845,10 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) TGPolygon poly = line_polys[i].get_poly(); #if 0 - if ( i == 627 ) { + if ( i == 282 ) { SG_LOG(SG_GENERAL, SG_INFO, "Problem poly: " << poly ); - tgChopNormalPolygon( "/home/pete", "Base", poly, false ); + tgChopNormalPolygon( "/home/pete/", "Base", poly, false ); verbose_triangulation = true; } else { verbose_triangulation = false; @@ -839,7 +870,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) { SG_LOG(SG_GENERAL, SG_INFO, "Problem poly: " << base_poly ); - tgChopNormalPolygon( "/home/pete", "Base", base_poly, false ); + tgChopNormalPolygon( "/home/pete/", "Base", base_poly, false ); verbose_triangulation = true; } #endif @@ -1099,7 +1130,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) // Now build the fitted airport surface ... // calculation min/max coordinates of airport area - SG_LOG(SG_GENERAL, SG_INFO, " calculation min/max coordinates of airport area"); + SG_LOG(SG_GENERAL, SG_DEBUG, " calculation min/max coordinates of airport area"); Point3D min_deg(9999.0, 9999.0, 0), max_deg(-9999.0, -9999.0, 0); for ( unsigned int j = 0; j < nodes.get_node_list().size(); ++j ) @@ -1127,7 +1158,7 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) } } - SG_LOG(SG_GENERAL, SG_INFO, "Before extending for lights: min = " << min_deg << " max = " << max_deg ); + SG_LOG(SG_GENERAL, SG_DEBUG, "Before extending for lights: min = " << min_deg << " max = " << max_deg ); // extend the min/max coordinates of airport area to cover all // lights as well @@ -1159,8 +1190,6 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) } } - SG_LOG(SG_GENERAL, SG_DEBUG, " done " ); - // Extend the area a bit so we don't have wierd things on the edges double dlon = max_deg.lon() - min_deg.lon(); double dlat = max_deg.lat() - min_deg.lat(); @@ -1571,4 +1600,6 @@ void Airport::BuildBtg(const string& root, const string_list& elev_src ) // write_boundary( holepath, b, hull, poly_index ); tgChopNormalPolygon( holepath, "Hole", divided_base, true ); tgChopNormalPolygon( holepath, "Airport", apt_clearing, false ); + + SG_LOG( SG_GENERAL, SG_ALERT, SGLOG_GREEN << "\nSUCCESS generating " << icao << SGLOG_NORMAL << "\n" ); } diff --git a/src/Airports/GenAirports850/airport.hxx b/src/Airports/GenAirports850/airport.hxx index 37fe5ad2..aaa851a5 100644 --- a/src/Airports/GenAirports850/airport.hxx +++ b/src/Airports/GenAirports850/airport.hxx @@ -58,6 +58,11 @@ public: } } + int NumFeatures( void ) + { + return features.size(); + } + void SetBoundary( ClosedPoly* bndry ) { boundary = bndry; diff --git a/src/Airports/GenAirports850/beznode.hxx b/src/Airports/GenAirports850/beznode.hxx index 257c9972..fc9f3e2e 100644 --- a/src/Airports/GenAirports850/beznode.hxx +++ b/src/Airports/GenAirports850/beznode.hxx @@ -9,6 +9,12 @@ #include // TEMP... +inline double LinearDistance( Point3D p0, Point3D p1 ) +{ + return sqrt( pow( fabs( p1.x()-p0.x() ), 2 ) + + pow( fabs( p1.y()-p0.y() ), 2 ) ); +} + inline Point3D CalculateLinearLocation( Point3D p0, Point3D p1, double t ) { Point3D result; @@ -21,6 +27,11 @@ inline Point3D CalculateLinearLocation( Point3D p0, Point3D p1, double t ) return result; } +inline double QuadraticDistance( Point3D p0, Point3D cp, Point3D p1 ) +{ + return LinearDistance( p0, cp ) + LinearDistance( cp, p1 ); +} + inline Point3D CalculateQuadraticLocation( Point3D p0, Point3D cp, Point3D p1, double t ) { Point3D result; @@ -34,7 +45,11 @@ inline Point3D CalculateQuadraticLocation( Point3D p0, Point3D cp, Point3D p1, d return result; } -// TODO: Should be in a math library +inline double CubicDistance( Point3D p0, Point3D cp0, Point3D cp1, Point3D p1 ) +{ + return LinearDistance( p0, cp0 ) + LinearDistance( cp0, cp1 ) + LinearDistance( cp1, p1 ); +} + inline Point3D CalculateCubicLocation( Point3D p0, Point3D cp0, Point3D cp1, Point3D p1, double t ) { Point3D result; diff --git a/src/Airports/GenAirports850/closedpoly.cxx b/src/Airports/GenAirports850/closedpoly.cxx index fd8d7f6e..a57c1042 100644 --- a/src/Airports/GenAirports850/closedpoly.cxx +++ b/src/Airports/GenAirports850/closedpoly.cxx @@ -100,6 +100,7 @@ void ClosedPoly::AddNode( BezNode* node ) } } +#if 0 void ClosedPoly::CreateConvexHull( void ) { TGPolygon convexHull; @@ -122,6 +123,7 @@ void ClosedPoly::CreateConvexHull( void ) SG_LOG(SG_GENERAL, SG_ALERT, "Boundary size too small: " << boundary->size() << ". Ignoring..." ); } } +#endif void ClosedPoly::CloseCurContour() { @@ -132,7 +134,7 @@ void ClosedPoly::CloseCurContour() if (cur_feature) { SG_LOG(SG_GENERAL, SG_DEBUG, "We still have an active linear feature - add the first node to close it"); - cur_feature->Finish(true); + cur_feature->Finish(true, features.size() ); features.push_back(cur_feature); cur_feature = NULL; @@ -145,7 +147,7 @@ void ClosedPoly::CloseCurContour() boundary = cur_contour; // generate the convex hull from the bezcontour node locations - CreateConvexHull(); + // CreateConvexHull(); cur_contour = NULL; } @@ -158,18 +160,18 @@ void ClosedPoly::CloseCurContour() void ClosedPoly::ConvertContour( BezContour* src, point_list *dst ) { - BezNode* prevNode; BezNode* curNode; BezNode* nextNode; - Point3D prevLoc; Point3D curLoc; Point3D nextLoc; Point3D cp1; Point3D cp2; int curve_type = CURVE_LINEAR; - unsigned int i; + double total_dist; + double meter_dist = 1.0f/96560.64f; + int num_segs = BEZIER_DETAIL; SG_LOG(SG_GENERAL, SG_DEBUG, "Creating a contour with " << src->size() << " nodes"); @@ -177,23 +179,11 @@ void ClosedPoly::ConvertContour( BezContour* src, point_list *dst ) dst->empty(); // iterate through each bezier node in the contour - for (i=0; i <= src->size()-1; i++) + for (unsigned int i = 0; i <= src->size()-1; i++) { SG_LOG(SG_GENERAL, SG_DEBUG, "\nHandling Node " << i << "\n\n"); - if (i == 0) - { - // set prev node to last in the contour, as all contours must be closed - prevNode = src->at( src->size()-1 ); - } - else - { - // otherwise, it's just the previous index - prevNode = src->at( i-1 ); - } - curNode = src->at(i); - if (i < src->size() - 1) { nextNode = src->at(i+1); @@ -204,37 +194,6 @@ void ClosedPoly::ConvertContour( BezContour* src, point_list *dst ) nextNode = src->at(0); } - // determine the type of curve from prev (just to get correct prev location) - // once we start drawing the curve from cur to next, we can just remember the prev loc - if (prevNode->HasNextCp()) - { - // curve from prev is cubic or quadratic - if(curNode->HasPrevCp()) - { - // curve from prev is cubic : calculate the last location on the curve - prevLoc = CalculateCubicLocation( prevNode->GetLoc(), prevNode->GetNextCp(), curNode->GetPrevCp(), curNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (BEZIER_DETAIL-1) ); - } - else - { - // curve from prev is quadratic : use prev node next cp - prevLoc = CalculateQuadraticLocation( prevNode->GetLoc(), prevNode->GetNextCp(), curNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (BEZIER_DETAIL-1) ); - } - } - else - { - // curve from prev is quadratic or linear - if( curNode->HasPrevCp() ) - { - // curve from prev is quadratic : calculate the last location on the curve - prevLoc = CalculateQuadraticLocation( prevNode->GetLoc(), curNode->GetPrevCp(), curNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (BEZIER_DETAIL-1) ); - } - else - { - // curve from prev is linear : just use prev node location - prevLoc = prevNode->GetLoc(); - } - } - // now determine how we will iterate from current node to next node if( curNode->HasNextCp() ) { @@ -245,12 +204,14 @@ void ClosedPoly::ConvertContour( BezContour* src, point_list *dst ) curve_type = CURVE_CUBIC; cp1 = curNode->GetNextCp(); cp2 = nextNode->GetPrevCp(); + total_dist = CubicDistance( curNode->GetLoc(), cp1, cp2, nextNode->GetLoc() ); } else { // curve is quadratic using current nodes cp as the cp curve_type = CURVE_QUADRATIC; cp1 = curNode->GetNextCp(); + total_dist = QuadraticDistance( curNode->GetLoc(), cp1, nextNode->GetLoc() ); } } else @@ -261,28 +222,94 @@ void ClosedPoly::ConvertContour( BezContour* src, point_list *dst ) // curve is quadratic using next nodes cp as the cp curve_type = CURVE_QUADRATIC; cp1 = nextNode->GetPrevCp(); + total_dist = QuadraticDistance( curNode->GetLoc(), cp1, nextNode->GetLoc() ); } else { // curve is linear curve_type = CURVE_LINEAR; + total_dist = LinearDistance( curNode->GetLoc(), nextNode->GetLoc() ); } } +#if 0 + double num_meters = total_dist / meter_dist; + + // If total distance is < 4 meters, then we need to modify num Segments so that each segment >= 0.5 meter + if (num_meters < 4.0f) + { + num_segs = ((int)num_meters + 1) * 2; + } + else if (num_meters > 800.0f) + { + num_segs = num_meters / 100.0f + 1; + } + else + { + num_segs = 8; + } +#endif + + double num_meters = total_dist / meter_dist; + if (num_meters < 4.0f) + { + if (curve_type != CURVE_LINEAR) + { + // If total distance is < 4 meters, then we need to modify num Segments so that each segment >= 1/2 meter + num_segs = ((int)num_meters + 1) * 2; + SG_LOG(SG_GENERAL, SG_DEBUG, "Segment from (" << curNode->GetLoc().x() << "," << curNode->GetLoc().y() << ") to (" << nextNode->GetLoc().x() << "," << nextNode->GetLoc().y() << ")" ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Distance is " << num_meters << " ( < 4.0) so num_segs is " << num_segs ); + } + else + { + num_segs = 1; + } + } + else if (num_meters > 800.0f) + { + // If total distance is > 800 meters, then we need to modify num Segments so that each segment <= 100 meters + num_segs = num_meters / 100.0f + 1; + SG_LOG(SG_GENERAL, SG_DEBUG, "Segment from (" << curNode->GetLoc().x() << "," << curNode->GetLoc().y() << ") to (" << nextNode->GetLoc().x() << "," << nextNode->GetLoc().y() << ")" ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Distance is " << num_meters << " ( > 100.0) so num_segs is " << num_segs ); + } + else + { + if (curve_type != CURVE_LINEAR) + { + num_segs = 8; + // num_segs = 16; + SG_LOG(SG_GENERAL, SG_DEBUG, "Segment from (" << curNode->GetLoc().x() << "," << curNode->GetLoc().y() << ") to (" << nextNode->GetLoc().x() << "," << nextNode->GetLoc().y() << ")" ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Distance is " << num_meters << " (OK) so num_segs is " << num_segs ); + } + else + { + // make sure linear segments don't got over 100m + //num_segs = 1; + num_segs = num_meters / 100.0f + 1; + } + } + + // if only one segment, revert to linear + if (num_segs == 1) + { + curve_type = CURVE_LINEAR; + } + + // initialize current location curLoc = curNode->GetLoc(); if (curve_type != CURVE_LINEAR) { - for (int p=0; pGetLoc(), cp1, nextNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (p+1) ); + nextLoc = CalculateQuadraticLocation( curNode->GetLoc(), cp1, nextNode->GetLoc(), (1.0f/num_segs) * (p+1) ); } else { - nextLoc = CalculateCubicLocation( curNode->GetLoc(), cp1, cp2, nextNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (p+1) ); + nextLoc = CalculateCubicLocation( curNode->GetLoc(), cp1, cp2, nextNode->GetLoc(), (1.0f/num_segs) * (p+1) ); } // add the pavement vertex @@ -298,22 +325,57 @@ void ClosedPoly::ConvertContour( BezContour* src, point_list *dst ) SG_LOG(SG_GENERAL, SG_DEBUG, " add bezier node (type " << curve_type << ") at (" << curLoc.x() << "," << curLoc.y() << ")"); } - // now set set prev and cur locations for the next iteration - prevLoc = curLoc; + // now set set cur location for the next iteration curLoc = nextLoc; } } else { - nextLoc = nextNode->GetLoc(); +// nextLoc = nextNode->GetLoc(); - // just add the one vertex - linear - dst->push_back( curLoc ); - SG_LOG(SG_GENERAL, SG_DEBUG, "adding Linear Anchor node at (" << curLoc.x() << "," << curLoc.y() << ")"); +// // just add the one vertex - linear +// dst->push_back( curLoc ); +// SG_LOG(SG_GENERAL, SG_DEBUG, "adding Linear Anchor node at (" << curLoc.x() << "," << curLoc.y() << ")"); + + if (num_segs > 1) + { + for (int p=0; pGetLoc(), nextNode->GetLoc(), (1.0f/num_segs) * (p+1) ); + + // add the feature vertex + dst->push_back( curLoc ); + + if (p==0) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "adding Linear anchor node at (" << curLoc.x() << "," << curLoc.y() << ")"); + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, " add linear node at (" << curLoc.x() << "," << curLoc.y() << ")"); + } + + // now set set prev and cur locations for the next iteration + curLoc = nextLoc; + } + } + else + { + nextLoc = nextNode->GetLoc(); + + // just add the one vertex - dist is small + dst->push_back( curLoc ); + + SG_LOG(SG_GENERAL, SG_DEBUG, "adding Linear Anchor node at (" << curLoc.x() << "," << curLoc.y() << ")"); + + curLoc = nextLoc; + } } } } +#if 0 void ExpandPoint( Point3D *prev, Point3D *cur, Point3D *next, double expand_by, double *heading, double *offset ) { double offset_dir; @@ -352,7 +414,9 @@ void ExpandPoint( Point3D *prev, Point3D *cur, Point3D *next, double expand_by, SG_LOG(SG_GENERAL, SG_DEBUG, "heading is " << *heading << " distance is " << *offset ); } +#endif +#if 0 void ClosedPoly::ExpandContour( point_list& src, TGPolygon& dst, double dist ) { point_list expanded_boundary; @@ -443,6 +507,7 @@ void ClosedPoly::ExpandContour( point_list& src, TGPolygon& dst, double dist ) dst.add_contour( expanded_boundary, 9 ); } +#endif // finish the poly - convert to TGPolygon, and tesselate void ClosedPoly::Finish() @@ -465,7 +530,7 @@ void ClosedPoly::Finish() // and add it to the geometry pre_tess.add_contour( dst_contour, 0 ); - // The convert the hole contours + // Then convert the hole contours for (unsigned int i=0; ipush_back( sp ); SG_LOG(SG_GENERAL, SG_DEBUG, "clipped = " << clipped.contours()); #if 1 - *accum = tgPolygonUnionClipper( pre_tess, *accum ); -#else *accum = tgPolygonUnion( pre_tess, *accum ); +#else + *accum = tgPolygonUnionClipper( pre_tess, *accum ); #endif tp = TGTexParams( pre_tess.get_pt(0,0), 5.0, 5.0, texture_heading ); texparams->push_back( tp ); diff --git a/src/Airports/GenAirports850/linearfeature.cxx b/src/Airports/GenAirports850/linearfeature.cxx index fa987bb2..4f864c65 100644 --- a/src/Airports/GenAirports850/linearfeature.cxx +++ b/src/Airports/GenAirports850/linearfeature.cxx @@ -6,26 +6,33 @@ #include +// for debugging clipping errors +#include + #include "beznode.hxx" #include "linearfeature.hxx" #include "math.h" void LinearFeature::ConvertContour( BezContour* src, bool closed ) { - BezNode* prevNode; BezNode* curNode; BezNode* nextNode; - Point3D prevLoc; Point3D curLoc; Point3D nextLoc; Point3D cp1; Point3D cp2; int curve_type = CURVE_LINEAR; + double total_dist, linear_dist; + double meter_dist = 1.0f/96560.64f; + double theta1, theta2; + int num_segs = BEZIER_DETAIL; + Marking* cur_mark = NULL; Lighting* cur_light = NULL; + SG_LOG(SG_GENERAL, SG_DEBUG, " LinearFeature::ConvertContour - Creating a contour with " << src->size() << " nodes"); // clear anything in the point list @@ -34,19 +41,6 @@ void LinearFeature::ConvertContour( BezContour* src, bool closed ) // iterate through each bezier node in the contour for (unsigned int i=0; i <= src->size()-1; i++) { - SG_LOG(SG_GENERAL, SG_DEBUG, " LinearFeature::ConvertContour: Handling Node " << i << "\n\n"); - - if (i == 0) - { - // set prev node to last in the contour, as all contours must be closed - prevNode = src->at( src->size()-1 ); - } - else - { - // otherwise, it's just the previous index - prevNode = src->at( i-1 ); - } - curNode = src->at(i); if (i < src->size() - 1) @@ -132,38 +126,6 @@ void LinearFeature::ConvertContour( BezContour* src, bool closed ) } //////////////////////////////////////////////////////////////////////////////////// - - // determine the type of curve from prev (just to get correct prev location) - // once we start drawing the curve from cur to next, we can just remember the prev loc - if (prevNode->HasNextCp()) - { - // curve from prev is cubic or quadratic - if(curNode->HasPrevCp()) - { - // curve from prev is cubic : calculate the last location on the curve - prevLoc = CalculateCubicLocation( prevNode->GetLoc(), prevNode->GetNextCp(), curNode->GetPrevCp(), curNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (BEZIER_DETAIL-1) ); - } - else - { - // curve from prev is quadratic : use prev node next cp - prevLoc = CalculateQuadraticLocation( prevNode->GetLoc(), prevNode->GetNextCp(), curNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (BEZIER_DETAIL-1) ); - } - } - else - { - // curve from prev is quadratic or linear - if( curNode->HasPrevCp() ) - { - // curve from prev is quadratic : calculate the last location on the curve - prevLoc = CalculateQuadraticLocation( prevNode->GetLoc(), curNode->GetPrevCp(), curNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (BEZIER_DETAIL-1) ); - } - else - { - // curve from prev is linear : just use prev node location - prevLoc = prevNode->GetLoc(); - } - } - // now determine how we will iterate from current node to next node if( curNode->HasNextCp() ) { @@ -174,12 +136,17 @@ void LinearFeature::ConvertContour( BezContour* src, bool closed ) curve_type = CURVE_CUBIC; cp1 = curNode->GetNextCp(); cp2 = nextNode->GetPrevCp(); + total_dist = CubicDistance( curNode->GetLoc(), cp1, cp2, nextNode->GetLoc() ); + theta1 = SGMiscd::rad2deg( CalculateTheta( curNode->GetLoc(), cp1, nextNode->GetLoc()) ); + theta2 = SGMiscd::rad2deg( CalculateTheta( curNode->GetLoc(), cp2, nextNode->GetLoc()) ); } else { // curve is quadratic using current nodes cp as the cp curve_type = CURVE_QUADRATIC; cp1 = curNode->GetNextCp(); + total_dist = QuadraticDistance( curNode->GetLoc(), cp1, nextNode->GetLoc() ); + theta1 = SGMiscd::rad2deg( CalculateTheta( curNode->GetLoc(), cp1, nextNode->GetLoc()) ); } } else @@ -190,28 +157,123 @@ void LinearFeature::ConvertContour( BezContour* src, bool closed ) // curve is quadratic using next nodes cp as the cp curve_type = CURVE_QUADRATIC; cp1 = nextNode->GetPrevCp(); + total_dist = QuadraticDistance( curNode->GetLoc(), cp1, nextNode->GetLoc() ); + theta1 = SGMiscd::rad2deg( CalculateTheta( curNode->GetLoc(), cp1, nextNode->GetLoc()) ); } else { // curve is linear curve_type = CURVE_LINEAR; + total_dist = LinearDistance( curNode->GetLoc(), nextNode->GetLoc() ); } } + // One more test - some people are using bezier curves to draw straight lines - this can cause a bit of havoc... + // Sometimes, the control point lies just beyond the final point. We try to make a 'hook' at the end, which makes some really bad polys + // Just convert the entire segment to linear + // this can be detected in quadratic curves (current issue in LFKJ) when the contol point lies within the line generated from point 1 to point 2 + // theat close to 180 at the control point to the cur node and next node + if ( curve_type == CURVE_QUADRATIC ) + { + if ( (abs(theta1 - 180.0) < 5.0 ) || (abs(theta1) < 5.0 ) || (isnan(theta1)) ) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: Quadtratic curve with cp in line : convert to linear: " << description << ": theta is " << theta1 ); + curve_type = CURVE_LINEAR; + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: Quadtratic curve withOUT cp in line : keep quadtratic: " << description << ": theta is " << theta1 ); + } + } + + if ( curve_type == CURVE_CUBIC ) + { + if ( (abs(theta1 - 180.0) < 5.0 ) || (abs(theta1) < 5.0 ) || (isnan(theta1)) ) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: Cubic curve with cp1 in line : " << description << ": theta is " << theta1 ); + + if ( (abs(theta2 - 180.0) < 5.0 ) || (abs(theta2) < 5.0 ) || (isnan(theta2)) ) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\n and cp2 in line : " << description << ": theta is " << theta2 << " CONVERTING TO LINEAR" ); + + curve_type = CURVE_LINEAR; + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\n BUT cp2 NOT in line : " << description << ": theta is " << theta2 ); + } + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: Cubic curve withOUT cp1 in line : keep quadtratic: " << description << ": theta is " << theta1 ); + + if ( (abs(theta2 - 180.0) < 5.0 ) || (abs(theta2) < 5.0 ) || (isnan(theta2)) ) + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\n BUT cp2 IS in line : " << description << ": theta is " << theta2 ); + } + else + { + SG_LOG(SG_GENERAL, SG_DEBUG, "\n AND cp2 NOT in line : " << description << ": theta is " << theta2 ); + } + } + } + + double num_meters = total_dist / meter_dist; + if (num_meters < 4.0f) + { + if (curve_type != CURVE_LINEAR) + { + // If total distance is < 4 meters, then we need to modify num Segments so that each segment >= 1/2 meter + num_segs = ((int)num_meters + 1) * 2; + SG_LOG(SG_GENERAL, SG_DEBUG, "Segment from (" << curNode->GetLoc().x() << "," << curNode->GetLoc().y() << ") to (" << nextNode->GetLoc().x() << "," << nextNode->GetLoc().y() << ")" ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Distance is " << num_meters << " ( < 4.0) so num_segs is " << num_segs ); + } + else + { + num_segs = 1; + } + } + else if (num_meters > 800.0f) + { + // If total distance is > 800 meters, then we need to modify num Segments so that each segment <= 100 meters + num_segs = num_meters / 100.0f + 1; + SG_LOG(SG_GENERAL, SG_DEBUG, "Segment from (" << curNode->GetLoc().x() << "," << curNode->GetLoc().y() << ") to (" << nextNode->GetLoc().x() << "," << nextNode->GetLoc().y() << ")" ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Distance is " << num_meters << " ( > 100.0) so num_segs is " << num_segs ); + } + else + { + if (curve_type != CURVE_LINEAR) + { + num_segs = 8; + SG_LOG(SG_GENERAL, SG_DEBUG, "Segment from (" << curNode->GetLoc().x() << "," << curNode->GetLoc().y() << ") to (" << nextNode->GetLoc().x() << "," << nextNode->GetLoc().y() << ")" ); + SG_LOG(SG_GENERAL, SG_DEBUG, " Distance is " << num_meters << " (OK) so num_segs is " << num_segs ); + } + else + { + num_segs = 1; + } + } + + // if only one segment, revert to linear + if (num_segs == 1) + { + curve_type = CURVE_LINEAR; + } + // initialize current location curLoc = curNode->GetLoc(); if (curve_type != CURVE_LINEAR) { - for (int p=0; pGetLoc(), cp1, nextNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (p+1) ); + nextLoc = CalculateQuadraticLocation( curNode->GetLoc(), cp1, nextNode->GetLoc(), (1.0f/num_segs) * (p+1) ); } else { - nextLoc = CalculateCubicLocation( curNode->GetLoc(), cp1, cp2, nextNode->GetLoc(), (1.0f/BEZIER_DETAIL) * (p+1) ); + nextLoc = CalculateCubicLocation( curNode->GetLoc(), cp1, cp2, nextNode->GetLoc(), (1.0f/num_segs) * (p+1) ); } // add the feature vertex @@ -227,26 +289,16 @@ void LinearFeature::ConvertContour( BezContour* src, bool closed ) } // now set set prev and cur locations for the next iteration - prevLoc = curLoc; curLoc = nextLoc; } } else { - // For linear features, sometime long linear lines confuse the tesselator. Add intermediate nodes to keep the rectangles from - // getting too long. - double az1 = 0.0f; - double az2 = 0.0f; - double dist = 0.0f; - // calculate linear distance to determine how many segments we want Point3D destLoc = nextNode->GetLoc(); - geo_inverse_wgs_84( curLoc.y(), curLoc.x(), destLoc.y(), destLoc.x(), &az1, &az2, &dist); - if (dist > 100.0) + if (num_segs > 1) { - int num_segs = (dist / 100.0f) + 1; - for (int p=0; py(), 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 { SG_LOG(SG_GENERAL, SG_DEBUG, "\nLinearFeature: (theta NOT close to 180) " << description << ": theta is " << theta ); @@ -498,7 +554,9 @@ Point3D midpoint( Point3D p0, Point3D p1 ) return Point3D( (p0.x() + p1.x()) / 2, (p0.y() + p1.y()) / 2, (p0.z() + p1.z()) / 2 ); } -int LinearFeature::Finish( bool closed ) +#define DEBUG_LF (0) + +int LinearFeature::Finish( bool closed, unsigned int idx ) { TGPolygon poly; TGPolygon normals_poly; @@ -516,6 +574,15 @@ int LinearFeature::Finish( bool closed ) double light_delta = 0; double pt_x = 0, pt_y = 0; +#if DEBUG_LF + void* ds_id; + void* l_id; + + // Create a datasource for each linear feature + char ds_name[128]; + sprintf(ds_name, "./lf_debug/%04d_%s", idx, description.c_str()); + ds_id = tgShapefileOpenDatasource( ds_name ); +#endif // create the inner and outer boundaries to generate polys // this generates 2 point lists for the contours, and remembers @@ -644,6 +711,13 @@ int LinearFeature::Finish( bool closed ) exit(1); } +#if DEBUG_LF + // Create a new layer in the datasource for each Mark + char layer_name[128]; + sprintf( layer_name, "%04d_%s", i, material.c_str() ); + l_id = tgShapefileOpenLayer( ds_id, layer_name ); +#endif + last_end_v = 0.0f; for (unsigned int j = marks[i]->start_idx; j <= marks[i]->end_idx; j++) { @@ -681,6 +755,12 @@ int LinearFeature::Finish( bool closed ) poly.add_node( 0, cur_outer ); poly.add_node( 0, cur_inner ); +#if DEBUG_LF + char feature_name[128]; + sprintf( feature_name, "%04d", j); + tgShapefileCreateFeature( ds_id, l_id, poly, feature_name ); +#endif + sp.erase(); sp.set_poly( poly ); sp.set_material( material ); @@ -699,6 +779,11 @@ int LinearFeature::Finish( bool closed ) } } +#if DEBUG_LF + // Close the datasource + tgShapefileCloseDatasource( ds_id ); +#endif + // now generate the supoerpoly list for lights with constant distance between lights (depending on feature type) for (unsigned int i=0; ipush_back( marking_polys[i] ); +#if LF_DEBUG + if ( (debug) && ( i == 78 ) ) { + void* ds_id; + void* l_id; + + SG_LOG(SG_GENERAL, SG_INFO, "Problem poly: " << poly ); + + char ds_name[128]; + sprintf(ds_name, "./lf_debug/problem"); + ds_id = tgShapefileOpenDatasource( ds_name ); + + char layer_name[128]; + sprintf( layer_name, "problem"); + l_id = tgShapefileOpenLayer( ds_id, layer_name ); + + char feature_name[128]; + sprintf( feature_name, "prob"); + tgShapefileCreateFeature( ds_id, l_id, poly, feature_name ); + + sprintf( layer_name, "accum"); + l_id = tgShapefileOpenLayer( ds_id, layer_name ); + + sprintf( feature_name, "accum"); + tgShapefileCreateFeature( ds_id, l_id, *line_accum, feature_name ); + + tgShapefileCloseDatasource( ds_id ); + } +#endif + + SG_LOG(SG_GENERAL, SG_DEBUG, "LinearFeature::BuildBtg: union poly " << i << " of " << marking_polys.size() ); + #if 1 *line_accum = tgPolygonUnion( poly, *line_accum ); #else diff --git a/src/Airports/GenAirports850/linearfeature.hxx b/src/Airports/GenAirports850/linearfeature.hxx index 20b3ebb2..e5218b48 100644 --- a/src/Airports/GenAirports850/linearfeature.hxx +++ b/src/Airports/GenAirports850/linearfeature.hxx @@ -91,8 +91,8 @@ public: contour.push_back( b ); } - int Finish( bool closed ); - int BuildBtg( float alt_m, superpoly_list* line_polys, texparams_list* line_tps, ClipPolyType* line_accum, superpoly_list* lights ); + int Finish( bool closed, unsigned int idx ); + int BuildBtg( float alt_m, superpoly_list* line_polys, texparams_list* line_tps, ClipPolyType* line_accum, superpoly_list* lights, bool debug ); int BuildBtg( float alt_m, superpoly_list* line_polys, texparams_list* line_tps, Polygons* line_accum, superpoly_list* lights ); private: diff --git a/src/Airports/GenAirports850/main.cxx b/src/Airports/GenAirports850/main.cxx index e0b05376..50a56a64 100644 --- a/src/Airports/GenAirports850/main.cxx +++ b/src/Airports/GenAirports850/main.cxx @@ -25,6 +25,7 @@ #include #include +#include #include "beznode.hxx" #include "closedpoly.hxx" @@ -50,10 +51,10 @@ using namespace std; static void usage( int argc, char **argv ) { SG_LOG(SG_GENERAL, SG_ALERT, "Usage " << argv[0] << " --input= " - << "--work= [ --start-id=abcd ] [ --nudge=n ] " + << "--work= [ --start-id=abcd ] [ --restart-id=abcd ] [ --nudge=n ] " << "[--min-lon=] [--max-lon=] [--min-lat=] [--max-lat=] " << "[--clear-dem-path] [--dem-path=] [--max-slope=] " - << "[ --airport=abcd ] [--tile=] [--chunk=] [--verbose] [--help]"); + << "[ --airport=abcd ] [--tile=] [--chunk=] [--verbose] [--help]"); } @@ -90,6 +91,8 @@ static void help( int argc, char **argv, const string_list& elev_src ) { cout << "a valid airport code eg. --airport-id=KORD, or a starting airport can be specified using --start-id=abcd \n"; cout << "where once again abcd is a valid airport code. In this case, all airports in the file subsequent to the \n"; cout << "start-id are done. This is convienient when re-starting after a previous error. \n"; + cout << "If you want to restart with the airport after a problam icao, use --restart-id=abcd, as this works the same as\n"; + cout << " with the exception that the airport abcd is skipped \n"; cout << "\nAn input area may be specified by lat and lon extent using min and max lat and lon. \n"; cout << "Alternatively, you may specify a chunk (10 x 10 degrees) or tile (1 x 1 degree) using a string \n"; cout << "such as eg. w080n40, e000s27. \n"; @@ -134,7 +137,9 @@ int main(int argc, char **argv) string work_dir = ""; string input_file = ""; string start_id = ""; + string restart_id = ""; string airport_id = ""; + string last_apt_file = "./last_apt.txt"; int arg_pos; for (arg_pos = 1; arg_pos < argc; arg_pos++) @@ -156,10 +161,18 @@ int main(int argc, char **argv) { start_id = arg.substr(11); } + else if ( arg.find("--restart-id=") == 0 ) + { + restart_id = arg.substr(13); + } else if ( arg.find("--nudge=") == 0 ) { nudge = atoi( arg.substr(8).c_str() ); } + else if ( arg.find("--last_apt_file=") == 0 ) + { + last_apt_file = arg.substr(16); + } else if ( arg.find("--min-lon=") == 0 ) { min_lon = atof( arg.substr(10).c_str() ); @@ -270,6 +283,9 @@ int main(int argc, char **argv) string counter_file = airportareadir+"/poly_counter"; poly_index_init( counter_file ); + // Initialize shapefile support (for debugging) + tgShapefileInit(); + sg_gzifstream in( input_file ); if ( !in.is_open() ) { @@ -289,7 +305,7 @@ int main(int argc, char **argv) SG_LOG(SG_GENERAL, SG_INFO, "Finished Adding airport - now parse"); // and start the parser - parser->Parse(); + parser->Parse( last_apt_file ); } else if ( start_id != "" ) { @@ -302,7 +318,23 @@ int main(int argc, char **argv) parser->AddAirports( position, min_lat, min_lon, max_lat, max_lon ); // parse all the airports that were found - parser->Parse(); + parser->Parse( last_apt_file ); + } + else if ( restart_id != "" ) + { + SG_LOG(SG_GENERAL, SG_INFO, "move forward airport after " << restart_id ); + + // scroll forward in datafile + position = parser->FindAirport( restart_id ); + + // add all remaining airports within boundary + parser->AddAirports( position, min_lat, min_lon, max_lat, max_lon ); + + // but remove the restart id - it's broken + parser->RemoveAirport( restart_id ); + + // parse all the airports that were found + parser->Parse( last_apt_file ); } else { @@ -310,12 +342,13 @@ int main(int argc, char **argv) parser->AddAirports( 0, min_lat, min_lon, max_lat, max_lon ); // and parser them - parser->Parse(); + parser->Parse( last_apt_file ); } delete parser; SG_LOG(SG_GENERAL, SG_INFO, "Done"); + exit(0); return 0; } diff --git a/src/Airports/GenAirports850/parser.cxx b/src/Airports/GenAirports850/parser.cxx index 8841324f..850e801d 100644 --- a/src/Airports/GenAirports850/parser.cxx +++ b/src/Airports/GenAirports850/parser.cxx @@ -306,7 +306,18 @@ void Parser::AddAirports( long start_pos, float min_lat, float min_lon, float ma } } -void Parser::Parse() +void Parser::RemoveAirport( string icao ) +{ + for (unsigned int i = 0; i < airport_icaos.size(); i++ ) { + if (airport_icaos[i] == icao) { + parse_positions.erase(parse_positions.begin()+i); + airport_icaos.erase(airport_icaos.begin()+i); + break; + } + } +} + +void Parser::Parse( string last_apt_file ) { char tmp[2048]; struct timeval parse_start; @@ -335,6 +346,11 @@ void Parser::Parse() SG_LOG( SG_GENERAL, SG_ALERT, "Parsing airport " << airport_icaos[i] << " at " << parse_positions[i] << " start time " << ctime(&parse_start.tv_sec) ); in.seekg(parse_positions[i], ios::beg); + // save the airport we are working on + char command[256]; + sprintf( command, "echo %s > %s", airport_icaos[i].c_str(), last_apt_file.c_str() ); + system( command ); + while ( !in.eof() && (cur_state != STATE_DONE ) ) { in.getline(tmp, 2048); @@ -760,7 +776,7 @@ int Parser::ParseLine(char* line) } if (cur_airport) { - cur_feat->Finish( true ); + cur_feat->Finish( true, cur_airport->NumFeatures() ); cur_airport->AddFeature( cur_feat ); } cur_feat = NULL; @@ -794,7 +810,7 @@ int Parser::ParseLine(char* line) } if (cur_airport) { - cur_feat->Finish( false ); + cur_feat->Finish( false, cur_airport->NumFeatures() ); cur_airport->AddFeature( cur_feat ); } } diff --git a/src/Airports/GenAirports850/parser.hxx b/src/Airports/GenAirports850/parser.hxx index 833bd2b9..12658634 100644 --- a/src/Airports/GenAirports850/parser.hxx +++ b/src/Airports/GenAirports850/parser.hxx @@ -88,7 +88,8 @@ public: long FindAirport( string icao ); void AddAirport( string icao ); void AddAirports( long start_pos, float min_lat, float min_lon, float max_lat, float max_lon ); - void Parse( void ); + void RemoveAirport( string icao ); + void Parse( string last_apt_file ); private: bool IsAirportDefinition( char* line, string icao ); diff --git a/src/Airports/GenAirports850/runway.cxx b/src/Airports/GenAirports850/runway.cxx index 662248ac..672f5ad2 100644 --- a/src/Airports/GenAirports850/runway.cxx +++ b/src/Airports/GenAirports850/runway.cxx @@ -163,6 +163,8 @@ int Runway::BuildBtg( float alt_m, superpoly_list* rwy_polys, texparams_list* te if (apt_base) { + // If we have shoulders, we need to grow even further... + // generate area around runways base = gen_runway_area_w_extend( 0.0, 20.0, -rwy.overrun[0], -rwy.overrun[1], 20.0); diff --git a/src/Airports/GenAirports850/runway.hxx b/src/Airports/GenAirports850/runway.hxx index 7f192072..e86c6951 100644 --- a/src/Airports/GenAirports850/runway.hxx +++ b/src/Airports/GenAirports850/runway.hxx @@ -50,7 +50,9 @@ public: void BuildShoulder( float alt_m, superpoly_list *rwy_polys, texparams_list *texparams, - ClipPolyType *accum ); + ClipPolyType *accum, + TGPolygon* apt_base, + TGPolygon* apt_clearing ); private: struct TGRunway { @@ -152,6 +154,7 @@ private: superpoly_list gen_ssalx( float alt_m, const string& kind, bool recip ); superpoly_list gen_malsx( float alt_m, const string& kind, bool recip ); }; + typedef std::vector RunwayList; diff --git a/src/Airports/GenAirports850/rwy_gen.cxx b/src/Airports/GenAirports850/rwy_gen.cxx index 10eb4c8b..ed2dbae8 100644 --- a/src/Airports/GenAirports850/rwy_gen.cxx +++ b/src/Airports/GenAirports850/rwy_gen.cxx @@ -23,6 +23,8 @@ #include #include #include + +#include "beznode.hxx" #include "runway.hxx" #include @@ -153,80 +155,78 @@ void Runway::gen_rwy( double alt_m, TGPolygon runway = gen_runway_w_mid( alt_m, 0, 0 ); TGPolygon runway_half; + for ( int rwhalf=0; rwhalf<2; ++rwhalf ){ -for ( int rwhalf=0; rwhalf<2; ++rwhalf ){ + if (rwhalf == 0) { - if (rwhalf == 0) { + //Create the first half of the runway (first entry in apt.dat) + runway_half.erase(); + runway_half.add_node( 0, runway.get_pt(0, 3) ); + runway_half.add_node( 0, runway.get_pt(0, 4) ); + runway_half.add_node( 0, runway.get_pt(0, 5) ); + runway_half.add_node( 0, runway.get_pt(0, 2) ); + } + + else if (rwhalf == 1) { - //Create the first half of the runway (first entry in apt.dat) - runway_half.erase(); - runway_half.add_node( 0, runway.get_pt(0, 3) ); - runway_half.add_node( 0, runway.get_pt(0, 4) ); - runway_half.add_node( 0, runway.get_pt(0, 5) ); - runway_half.add_node( 0, runway.get_pt(0, 2) ); - } + //Create the second runway half from apt.dat + runway_half.erase(); + runway_half.add_node( 0, runway.get_pt(0, 0) ); + runway_half.add_node( 0, runway.get_pt(0, 1) ); + runway_half.add_node( 0, runway.get_pt(0, 2) ); + runway_half.add_node( 0, runway.get_pt(0, 5) ); + } - else if (rwhalf == 1) { + Point3D p; + SG_LOG(SG_GENERAL, SG_DEBUG, "raw runway half pts (run " << rwhalf << ")"); + for ( i = 0; i < runway_half.contour_size( 0 ); ++i ) { + p = runway_half.get_pt(0, i); + SG_LOG(SG_GENERAL, SG_DEBUG, " point = " << p); + } - //Create the second runway half from apt.dat - runway_half.erase(); - runway_half.add_node( 0, runway.get_pt(0, 0) ); - runway_half.add_node( 0, runway.get_pt(0, 1) ); - runway_half.add_node( 0, runway.get_pt(0, 2) ); - runway_half.add_node( 0, runway.get_pt(0, 5) ); - } - - Point3D p; - SG_LOG(SG_GENERAL, SG_DEBUG, "raw runway half pts (run " << rwhalf << ")"); - for ( i = 0; i < runway_half.contour_size( 0 ); ++i ) { - p = runway_half.get_pt(0, i); - SG_LOG(SG_GENERAL, SG_DEBUG, " point = " << p); - } - - double length = rwy.length / 2.0; - if ( length < 3075 * SG_FEET_TO_METER ) { - SG_LOG( SG_GENERAL, SG_DEBUG, - "Runway " << rwy.rwnum[0] << " is not long enough (" + double length = rwy.length / 2.0; + if ( length < 3075 * SG_FEET_TO_METER ) { + SG_LOG( SG_GENERAL, SG_DEBUG, + "Runway " << rwy.rwnum[0] << " is not long enough (" << rwy.length << ") for precision markings!"); - } + } - double start1_pct = 0.0; - double end1_pct = 0.0; - double heading = 0.0; - string rwname; + double start1_pct = 0.0; + double end1_pct = 0.0; + double heading = 0.0; + string rwname; + // + // Displaced threshold if it exists + // - // - // Displaced threshold if it exists - // - - if (rwhalf == 0) { - heading = rwy.heading + 180.0; + if (rwhalf == 0) { + heading = rwy.heading + 180.0; rwname = rwy.rwnum[0]; - } - else if (rwhalf == 1) { + } + else if (rwhalf == 1) { heading = rwy.heading; rwname = rwy.rwnum[1]; - } - SG_LOG( SG_GENERAL, SG_DEBUG, "runway marking = " << rwy.marking[rwhalf] ); - if ( rwy.threshold[rwhalf] > 0.0 ) { - SG_LOG( SG_GENERAL, SG_DEBUG, "Displaced threshold for RW side " << rwhalf << " is " + } + SG_LOG( SG_GENERAL, SG_DEBUG, "runway marking = " << rwy.marking[rwhalf] ); + if ( rwy.threshold[rwhalf] > 0.0 ) { + SG_LOG( SG_GENERAL, SG_DEBUG, "Displaced threshold for RW side " << rwhalf << " is " << rwy.threshold[rwhalf] ); - // reserve 90' for final arrows - double thresh = rwy.threshold[rwhalf] - 90.0 * SG_FEET_TO_METER; + // reserve 90' for final arrows + double thresh = rwy.threshold[rwhalf] - 90.0 * SG_FEET_TO_METER; - // number of full center arrows - int count = (int)(thresh / 200.0 * SG_FEET_TO_METER); + // number of full center arrows + int count = (int)(thresh / 200.0 * SG_FEET_TO_METER); - // length of starting partial arrow - double part_len = thresh - ( count * 200.0 * SG_FEET_TO_METER); - double tex_pct = (200.0 * SG_FEET_TO_METER - part_len) / 200.0 * SG_FEET_TO_METER; + // length of starting partial arrow + double part_len = thresh - ( count * 200.0 * SG_FEET_TO_METER); + double tex_pct = (200.0 * SG_FEET_TO_METER - part_len) / 200.0 * SG_FEET_TO_METER; - // starting (possibly partial chunk) - start1_pct = end1_pct; - end1_pct = start1_pct + ( part_len / length ); - gen_runway_section( runway_half, + // starting (possibly partial chunk) + start1_pct = end1_pct; + end1_pct = start1_pct + ( part_len / length ); + gen_runway_section( runway_half, start1_pct, end1_pct, 0.0, 1.0, 0.0, 1.0, tex_pct, 1.0, @@ -234,72 +234,70 @@ for ( int rwhalf=0; rwhalf<2; ++rwhalf ){ material, "dspl_thresh", rwy_polys, texparams, accum ); - // main chunks - for ( i = 0; i < count; ++i ) { - start1_pct = end1_pct; - end1_pct = start1_pct + ( 200.0 * SG_FEET_TO_METER / length ); - gen_runway_section( runway_half, + // main chunks + for ( i = 0; i < count; ++i ) { + start1_pct = end1_pct; + end1_pct = start1_pct + ( 200.0 * SG_FEET_TO_METER / length ); + gen_runway_section( runway_half, start1_pct, end1_pct, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, heading, material, "dspl_thresh", rwy_polys, texparams, accum ); - } + } - // final arrows - start1_pct = end1_pct; - end1_pct = start1_pct + ( 90.0 * SG_FEET_TO_METER / length ); - gen_runway_section( runway_half, + // final arrows + start1_pct = end1_pct; + end1_pct = start1_pct + ( 90.0 * SG_FEET_TO_METER / length ); + gen_runway_section( runway_half, start1_pct, end1_pct, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, heading, material, "dspl_arrows", rwy_polys, texparams, accum ); + } + + + if (rwy.marking[rwhalf] == 0){ + + // No threshold + + start1_pct = end1_pct; + end1_pct = start1_pct + ( 10 / length ); + gen_runway_section( runway_half, + start1_pct, end1_pct, + 0.0, 1.0, + 0.0, 1.0, 0.0, 1.0, + heading, + material, "no_threshold", + rwy_polys, texparams, accum ); + } else { + + // Thresholds for all others + + start1_pct = end1_pct; + end1_pct = start1_pct + ( 202.0 * SG_FEET_TO_METER / length ); + gen_runway_section( runway_half, + start1_pct, end1_pct, + 0.0, 1.0, + 0.0, 1.0, 0.0, 1.0, + heading, + material, "threshold", + rwy_polys, texparams, accum ); } - - if (rwy.marking[rwhalf] == 0){ - - // No threshold - - start1_pct = end1_pct; - end1_pct = start1_pct + ( 10 / length ); - gen_runway_section( runway_half, - start1_pct, end1_pct, - 0.0, 1.0, - 0.0, 1.0, 0.0, 1.0, - heading, - material, "no_threshold", - rwy_polys, texparams, accum ); - } else { - - // Thresholds for all others - - start1_pct = end1_pct; - end1_pct = start1_pct + ( 202.0 * SG_FEET_TO_METER / length ); - gen_runway_section( runway_half, - start1_pct, end1_pct, - 0.0, 1.0, - 0.0, 1.0, 0.0, 1.0, - heading, - material, "threshold", - rwy_polys, texparams, accum ); - } - // Runway designation block gen_rw_designation( material, runway_half, heading, rwname, start1_pct, end1_pct, rwy_polys, texparams, accum ); - - if (rwy.marking[rwhalf] > 1){ - // Generate remaining markings depending on type of runway - gen_rw_marking( runway_half, - start1_pct, end1_pct, - heading, material, - rwy_polys, texparams, accum, rwy.marking[rwhalf] ); + // Generate remaining markings depending on type of runway + gen_rw_marking( runway_half, + start1_pct, end1_pct, + heading, material, + rwy_polys, texparams, accum, rwy.marking[rwhalf] ); } // @@ -314,48 +312,49 @@ for ( int rwhalf=0; rwhalf<2; ++rwhalf ){ double rest1_inc = (1.0 - end1_pct) / divs; while ( end1_pct < 1.0 ) { - start1_pct = end1_pct; - end1_pct = start1_pct + rest1_inc; + start1_pct = end1_pct; + end1_pct = start1_pct + rest1_inc; - gen_runway_section( runway_half, + gen_runway_section( runway_half, start1_pct, end1_pct, - 0.0, 1.0, - 0.0, 1.0, 0.0, 1.0, - heading, - material, "rest", - rwy_polys, texparams, accum ); - } + 0.0, 1.0, + 0.0, 1.0, 0.0, 1.0, + heading, + material, "rest", + rwy_polys, texparams, accum ); + } - gen_runway_overrun( runway_half, rwhalf, + gen_runway_overrun( runway_half, rwhalf, material, rwy_polys, texparams, accum ); - - -} + } } void Runway::BuildShoulder( float alt_m, superpoly_list *rwy_polys, texparams_list *texparams, - ClipPolyType *accum ) + ClipPolyType *accum, + TGPolygon* apt_base, + TGPolygon* apt_clearing ) { + TGPolygon base, safe_base; + string shoulder_surface = ""; double shoulder_width = 0.0f; - if (rwy.shoulder > 0){ // Add a shoulder to the runway - //shoulder_width = rwy.width * 0.15; - //if (shoulder_width > 11.0){ - // shoulder_width = 11.0; - //} + if (rwy.shoulder > 0){ + // Add a shoulder to the runway shoulder_width = 11.0f; if (rwy.shoulder == 1){ shoulder_surface = "pa_shoulder"; } else if (rwy.shoulder == 2){ shoulder_surface = "pc_shoulder"; - } else SG_LOG(SG_GENERAL, SG_ALERT, "Unknown shoulder surface code = " << rwy.shoulder ); - - } else if (rwy.shoulder == 0){ // We add a fake shoulder if the runway has an asphalt or concrete surface + } else { + SG_LOG(SG_GENERAL, SG_ALERT, "Unknown shoulder surface code = " << rwy.shoulder ); + } + } else { + // We add a fake shoulder if the runway has an asphalt or concrete surface shoulder_width = 1.0; if (rwy.surface == 1){ shoulder_surface = "pa_shoulder_f"; @@ -363,48 +362,104 @@ void Runway::BuildShoulder( float alt_m, shoulder_surface = "pc_shoulder_f"; } } + SG_LOG(SG_GENERAL, SG_DEBUG, "Shoulder width = " << shoulder_width ); SG_LOG(SG_GENERAL, SG_DEBUG, "Shoulder surface is: " << shoulder_surface ); if (shoulder_width > 0.0f) { - // we need to break these shoulders up into managable pieces, as we're asked to triangulate // 3-4 km long by 1m wide strips - jrs-can't handle it. double max_dist = (double)shoulder_width * 25.0f; - int numSegs = (rwy.length / max_dist) + 1; - double dist = rwy.length / (double)numSegs; + if (max_dist > 100.0f) { + max_dist = 100.0f; + } + int numSegs = (rwy.length / max_dist) + 1; + double dist = rwy.length / numSegs; - // Create both shoulder sides - for (int i=0; i<2; ++i){ - double step; - double lat = 0.0f, lon = 0.0f, r; + TGPolygon poly; + TGSuperPoly sp; + TGTexParams tp; - /* nudge the shoulders so the really long lines overlap the runway a bit */ - /* If the are 'equal' there's a good chance roundoff error can create a */ - /* REALY thin long polygon, which causes a segfault */ - if (i == 0){ - step= (rwy.width + shoulder_width) * 0.5; - } else if (i == 1) { - step= -(rwy.width + shoulder_width) * 0.5; + double lat = 0.0f; + double lon = 0.0f; + double r = 0.0f; + Point3D inner_start, inner_end; + Point3D outer_start, outer_end; + Point3D curInnerLoc, nextInnerLoc; + Point3D curOuterLoc, nextOuterLoc; + + // Create two paralell lines from start position to end position, and interpolate in between + // many airports line the taxiway directly to the corner of the runway. This can create problems, + // so extend the shoulders 0.5 meters past each end of the runway + for (int i=0; i<2; i++) { + double rev_hdg = rwy.heading - 180.0; + if ( rev_hdg < 0 ) { rev_hdg += 360.0; } + + if (i == 0) { + // left side + double left_hdg = rwy.heading - 90.0; + if ( left_hdg < 0 ) { left_hdg += 360.0; } + + geo_direct_wgs_84 ( 0, rwy.lat[0], rwy.lon[0], left_hdg, rwy.width*.5, &lat, &lon, &r ); + geo_direct_wgs_84 ( 0, lat, lon, rev_hdg, 0.5f, &lat, &lon, &r ); + + inner_start = Point3D( lon, lat, 0.0f ); + geo_direct_wgs_84 ( 0, lat, lon, rwy.heading, rwy.length+1.0f, &lat, &lon, &r ); + inner_end = Point3D( lon, lat, 0.0f ); + + geo_direct_wgs_84 ( 0, rwy.lat[0], rwy.lon[0], left_hdg, rwy.width*.5 + shoulder_width, &lat, &lon, &r ); + geo_direct_wgs_84 ( 0, lat, lon, rev_hdg, 0.5f, &lat, &lon, &r ); + + outer_start = Point3D( lon, lat, 0.0f ); + geo_direct_wgs_84 ( 0, lat, lon, rwy.heading, rwy.length+1.0f, &lat, &lon, &r ); + outer_end = Point3D( lon, lat, 0.0f ); + } else { + // right side + double right_hdg = rwy.heading + 90.0; + if ( right_hdg > 360 ) { right_hdg -= 360.0; } + + geo_direct_wgs_84 ( 0, rwy.lat[0], rwy.lon[0], right_hdg, rwy.width*.5, &lat, &lon, &r ); + geo_direct_wgs_84 ( 0, lat, lon, rev_hdg, 0.5f, &lat, &lon, &r ); + + inner_start = Point3D( lon, lat, 0.0f ); + geo_direct_wgs_84 ( 0, lat, lon, rwy.heading, rwy.length+1.0f, &lat, &lon, &r ); + inner_end = Point3D( lon, lat, 0.0f ); + + geo_direct_wgs_84 ( 0, rwy.lat[0], rwy.lon[0], right_hdg, rwy.width*.5 + shoulder_width, &lat, &lon, &r ); + geo_direct_wgs_84 ( 0, lat, lon, rev_hdg, 0.5f, &lat, &lon, &r ); + + outer_start = Point3D( lon, lat, 0.0f ); + geo_direct_wgs_84 ( 0, lat, lon, rwy.heading, rwy.length+1.0f, &lat, &lon, &r ); + outer_end = Point3D( lon, lat, 0.0f ); } - double left_hdg = rwy.heading - 90.0; - if ( left_hdg < 0 ) { left_hdg += 360.0; } + curInnerLoc = inner_start; + curOuterLoc = outer_start; - geo_direct_wgs_84 ( alt_m, rwy.lat[0], rwy.lon[0], left_hdg, step, &lat, &lon, &r ); - Point3D ref = Point3D( lon, lat, 0.0f ); - - for (int j=0; jpush_back( sp ); #if 1 - *accum = tgPolygonUnion( shoulderSegment, *accum ); + *accum = tgPolygonUnion( poly, *accum ); #else - *accum = tgPolygonUnionClipper( shoulderSegment, *accum ); + *accum = tgPolygonUnionClipper( poly, *accum ); #endif - tp = TGTexParams( shoulderSegment.get_pt(0,0), -shoulder_width, dist, rwy.heading ); - if (i == 0){ + tp = TGTexParams( poly.get_pt(0,0), shoulder_width, dist, rwy.heading ); + tp.set_maxv(dist); + // reverse u direction for right side + if ( i == 1 ) { tp.set_maxu(0); tp.set_minu(1); } texparams->push_back( tp ); + + // Add to base / safe base + base = tgPolygonExpand( poly, 20.0f ); + safe_base = tgPolygonExpand( poly, 50.0f ); + + // add this to the airport clearing + *apt_clearing = tgPolygonUnion(safe_base, *apt_clearing); + + // and add the clearing to the base + *apt_base = tgPolygonUnion( base, *apt_base ); + + // now set cur locations for the next iteration + curInnerLoc = nextInnerLoc; + curOuterLoc = nextOuterLoc; } - } + } } } diff --git a/src/BuildTiles/Main/CMakeLists.txt b/src/BuildTiles/Main/CMakeLists.txt index c903c0ef..7eba5e69 100644 --- a/src/BuildTiles/Main/CMakeLists.txt +++ b/src/BuildTiles/Main/CMakeLists.txt @@ -14,6 +14,7 @@ target_link_libraries(fgfs-construct Polygon Geometry Array Optimize Output landcover poly2tri TriangleJRS + ${GDAL_LIBRARY} ${SIMGEAR_CORE_LIBRARIES} ${SIMGEAR_CORE_LIBRARY_DEPENDENCIES} ${GPC_LIBRARY}) diff --git a/src/Lib/Geometry/CMakeLists.txt b/src/Lib/Geometry/CMakeLists.txt index 97019880..05fbbdd5 100644 --- a/src/Lib/Geometry/CMakeLists.txt +++ b/src/Lib/Geometry/CMakeLists.txt @@ -1,4 +1,6 @@ - +if (GDAL_FOUND) +include_directories(${GDAL_INCLUDE_DIR}) +endif (GDAL_FOUND) add_library(Geometry STATIC contour_tree.cxx diff --git a/src/Lib/Geometry/poly_support.cxx b/src/Lib/Geometry/poly_support.cxx index faac9b92..6ff8acc7 100644 --- a/src/Lib/Geometry/poly_support.cxx +++ b/src/Lib/Geometry/poly_support.cxx @@ -1452,25 +1452,161 @@ TGPolygon remove_bad_contours( const TGPolygon &poly ) { // remove any too small contours TGPolygon remove_tiny_contours( const TGPolygon &poly ) { + double min_area = SG_EPSILON*SG_EPSILON; TGPolygon result; result.erase(); for ( int i = 0; i < poly.contours(); ++i ) { - point_list contour = poly.get_contour( i ); + point_list contour = poly.get_contour( i ); - double area=poly.area_contour(i); - - if (-SG_EPSILON*SG_EPSILON +const char* format_name="ESRI Shapefile"; + +void tgShapefileInit( void ) +{ + OGRRegisterAll(); +} + +void* tgShapefileOpenDatasource( const char* datasource_name ) +{ + OGRDataSource *datasource; + OGRSFDriver *ogrdriver; + + ogrdriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(format_name); + if (!ogrdriver) { + SG_LOG(SG_GENERAL, SG_ALERT, "Unknown datasource format driver: " << format_name); + exit(1); + } + + datasource = ogrdriver->CreateDataSource(datasource_name, NULL); + if (!datasource) { + SG_LOG(SG_GENERAL, SG_ALERT, "Unable to create datasource: " << datasource_name); + exit(1); + } + + return (void*)datasource; +} + +void* tgShapefileOpenLayer( void* ds_id, const char* layer_name ) { + OGRDataSource* datasource = (OGRDataSource *)ds_id; + OGRLayer* layer; + + OGRSpatialReference srs; + srs.SetWellKnownGeogCS("WGS84"); + layer = datasource->CreateLayer( layer_name, &srs, wkbPolygon25D, NULL); + if (!layer) { + SG_LOG(SG_GENERAL, SG_ALERT, "Creation of layer '" << layer_name << "' failed"); + return NULL; + } + + OGRFieldDefn descriptionField("ID", OFTString); + descriptionField.SetWidth(128); + + if( layer->CreateField( &descriptionField ) != OGRERR_NONE ) { + SG_LOG(SG_GENERAL, SG_ALERT, "Creation of field 'Description' failed"); + } + + return (void*)layer; +} + +void tgShapefileCreateFeature( void* ds_id, void* l_id, const TGPolygon &poly, const char* description ) +{ + OGRDataSource* datasource = (OGRDataSource *)ds_id; + OGRLayer* layer = (OGRLayer*)l_id; + OGRPolygon* polygon = new OGRPolygon(); + + for ( int i = 0; i < poly.contours(); i++ ) { + bool skip_ring=false; + point_list contour = poly.get_contour( i ); + + if (contour.size()<3) { + SG_LOG(SG_GENERAL, SG_DEBUG, "Polygon with less than 3 points"); + skip_ring=true; + } + + // FIXME: Current we ignore the hole-flag and instead assume + // that the first ring is not a hole and the rest + // are holes + OGRLinearRing *ring=new OGRLinearRing(); + for (unsigned int pt = 0; pt < contour.size(); pt++) { + OGRPoint *point=new OGRPoint(); + + point->setX( contour[pt].x() ); + point->setY( contour[pt].y() ); + point->setZ( 0.0 ); + ring->addPoint(point); + } + ring->closeRings(); + + if (!skip_ring) { + polygon->addRingDirectly(ring); + } + + OGRFeature* feature = NULL; + feature = new OGRFeature( layer->GetLayerDefn() ); + feature->SetField("ID", description); + feature->SetGeometry(polygon); + if( layer->CreateFeature( feature ) != OGRERR_NONE ) + { + SG_LOG(SG_GENERAL, SG_ALERT, "Failed to create feature in shapefile"); + } + OGRFeature::DestroyFeature(feature); + } +} + +void tgShapefileCloseLayer( void* l_id ) +{ + OGRLayer* layer = (OGRLayer *)l_id; + + //OGRLayer::DestroyLayer( layer ); +} + +void tgShapefileCloseDatasource( void* ds_id ) +{ + OGRDataSource* datasource = (OGRDataSource *)ds_id; + + OGRDataSource::DestroyDataSource( datasource ); +} diff --git a/src/Lib/Geometry/poly_support.hxx b/src/Lib/Geometry/poly_support.hxx index 59da8a00..99f00966 100644 --- a/src/Lib/Geometry/poly_support.hxx +++ b/src/Lib/Geometry/poly_support.hxx @@ -107,6 +107,16 @@ TGPolygon remove_bad_contours( const TGPolygon &poly ); // remove any too small contours TGPolygon remove_tiny_contours( const TGPolygon &poly ); +TGPolygon remove_small_contours( const TGPolygon &poly ); + + +// Write Polygons to Shapefile Support +void tgShapefileInit( void ); +void* tgShapefileOpenDatasource( const char* datasource_name ); +void* tgShapefileOpenLayer( void* ds_id, const char* layer_name ); +void tgShapefileCreateFeature( void* ds_id, void* l_id, const TGPolygon &poly, const char* feature_name ); +void tgShapefileCloseLayer( void* l_id ); +void tgShapefileCloseDatasource( void* ds_id ); #endif // _POLY_SUPPORT_HXX