diff --git a/src/Airports/GenAirports850/CMakeLists.txt b/src/Airports/GenAirports850/CMakeLists.txt index d5154a52..fd02c58c 100644 --- a/src/Airports/GenAirports850/CMakeLists.txt +++ b/src/Airports/GenAirports850/CMakeLists.txt @@ -4,7 +4,6 @@ include_directories(${PROJECT_SOURCE_DIR}/src/Lib) add_executable(genapts850 airport.hxx airport.cxx apt_math.hxx apt_math.cxx - apt_surface.hxx apt_surface.cxx beznode.hxx closedpoly.hxx closedpoly.cxx debug.hxx debug.cxx diff --git a/src/Airports/GenAirports850/airport.cxx b/src/Airports/GenAirports850/airport.cxx index 99eb999b..dcfaa1a2 100644 --- a/src/Airports/GenAirports850/airport.cxx +++ b/src/Airports/GenAirports850/airport.cxx @@ -213,7 +213,7 @@ bool Airport::isDebugFeature( int feat ) // TODO : Add somewhere // Determine node elevations of a point_list based on the provided // TGAptSurface. Offset is added to the final elevation -static std::vector calc_elevations( TGAptSurface &surf, const std::vector& geod_nodes, double offset ) +static std::vector calc_elevations( const tgSurface& surf, const std::vector& geod_nodes, double offset ) { std::vector result = geod_nodes; for ( unsigned int i = 0; i < result.size(); ++i ) { @@ -225,9 +225,7 @@ static std::vector calc_elevations( TGAptSurface &surf, const std::vecto } -static tgContour calc_elevations( TGAptSurface &surf, - const tgContour& geod_nodes, - double offset ) +static tgContour calc_elevations( const tgSurface& surf, const tgContour& geod_nodes, double offset ) { tgContour result = geod_nodes; for ( unsigned int i = 0; i < result.GetSize(); ++i ) { @@ -240,9 +238,7 @@ static tgContour calc_elevations( TGAptSurface &surf, return result; } -static double calc_elevation( TGAptSurface &surf, - const SGGeod& node, - double offset ) +static double calc_elevation( const tgSurface& surf, const SGGeod& node, double offset ) { double elev = surf.query( node ); elev += offset; @@ -253,11 +249,10 @@ static double calc_elevation( TGAptSurface &surf, // Determine node elevations of each node of a TGPolygon based on the // provided TGAptSurface. Offset is added to the final elevation -static tgPolygon calc_elevations( TGAptSurface &surf, - const tgPolygon& poly, - double offset ) +static tgPolygon calc_elevations( const tgSurface& surf, const tgPolygon& poly, double offset ) { tgPolygon result; + for ( unsigned int i = 0; i < poly.Contours(); ++i ) { tgContour contour = poly.GetContour( i ); tgContour elevated = calc_elevations( surf, contour, offset ); @@ -1057,9 +1052,9 @@ void Airport::BuildBtg(const std::string& root, const string_list& elev_src ) GENAPT_LOG(SG_GENERAL, SG_DEBUG, " max: " << max_deg ); GENAPT_LOG(SG_GENERAL, SG_DEBUG, " average: " << average ); - tg::Rectangle aptBounds(min_deg, max_deg); - - TGAptSurface apt_surf( root, elev_src, aptBounds, average ); + // TODO elevation queries should be performed as member functions of surface + tgRectangle aptBounds(min_deg, max_deg); + tgSurface apt_surf( root, elev_src, aptBounds, average, slope_max, slope_eps ); GENAPT_LOG(SG_GENERAL, SG_DEBUG, "Airport surface created"); // add light points diff --git a/src/Airports/GenAirports850/debug.hxx b/src/Airports/GenAirports850/debug.hxx index cae65314..277be825 100644 --- a/src/Airports/GenAirports850/debug.hxx +++ b/src/Airports/GenAirports850/debug.hxx @@ -4,6 +4,7 @@ #include #include +#include #ifndef __DEBUG_HXX__ #define __DEBUG_HXX__ diff --git a/src/Airports/GenAirports850/elevations.cxx b/src/Airports/GenAirports850/elevations.cxx index cbc6f281..63d5511a 100644 --- a/src/Airports/GenAirports850/elevations.cxx +++ b/src/Airports/GenAirports850/elevations.cxx @@ -31,7 +31,7 @@ #include #include "global.hxx" -#include "apt_surface.hxx" +#include "debug.hxx" // lookup node elevations for each point in the SGGeod list. Returns @@ -126,138 +126,4 @@ double tgAverageElevation( const std::string &root, const string_list elev_src, GENAPT_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of point list = " << average); return average; -} - - -// lookup node elevations for each point in the specified simple -// matrix. Returns average of all points. -void tgCalcElevations( const std::string &root, const string_list elev_src, - SimpleMatrix &Pts, const double average ) -{ - bool done = false; - int i, j; - TGArray array; - - // just bail if no work to do - if ( Pts.rows() == 0 || Pts.cols() == 0 ) { - return; - } - - // set all elevations to -9999 - for ( j = 0; j < Pts.rows(); ++j ) { - for ( i = 0; i < Pts.cols(); ++i ) { - SGGeod p = Pts.element(i, j); - p.setElevationM(-9999.0); - Pts.set(i, j, p); - } - } - - while ( !done ) { - // find first node with -9999 elevation - SGGeod first = SGGeod(); - bool found_one = false; - for ( j = 0; j < Pts.rows(); ++j ) { - for ( i = 0; i < Pts.cols(); ++i ) { - SGGeod p = Pts.element(i,j); - if ( p.getElevationM() < -9000.0 && !found_one ) { - first = p; - found_one = true; - } - } - } - - if ( found_one ) { - SGBucket b( first ); - std::string base = b.gen_base_path(); - - // try the various elevation sources - j = 0; - bool found_file = false; - while ( !found_file && j < (int)elev_src.size() ) { - std::string array_path = root + "/" + elev_src[j] + "/" + base + "/" + b.gen_index_str(); - - if ( array.open(array_path) ) { - found_file = true; - GENAPT_LOG( SG_GENERAL, SG_DEBUG, "Using array_path = " << array_path ); - } - j++; - } - - // this will fill in a zero structure if no array data - // found/opened - array.parse( b ); - - // this will do a hasty job of removing voids by inserting - // data from the nearest neighbor (sort of) - array.remove_voids(); - - // update all the non-updated elevations that are inside - // this array file - double elev; - done = true; - for ( j = 0; j < Pts.rows(); ++j ) { - for ( i = 0; i < Pts.cols(); ++i ) { - SGGeod p = Pts.element(i,j); - if ( p.getElevationM() < -9000.0 ) { - done = false; - elev = array.altitude_from_grid( p.getLongitudeDeg() * 3600.0, - p.getLatitudeDeg() * 3600.0 ); - if ( elev > -9000 ) { - p.setElevationM( elev ); - Pts.set(i, j, p); - } - } - } - } - - array.close(); - - } else { - done = true; - } - } - -#ifdef DEBUG - // do some post processing for sanity's sake - // find the average height of the queried points - double total = 0.0; - int count = 0; - for ( j = 0; j < Pts.rows(); ++j ) { - for ( i = 0; i < Pts.cols(); ++i ) { - SGGeod p = Pts.element(i,j); - total += p.getElevationM(); - count++; - } - } - double grid_average = total / (double) count; - GENAPT_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of matrix = " << grid_average); -#endif -} - -// clamp all elevations to the specified range -void tgClampElevations( SimpleMatrix &Pts, - double center_m, - double max_clamp_m ) -{ - int i, j; - - // go through the elevations and clamp all elevations to within - // +/-max_m of the center_m elevation. - for ( j = 0; j < Pts.rows(); ++j ) { - for ( i = 0; i < Pts.cols(); ++i ) { - SGGeod p = Pts.element(i,j); - if ( p.getElevationM() < center_m - max_clamp_m ) { - GENAPT_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.getElevationM() - << " to " << center_m - max_clamp_m ); - p.setElevationM( center_m - max_clamp_m ); - Pts.set(i, j, p); - } - if ( p.getElevationM() > center_m + max_clamp_m ) { - GENAPT_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.getElevationM() - << " to " << center_m + max_clamp_m ); - p.setElevationM( center_m + max_clamp_m ); - Pts.set(i, j, p); - } - } - } -} +} \ No newline at end of file diff --git a/src/Airports/GenAirports850/elevations.hxx b/src/Airports/GenAirports850/elevations.hxx index c0290022..e992a218 100644 --- a/src/Airports/GenAirports850/elevations.hxx +++ b/src/Airports/GenAirports850/elevations.hxx @@ -20,8 +20,7 @@ // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. // -#include "apt_surface.hxx" - +#include // lookup node elevations for each point in the SGGeod list. Returns // average of all points. Doesn't modify the original list. @@ -30,7 +29,7 @@ double tgAverageElevation( const std::string &root, const string_list elev_src, // lookup node elevations for each point in the specified nurbs++ // matrix. -void tgCalcElevations( const std::string &root, const string_list elev_src, SimpleMatrix &Pts, double average ); +void tgCalcElevations( const std::string &root, const string_list elev_src, tgMatrix& Pts, double average ); // clamp all elevations to the specified range -void tgClampElevations( SimpleMatrix &Pts, double center_m, double max_clamp_m ); +void tgClampElevations( tgMatrix& Pts, double center_m, double max_clamp_m ); diff --git a/src/Airports/GenAirports850/global.hxx b/src/Airports/GenAirports850/global.hxx index ddcd615f..e0e1c496 100644 --- a/src/Airports/GenAirports850/global.hxx +++ b/src/Airports/GenAirports850/global.hxx @@ -31,15 +31,7 @@ extern int nudge; // Each polygon vertex is snapped to a grid with this resolution (~1cm by default) extern double gSnap; -// Final grid size for airport surface (in meters) -const double coarse_grid = 300.0; - -// compared to the average surface elevation, clamp all values within -// this many meters of the average -const double max_clamp = 100.0; - -// maximum slope (rise/run) allowed on an airport surface -extern double slope_max; // = 0.02; -const double slope_eps = 0.00001; +extern double slope_max; +extern double slope_eps; #endif diff --git a/src/Airports/GenAirports850/main.cxx b/src/Airports/GenAirports850/main.cxx index ad9bc3b8..75103312 100644 --- a/src/Airports/GenAirports850/main.cxx +++ b/src/Airports/GenAirports850/main.cxx @@ -98,8 +98,9 @@ static void help( int argc, char **argv, const string_list& elev_src ) { // TODO: where do these belong int nudge = 10; -double slope_max = 0.02; double gSnap = 0.00000001; // approx 1 mm +double slope_max = 0.02; +double slope_eps = 0.00001; int main(int argc, char **argv) { @@ -284,7 +285,7 @@ int main(int argc, char **argv) std::string lastaptfile = work_dir+"/last_apt"; - tg::Rectangle boundingBox(min, max); + tgRectangle boundingBox(min, max); boundingBox.sanify(); if ( work_dir == "" ) diff --git a/src/Airports/GenAirports850/scheduler.cxx b/src/Airports/GenAirports850/scheduler.cxx index c43c11b4..382a2164 100644 --- a/src/Airports/GenAirports850/scheduler.cxx +++ b/src/Airports/GenAirports850/scheduler.cxx @@ -316,7 +316,7 @@ void Scheduler::RetryAirport( AirportInfo* pai ) // retryList.push_back( *pai ); } -bool Scheduler::AddAirports( long start_pos, tg::Rectangle* boundingBox ) +bool Scheduler::AddAirports( long start_pos, tgRectangle* boundingBox ) { char line[2048]; char* def; diff --git a/src/Airports/GenAirports850/scheduler.hxx b/src/Airports/GenAirports850/scheduler.hxx index 8279301c..c68917dd 100644 --- a/src/Airports/GenAirports850/scheduler.hxx +++ b/src/Airports/GenAirports850/scheduler.hxx @@ -101,7 +101,7 @@ public: long FindAirport( std::string icao ); void AddAirport( std::string icao ); - bool AddAirports( long start_pos, tg::Rectangle* boundingBox ); + bool AddAirports( long start_pos, tgRectangle* boundingBox ); void RetryAirport( AirportInfo* pInfo ); void Schedule( int num_threads, std::string& summaryfile ); diff --git a/src/BuildTiles/Main/tgconstruct.hxx b/src/BuildTiles/Main/tgconstruct.hxx index 5384fddd..3ad9bcb4 100644 --- a/src/BuildTiles/Main/tgconstruct.hxx +++ b/src/BuildTiles/Main/tgconstruct.hxx @@ -193,11 +193,11 @@ public: void ConstructBucketStage3(); int load_landcover (); -// double measure_roughness( TGPolygon &poly ); + double measure_roughness( tgContour &contour ); AreaType get_landcover_type (const LandCover &cover, double xpos, double ypos, double dx, double dy); -// void make_area( const LandCover &cover, TGPolygon *polys, -// double x1, double y1, double x2, double y2, -// double half_dx, double half_dy ); + void make_area( const LandCover &cover, tgpolygon_list& polys, + double x1, double y1, double x2, double y2, + double half_dx, double half_dy ); // land cover file inline std::string get_cover () const { return cover; } diff --git a/src/BuildTiles/Main/tgconstruct_cleanup.cxx b/src/BuildTiles/Main/tgconstruct_cleanup.cxx index 349b704b..351afcb2 100644 --- a/src/BuildTiles/Main/tgconstruct_cleanup.cxx +++ b/src/BuildTiles/Main/tgconstruct_cleanup.cxx @@ -25,19 +25,15 @@ # include #endif -//#include - -//#include +#include #include #include "tgconstruct.hxx" -using std::string; - void TGConstruct::FixTJunctions( void ) { int before, after; std::vector points; - tg::Rectangle bb; + tgRectangle bb; // traverse each poly, and add intermediate nodes for ( unsigned int i = 0; i < TG_MAX_AREA_TYPES; ++i ) { diff --git a/src/BuildTiles/Main/tgconstruct_elevation.cxx b/src/BuildTiles/Main/tgconstruct_elevation.cxx index 35dc5081..52b10530 100644 --- a/src/BuildTiles/Main/tgconstruct_elevation.cxx +++ b/src/BuildTiles/Main/tgconstruct_elevation.cxx @@ -50,14 +50,14 @@ void TGConstruct::LoadElevationArray( bool add_nodes ) { array.remove_voids( ); if ( add_nodes ) { - point_list corner_list = array.get_corner_list(); + std::vector const& corner_list = array.get_corner_list(); for (unsigned int i=0; i const& fit_list = array.get_fitted_list(); for (unsigned int i=0; i +#include #include "tgconstruct.hxx" #include "usgs.hxx" @@ -46,66 +47,61 @@ static const double quarter_cover_size = cover_size * 0.25; // make the area specified area, look up the land cover type, and add // it to polys -#if 0 -void TGConstruct::make_area( const LandCover &cover, TGPolygon *polys, - double x1, double y1, double x2, double y2, - double half_dx, double half_dy ) +void TGConstruct::make_area( const LandCover &cover, tgpolygon_list& polys, + double x1, double y1, double x2, double y2, + double half_dx, double half_dy ) { const double fudge = 0.0001; // (0.0001 degrees =~ 10 meters) AreaType area = get_landcover_type( cover, - x1 + half_dx, y1 + half_dy, - x2 - x1, y2 - y1 ); + x1 + half_dx, y1 + half_dy, + x2 - x1, y2 - y1 ); if ( area != get_default_area_type() ) { // Create a square polygon and merge it into the list. - TGPolygon poly; - poly.erase(); - poly.add_node(0, Point3D(x1 - fudge, y1 - fudge, 0.0)); - poly.add_node(0, Point3D(x1 - fudge, y2 + fudge, 0.0)); - poly.add_node(0, Point3D(x2 + fudge, y2 + fudge, 0.0)); - poly.add_node(0, Point3D(x2 + fudge, y1 - fudge, 0.0)); + tgContour contour; - if ( measure_roughness( poly ) < 1.0 ) { - // Add this poly to the area accumulator - if ( polys[area].contours() > 0 ) { - polys[area] = tgPolygonUnion( polys[area], poly ); + contour.Erase(); + contour.AddNode( SGGeod::fromDeg(x1 - fudge, y1 - fudge) ); + contour.AddNode( SGGeod::fromDeg(x1 - fudge, y2 + fudge) ); + contour.AddNode( SGGeod::fromDeg(x2 + fudge, y2 + fudge) ); + contour.AddNode( SGGeod::fromDeg(x2 + fudge, y1 - fudge) ); + contour.SetHole( false ); + + if ( measure_roughness( contour ) < 1.0 ) { + // Add this contour to the area accumulator + if ( polys[area].Contours() > 0 ) { + polys[area] = tgContour::Union( contour, polys[area] ); } else { + tgPolygon poly; + poly.AddContour( contour ); polys[area] = poly; } } } } -#endif // Come up with a "rough" metric for the roughness of the terrain // coverted by a polygon -#if 0 -double TGConstruct::measure_roughness( TGPolygon &poly ) { - int i; - unsigned int j; - +double TGConstruct::measure_roughness( tgContour &contour ) { // find the elevation range double max_z = -9999.0; - double min_z = 9999.0; + double min_z = 9999.0; - for ( i = 0; i < poly.contours(); ++i ) { - point_list points = poly.get_contour( i ); - for ( j = 0; j < points.size(); ++j ) { - double z; - z = array.altitude_from_grid( points[j].x() * 3600.0, - points[j].y() * 3600.0 ); - if ( z < -9000 ) { - z = array.closest_nonvoid_elev( points[j].x() * 3600.0, - points[j].y() * 3600.0 ); - } + for ( unsigned int i = 0; i < contour.GetSize(); i++ ) { + double z; + z = array.altitude_from_grid( contour[i].getLongitudeDeg() * 3600.0, + contour[i].getLatitudeDeg() * 3600.0 ); + if ( z < -9000 ) { + z = array.closest_nonvoid_elev( contour[i].getLongitudeDeg() * 3600.0, + contour[i].getLatitudeDeg() * 3600.0 ); + } - if ( z < min_z ) { - min_z = z; - } - if ( z > max_z ) { - max_z = z; - } + if ( z < min_z ) { + min_z = z; + } + if ( z > max_z ) { + max_z = z; } } @@ -119,7 +115,6 @@ double TGConstruct::measure_roughness( TGPolygon &poly ) { return diff / 50.0; } -#endif AreaType TGConstruct::get_landcover_type (const LandCover &cover, double xpos, double ypos, double dx, double dy) { diff --git a/src/BuildTiles/Main/tgconstruct_tesselate.cxx b/src/BuildTiles/Main/tgconstruct_tesselate.cxx index 0d4d7a7f..78a0a1fb 100644 --- a/src/BuildTiles/Main/tgconstruct_tesselate.cxx +++ b/src/BuildTiles/Main/tgconstruct_tesselate.cxx @@ -43,7 +43,7 @@ void TGConstruct::TesselatePolys( void ) tgPolygon::ToShapefile( poly, ds_name, "preteselate", "" ); } - tg::Rectangle rect = poly.GetBoundingBox(); + tgRectangle rect = poly.GetBoundingBox(); nodes.get_geod_inside( rect.getMin(), rect.getMax(), poly_extra ); SG_LOG( SG_CLIPPER, SG_DEBUG, "Tesselating " << get_area_name( (AreaType)area ) << "(" << area << "): " << @@ -59,7 +59,6 @@ void TGConstruct::TesselatePolys( void ) } } - for (unsigned int area = 0; area < TG_MAX_AREA_TYPES; area++) { for (unsigned int p = 0; p < polys_clipped.area_size(area); p++ ) { tgPolygon& poly = polys_clipped.get_poly(area, p ); diff --git a/src/Lib/Array/array.cxx b/src/Lib/Array/array.cxx index 737ebcdc..599a02ff 100644 --- a/src/Lib/Array/array.cxx +++ b/src/Lib/Array/array.cxx @@ -140,8 +140,8 @@ TGArray::parse( SGBucket& b ) { *fitted_in >> fitted_size; for ( int i = 0; i < fitted_size; ++i ) { *fitted_in >> x >> y >> z; - fitted_list.push_back( Point3D(x, y, z) ); - SG_LOG(SG_GENERAL, SG_DEBUG, " loading fitted = " << Point3D(x, y, z) ); + fitted_list.push_back( SGGeod::fromDegM(x, y, z) ); + SG_LOG(SG_GENERAL, SG_DEBUG, " loading fitted = " << SGGeod::fromDegM(x, y, z) ); } } @@ -298,18 +298,16 @@ void TGArray::remove_voids( ) { double TGArray::closest_nonvoid_elev( double lon, double lat ) const { double mindist = 99999999999.9; double minelev = -9999.0; - Point3D p0( lon, lat, 0.0 ); + SGGeod p0 = SGGeod::fromDeg( lon, lat ); for ( int row = 0; row < rows; row++ ) { for ( int col = 0; col < cols; col++ ) { - Point3D p1(originx + col * col_step, originy + row * row_step, 0.0); - double dist = p0.distance3D( p1 ); + SGGeod p1 = SGGeod::fromDeg( originx + col * col_step, originy + row * row_step ); + double dist = SGGeodesy::distanceM( p0, p1 ); double elev = get_array_elev(col, row); if ( dist < mindist && elev > -9000 ) { mindist = dist; minelev = elev; - // cout << "dist = " << mindist; - // cout << " elev = " << elev << endl; } } } diff --git a/src/Lib/Array/array.hxx b/src/Lib/Array/array.hxx index 2a249e4a..85514393 100644 --- a/src/Lib/Array/array.hxx +++ b/src/Lib/Array/array.hxx @@ -28,8 +28,6 @@ #include #include -#include - class TGArray { private: @@ -51,8 +49,8 @@ private: short *in_data; // output nodes - point_list corner_list; - point_list fitted_list; + std::vector corner_list; + std::vector fitted_list; void parse_bin(); public: @@ -99,12 +97,11 @@ public: inline double get_col_step() const { return col_step; } inline double get_row_step() const { return row_step; } - inline point_list get_corner_list() const { return corner_list; } - inline point_list get_fitted_list() const { return fitted_list; } + inline std::vector const& get_corner_list() const { return corner_list; } + inline std::vector const& get_fitted_list() const { return fitted_list; } int get_array_elev( int col, int row ) const; void set_array_elev( int col, int row, int val ); }; - #endif // _ARRAY_HXX diff --git a/src/Lib/Polygon/CMakeLists.txt b/src/Lib/Polygon/CMakeLists.txt index 7791b411..4d5ba863 100644 --- a/src/Lib/Polygon/CMakeLists.txt +++ b/src/Lib/Polygon/CMakeLists.txt @@ -9,6 +9,8 @@ add_library(Polygon STATIC polygon.hxx rectangle.cxx rectangle.hxx + surface.cxx + surface.hxx tg_nodes.cxx tg_nodes.hxx tg_unique_geod.hxx @@ -16,5 +18,4 @@ add_library(Polygon STATIC tg_unique_vec2f.hxx tg_unique_vec3d.hxx tg_unique_vec3f.hxx - point3d.hxx ) \ No newline at end of file diff --git a/src/Lib/Polygon/point3d.hxx b/src/Lib/Polygon/point3d.hxx deleted file mode 100644 index 92f3730a..00000000 --- a/src/Lib/Polygon/point3d.hxx +++ /dev/null @@ -1,569 +0,0 @@ -/** - * \file point3d.hxx - * A 3d point class (depricated). This class is depricated and we are - * in the process of removing all usage of it in favor of plib's "sg" - * library of point, vector, and math routines. Plib's sg lib is less - * object oriented, but integrates more seamlessly with opengl. - * - * Adapted from algebra3 by Jean-Francois Doue, started October 1998. - */ - -// Copyright (C) 1998 Curtis L. Olson - http://www.flightgear.org/~curt -// -// This library is free software; you can redistribute it and/or -// modify it under the terms of the GNU Library General Public -// License as published by the Free Software Foundation; either -// version 2 of the License, or (at your option) any later version. -// -// This library 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 -// Library 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 Street, Fifth Floor, Boston, MA 02110-1301, USA. -// -// $Id$ - - -#ifndef _POINT3D_HXX -#define _POINT3D_HXX - - -#ifndef __cplusplus -# error This library requires C++ -#endif - -#include - -#include -#include -#include -#include -#include -#include - -#include - -#include - -//const double fgPoint3_Epsilon = 0.0000001; -const double fgPoint3_Epsilon = 0.000001; - -enum {PX, PY, PZ}; // axes - -// Kludge for msvc++ 6.0 - requires forward decls of friend functions. -class Point3D; -std::istream& operator>> ( std::istream&, Point3D& ); -std::ostream& operator<< ( std::ostream&, const Point3D& ); -Point3D operator- (const Point3D& p); // -p1 -bool operator== (const Point3D& a, const Point3D& b); // p1 == p2? - - -/** - * 3D Point class. - */ - -class Point3D { - -protected: - - double n[3]; - -public: - - /** Default constructor */ - Point3D(); - Point3D(const double x, const double y, const double z); - explicit Point3D(const double d); - Point3D(const Point3D &p); - - static Point3D fromSGGeod(const SGGeod& geod); - static Point3D fromSGGeoc(const SGGeoc& geoc); - static Point3D fromSGVec3(const SGVec3& cart); - static Point3D fromSGVec3(const SGVec3& cart); - static Point3D fromSGVec2(const SGVec2& cart); - static Point3D fromSGVec2(const SGVec2& cart); - - // Assignment operators - - Point3D& operator = ( const Point3D& p ); // assignment of a Point3D - Point3D& operator += ( const Point3D& p ); // incrementation by a Point3D - Point3D& operator -= ( const Point3D& p ); // decrementation by a Point3D - Point3D& operator *= ( const double d ); // multiplication by a constant - Point3D& operator /= ( const double d ); // division by a constant - - bool operator > ( const Point3D& p ); - bool operator < ( const Point3D& p ); - - void setx(const double x); - void sety(const double y); - void setz(const double z); - void setlon(const double x); - void setlat(const double y); - void setradius(const double z); - void setelev(const double z); - - void snap( double grid ); - - // Queries - - double& operator [] ( int i); // indexing - double operator[] (int i) const; // read-only indexing - - inline const double *get_n() const { return n; }; - double x() const; // cartesian x - double y() const; // cartesian y - double z() const; // cartesian z - - double lon() const; // polar longitude - double lat() const; // polar latitude - double radius() const; // polar radius - double elev() const; // geodetic elevation (if specifying a surface point) - - SGGeod toSGGeod(void) const; - SGGeoc toSGGeoc(void) const; - - SGVec3d toSGVec3d(void) const; - SGVec3f toSGVec3f(void) const; - SGVec2f toSGVec2f(void) const; - - // friends - friend Point3D operator - (const Point3D& p); // -p1 - friend bool operator == (const Point3D& a, const Point3D& b); // p1 == p2? - friend std::istream& operator>> ( std::istream&, Point3D& ); - friend std::ostream& operator<< ( std::ostream&, const Point3D& ); - - // Special functions - double distance3D(const Point3D& a) const; // distance between - double distance3Dsquared(const Point3D& a) const; // distance between ^ 2 - - bool IsEqual2D(const Point3D& a) const; // equality check in X,Y only - bool HasElevation() const; // does point have elevation data? - - bool IsWithin( Point3D min, Point3D max ) const; - bool IsWithin( double xmin, double xmax, double ymin, double ymax ) const; - bool IsAlmostWithin( Point3D min, Point3D max ) const; - bool IsAlmostWithin( double xmin, double xmax, double ymin, double ymax ) const; - -}; - - -// input from stream -inline std::istream& operator >> ( std::istream& in, Point3D& p) -{ - in >> p.n[PX]; - in >> p.n[PY]; - in >> p.n[PZ]; - - return in; -} - -inline std::ostream& operator<< ( std::ostream& out, const Point3D& p ) -{ - out << p.n[PX] << " "; - out << p.n[PY] << " "; - out << p.n[PZ] << "\n"; - - return out; -} - -inline void sgWritePoint3D ( gzFile fd, const Point3D& var ) { - sgWriteDouble ( fd, 3, var.get_n() ) ; -} - -inline void sgReadPoint3D ( gzFile fd, Point3D& var ) { - double data[3]; - sgReadDouble ( fd, 3, data ); - var = Point3D(data[0], data[1], data[2]); -} - -/////////////////////////// -// -// Point3D Member functions -// -/////////////////////////// - -// CONSTRUCTORS - -inline Point3D::Point3D() -{ - n[PX] = n[PY] = 0.0; - n[PZ] = -9999.0; -} - -inline Point3D::Point3D(const double x, const double y, const double z) -{ - n[PX] = x; n[PY] = y; n[PZ] = z; -} - -inline Point3D::Point3D(const double d) -{ - n[PX] = n[PY] = n[PZ] = d; -} - -inline Point3D::Point3D(const Point3D& p) -{ - n[PX] = p.n[PX]; n[PY] = p.n[PY]; n[PZ] = p.n[PZ]; -} - -inline Point3D Point3D::fromSGGeod(const SGGeod& geod) -{ - Point3D pt; - pt.setlon(geod.getLongitudeDeg()); - pt.setlat(geod.getLatitudeDeg()); - pt.setelev(geod.getElevationM()); - - return pt; -} - -inline Point3D Point3D::fromSGGeoc(const SGGeoc& geoc) -{ - Point3D pt; - pt.setlon(geoc.getLongitudeRad()); - pt.setlat(geoc.getLatitudeRad()); - pt.setradius(geoc.getRadiusM()); - - return pt; -} - -inline Point3D Point3D::fromSGVec3(const SGVec3& cart) -{ - Point3D pt; - pt.setx(cart.x()); - pt.sety(cart.y()); - pt.setz(cart.z()); - - return pt; -} - -inline Point3D Point3D::fromSGVec3(const SGVec3& cart) -{ - Point3D pt; - pt.setx(cart.x()); - pt.sety(cart.y()); - pt.setz(cart.z()); - - return pt; -} - -inline Point3D Point3D::fromSGVec2(const SGVec2& cart) -{ - Point3D pt; - pt.setx(cart.x()); - pt.sety(cart.y()); - pt.setz(0); - - return pt; -} - -inline Point3D Point3D::fromSGVec2(const SGVec2& cart) -{ - Point3D pt; - pt.setx(cart.x()); - pt.sety(cart.y()); - pt.setz(0); - - return pt; -} - - -// ASSIGNMENT OPERATORS - -inline Point3D& Point3D::operator = (const Point3D& p) -{ - n[PX] = p.n[PX]; n[PY] = p.n[PY]; n[PZ] = p.n[PZ]; - - return *this; -} - -inline Point3D& Point3D::operator += ( const Point3D& p ) -{ - n[PX] += p.n[PX]; n[PY] += p.n[PY]; n[PZ] += p.n[PZ]; - - return *this; -} - -inline Point3D& Point3D::operator -= ( const Point3D& p ) -{ - n[PX] -= p.n[PX]; n[PY] -= p.n[PY]; n[PZ] -= p.n[PZ]; - - return *this; -} - -inline Point3D& Point3D::operator *= ( const double d ) -{ - n[PX] *= d; n[PY] *= d; n[PZ] *= d; - - return *this; -} - -inline Point3D& Point3D::operator /= ( const double d ) -{ - double d_inv = 1./d; n[PX] *= d_inv; n[PY] *= d_inv; n[PZ] *= d_inv; - - return *this; -} - -inline bool Point3D::operator < ( const Point3D& p ) -{ - return ( (n[PX] < p.n[PX]) && (n[PY] < p.n[PY]) ); -} - -inline bool Point3D::operator > ( const Point3D& p ) -{ - return ( (n[PX] > p.n[PX]) && (n[PY] > p.n[PY]) ); -} - -inline void Point3D::setx(const double x) { - n[PX] = x; -} - -inline void Point3D::sety(const double y) { - n[PY] = y; -} - -inline void Point3D::setz(const double z) { - n[PZ] = z; -} - -inline void Point3D::setlon(const double x) { - n[PX] = x; -} - -inline void Point3D::setlat(const double y) { - n[PY] = y; -} - -inline void Point3D::setradius(const double z) { - n[PZ] = z; -} - -inline void Point3D::setelev(const double z) { - n[PZ] = z; -} - -inline void Point3D::snap( double grid ) -{ - n[PX] = grid * SGMisc::round( n[PX]/grid ); - n[PY] = grid * SGMisc::round( n[PY]/grid ); - n[PZ] = grid * SGMisc::round( n[PZ]/grid ); -} - -// QUERIES - -inline double& Point3D::operator [] ( int i) -{ - assert(! (i < PX || i > PZ)); - return n[i]; -} - -inline double Point3D::operator [] ( int i) const { - assert(! (i < PX || i > PZ)); - return n[i]; -} - - -inline double Point3D::x() const { return n[PX]; } - -inline double Point3D::y() const { return n[PY]; } - -inline double Point3D::z() const { return n[PZ]; } - -inline double Point3D::lon() const { return n[PX]; } - -inline double Point3D::lat() const { return n[PY]; } - -inline double Point3D::radius() const { return n[PZ]; } - -inline double Point3D::elev() const { return n[PZ]; } - -inline SGGeod Point3D::toSGGeod(void) const -{ - SGGeod geod; - geod.setLongitudeDeg(lon()); - geod.setLatitudeDeg(lat()); - geod.setElevationM(elev()); - return geod; -} - -inline SGGeoc Point3D::toSGGeoc(void) const -{ - SGGeoc geoc; - geoc.setLongitudeRad(lon()); - geoc.setLatitudeRad(lat()); - geoc.setRadiusM(radius()); - return geoc; -} - -inline SGVec3d Point3D::toSGVec3d(void) const -{ - return SGVec3d(x(), y(), z()); -} - -inline SGVec3f Point3D::toSGVec3f(void) const -{ - return SGVec3f(x(), y(), z()); -} - -inline SGVec2f Point3D::toSGVec2f(void) const -{ - return SGVec2f(x(), y()); -} - -inline bool Point3D::IsEqual2D(const Point3D& a) const -{ - return - fabs(a.n[PX] - n[PX]) < fgPoint3_Epsilon && - fabs(a.n[PY] - n[PY]) < fgPoint3_Epsilon; -} - -inline bool Point3D::HasElevation() const -{ - return ( fabs( n[PZ] + 9999.0 ) > fgPoint3_Epsilon ); -} - -inline bool Point3D::IsWithin( Point3D min, Point3D max ) const -{ - return ( (min.n[PX] <= n[PX]) && (min.n[PY] <= n[PY]) && - (max.n[PX] >= n[PX]) && (max.n[PY] >= n[PY]) ); -} - -inline bool Point3D::IsWithin( double xmin, double xmax, double ymin, double ymax ) const -{ - // make sure we take epsilon into account - xmin -= fgPoint3_Epsilon; - ymin -= fgPoint3_Epsilon; - - xmax += fgPoint3_Epsilon; - ymax += fgPoint3_Epsilon; - - return ( (xmin <= n[PX]) && (ymin <= n[PY]) && - (xmax >= n[PX]) && (ymax >= n[PY]) ); -} - -inline bool Point3D::IsAlmostWithin( Point3D min, Point3D max ) const -{ - // make sure we take epsilon into account - min.n[PX] -= fgPoint3_Epsilon; - min.n[PY] -= fgPoint3_Epsilon; - - max.n[PX] += fgPoint3_Epsilon; - max.n[PY] += fgPoint3_Epsilon; - - return ( IsWithin(min, max) ); -} - -inline bool Point3D::IsAlmostWithin( double xmin, double xmax, double ymin, double ymax ) const -{ - // make sure we take epsilon into account - xmin -= fgPoint3_Epsilon; - ymin -= fgPoint3_Epsilon; - - xmax += fgPoint3_Epsilon; - ymax += fgPoint3_Epsilon; - - return ( IsWithin( xmin, xmax, ymin,ymax ) ); -} - -// FRIENDS - -inline Point3D operator - (const Point3D& a) -{ - return Point3D(-a.n[PX],-a.n[PY],-a.n[PZ]); -} - -inline Point3D operator + (const Point3D& a, const Point3D& b) -{ - return Point3D(a) += b; -} - -inline Point3D operator - (const Point3D& a, const Point3D& b) -{ - return Point3D(a) -= b; -} - -inline Point3D operator * (const Point3D& a, const double d) -{ - return Point3D(a) *= d; -} - -inline Point3D operator * (const double d, const Point3D& a) -{ - Point3D pt = a*d; - - return pt; -} - -inline Point3D operator / (const Point3D& a, const double d) -{ - Point3D pt = Point3D(a) *= (1.0 / d ); - - return pt; -} - -inline bool operator == (const Point3D& a, const Point3D& b) -{ - return - fabs(a.n[PX] - b.n[PX]) < fgPoint3_Epsilon && - fabs(a.n[PY] - b.n[PY]) < fgPoint3_Epsilon && - fabs(a.n[PZ] - b.n[PZ]) < fgPoint3_Epsilon; -} - -inline bool operator != (const Point3D& a, const Point3D& b) -{ - return !(a == b); -} - -// Special functions - -inline double -Point3D::distance3D(const Point3D& a ) const -{ - double x, y, z; - - x = n[PX] - a.n[PX]; - y = n[PY] - a.n[PY]; - z = n[PZ] - a.n[PZ]; - - return sqrt(x*x + y*y + z*z); -} - - -inline double -Point3D::distance3Dsquared(const Point3D& a ) const -{ - double x, y, z; - - x = n[PX] - a.n[PX]; - y = n[PY] - a.n[PY]; - z = n[PZ] - a.n[PZ]; - - return(x*x + y*y + z*z); -} - - -typedef std::vector< Point3D > point_list; - -typedef point_list::iterator point_list_iterator; - -typedef point_list::const_iterator const_point_list_iterator; - -inline Point3D sgCartToGeod( const Point3D& p ) -{ - SGGeod geod; - SGGeodesy::SGCartToGeod(SGVec3d(p.x(), p.y(), p.z()), geod); - return Point3D::fromSGGeod(geod); -} - -inline Point3D sgGeodToCart(const Point3D& geod) -{ - SGVec3 cart; - SGGeodesy::SGGeodToCart(SGGeod::fromRadM(geod.lon(), geod.lat(), geod.elev()), cart); - return Point3D::fromSGVec3(cart); -} - -#endif // _POINT3D_HXX - - diff --git a/src/Lib/Polygon/polygon.cxx b/src/Lib/Polygon/polygon.cxx index fe20481b..d11618c7 100644 --- a/src/Lib/Polygon/polygon.cxx +++ b/src/Lib/Polygon/polygon.cxx @@ -30,14 +30,13 @@ #include #include #include +#include #include #include #include #include #include -#include - #include "polygon.hxx" #ifdef _MSC_VER @@ -159,87 +158,6 @@ double SGGeod_CalculateTheta( const SGGeod& p0, const SGGeod& p1, const SGGeod& return acos( uv_dot / (udist * vdist) ); } -bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, - const point_list& nodes, SGGeod& result, - double bbEpsilon, double errEpsilon ) -{ - bool found_node = false; - double m, m1, b, b1, y_err, x_err, y_err_min, x_err_min; - - SGGeod p0 = start; - SGGeod p1 = end; - - double xdist = fabs(p0.getLongitudeDeg() - p1.getLongitudeDeg()); - double ydist = fabs(p0.getLatitudeDeg() - p1.getLatitudeDeg()); - - x_err_min = xdist + 1.0; - y_err_min = ydist + 1.0; - - if ( xdist > ydist ) { - // sort these in a sensible order - SGGeod p_min, p_max; - if ( p0.getLongitudeDeg() < p1.getLongitudeDeg() ) { - p_min = p0; - p_max = p1; - } else { - p_min = p1; - p_max = p0; - } - - m = (p_min.getLatitudeDeg() - p_max.getLatitudeDeg()) / (p_min.getLongitudeDeg() - p_max.getLongitudeDeg()); - b = p_max.getLatitudeDeg() - m * p_max.getLongitudeDeg(); - - for ( int i = 0; i < (int)nodes.size(); ++i ) { - // cout << i << endl; - SGGeod current = nodes[i].toSGGeod(); - - if ( (current.getLongitudeDeg() > (p_min.getLongitudeDeg() + (bbEpsilon))) && (current.getLongitudeDeg() < (p_max.getLongitudeDeg() - (bbEpsilon))) ) { - y_err = fabs(current.getLatitudeDeg() - (m * current.getLongitudeDeg() + b)); - - if ( y_err < errEpsilon ) { - found_node = true; - if ( y_err < y_err_min ) { - result = current; - y_err_min = y_err; - } - } - } - } - } else { - // sort these in a sensible order - SGGeod p_min, p_max; - if ( p0.getLatitudeDeg() < p1.getLatitudeDeg() ) { - p_min = p0; - p_max = p1; - } else { - p_min = p1; - p_max = p0; - } - - m1 = (p_min.getLongitudeDeg() - p_max.getLongitudeDeg()) / (p_min.getLatitudeDeg() - p_max.getLatitudeDeg()); - b1 = p_max.getLongitudeDeg() - m1 * p_max.getLatitudeDeg(); - - for ( int i = 0; i < (int)nodes.size(); ++i ) { - SGGeod current = nodes[i].toSGGeod(); - - if ( (current.getLatitudeDeg() > (p_min.getLatitudeDeg() + (bbEpsilon))) && (current.getLatitudeDeg() < (p_max.getLatitudeDeg() - (bbEpsilon))) ) { - - x_err = fabs(current.getLongitudeDeg() - (m1 * current.getLatitudeDeg() + b1)); - - if ( x_err < errEpsilon ) { - found_node = true; - if ( x_err < x_err_min ) { - result = current; - x_err_min = x_err; - } - } - } - } - } - - return found_node; -} - bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, const std::vector& nodes, SGGeod& result, double bbEpsilon, double errEpsilon ) @@ -321,24 +239,6 @@ bool FindIntermediateNode( const SGGeod& start, const SGGeod& end, return found_node; } -void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, point_list& tmp_nodes, tgContour& result, double bbEpsilon, double errEpsilon ) -{ - SGGeod new_pt; - - SG_LOG(SG_GENERAL, SG_BULK, " " << p0 << " <==> " << p1 ); - - bool found_extra = FindIntermediateNode( p0, p1, tmp_nodes, new_pt, bbEpsilon, errEpsilon ); - - if ( found_extra ) { - AddIntermediateNodes( p0, new_pt, tmp_nodes, result, bbEpsilon, errEpsilon ); - - result.AddNode( new_pt ); - SG_LOG(SG_GENERAL, SG_BULK, " adding = " << new_pt); - - AddIntermediateNodes( new_pt, p1, tmp_nodes, result, bbEpsilon, errEpsilon ); - } -} - void AddIntermediateNodes( const SGGeod& p0, const SGGeod& p1, std::vector& nodes, tgContour& result, double bbEpsilon, double errEpsilon ) { SGGeod new_pt; @@ -385,7 +285,7 @@ static double Dist_ToClipper( double dist ) return ( dist * ( CLIPPER_FIXEDPT / CLIPPER_FIXED1M ) ); } -static tg::Rectangle BoundingBox_FromClipper( const ClipperLib::Polygons& subject ) +static tgRectangle BoundingBox_FromClipper( const ClipperLib::Polygons& subject ) { ClipperLib::IntPoint min_pt, max_pt; SGGeod min, max; @@ -417,7 +317,7 @@ static tg::Rectangle BoundingBox_FromClipper( const ClipperLib::Polygons& subjec min = SGGeod_FromClipper( min_pt ); max = SGGeod_FromClipper( max_pt ); - return tg::Rectangle( min, max ); + return tgRectangle( min, max ); } @@ -848,7 +748,7 @@ tgContour tgContour::Expand( const tgContour& subject, double offset ) return result; } -tg::Rectangle tgContour::GetBoundingBox( void ) const +tgRectangle tgContour::GetBoundingBox( void ) const { SGGeod min, max; @@ -868,7 +768,33 @@ tg::Rectangle tgContour::GetBoundingBox( void ) const min = SGGeod::fromDeg( minx, miny ); max = SGGeod::fromDeg( maxx, maxy ); - return tg::Rectangle( min, max ); + return tgRectangle( min, max ); +} + +tgPolygon tgContour::Union( const tgContour& subject, tgPolygon& clip ) +{ + tgPolygon result; + UniqueSGGeodSet all_nodes; + + /* before diff - gather all nodes */ + for ( unsigned int i = 0; i < subject.GetSize(); ++i ) { + all_nodes.add( subject.GetNode(i) ); + } + + ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); + ClipperLib::Polygons clipper_clip = tgPolygon::ToClipper( clip ); + ClipperLib::Polygons clipper_result; + + ClipperLib::Clipper c; + c.Clear(); + c.AddPolygon(clipper_subject, ClipperLib::ptSubject); + c.AddPolygons(clipper_clip, ClipperLib::ptClip); + c.Execute(ClipperLib::ctUnion, clipper_result, ClipperLib::pftEvenOdd, ClipperLib::pftEvenOdd); + + result = tgPolygon::FromClipper( clipper_result ); + result = tgPolygon::AddColinearNodes( result, all_nodes ); + + return result; } tgPolygon tgContour::Diff( const tgContour& subject, tgPolygon& clip ) @@ -1720,7 +1646,7 @@ tgPolygon tgPolygon::Intersect( const tgPolygon& subject, const tgPolygon& clip return result; } -tg::Rectangle tgPolygon::GetBoundingBox( void ) const +tgRectangle tgPolygon::GetBoundingBox( void ) const { SGGeod min, max; @@ -1742,7 +1668,7 @@ tg::Rectangle tgPolygon::GetBoundingBox( void ) const min = SGGeod::fromDeg( minx, miny ); max = SGGeod::fromDeg( maxx, maxy ); - return tg::Rectangle( min, max ); + return tgRectangle( min, max ); } void clipperToShapefile( ClipperLib::Polygons polys, const std::string& path, const std::string&layer, const std::string& name ) @@ -2460,7 +2386,7 @@ tgPolygon tgAccumulator::Diff( const tgContour& subject ) } unsigned int num_hits = 0; - tg::Rectangle box1 = subject.GetBoundingBox(); + tgRectangle box1 = subject.GetBoundingBox(); ClipperLib::Polygon clipper_subject = tgContour::ToClipper( subject ); ClipperLib::Polygons clipper_result; @@ -2472,7 +2398,7 @@ tgPolygon tgAccumulator::Diff( const tgContour& subject ) // clip result against all polygons in the accum that intersect our bb for (unsigned int i=0; i < accum.size(); i++) { - tg::Rectangle box2 = BoundingBox_FromClipper( accum[i] ); + tgRectangle box2 = BoundingBox_FromClipper( accum[i] ); if ( box2.intersects(box1) ) { @@ -2529,7 +2455,7 @@ tgPolygon tgAccumulator::Diff( const tgPolygon& subject ) } unsigned int num_hits = 0; - tg::Rectangle box1 = subject.GetBoundingBox(); + tgRectangle box1 = subject.GetBoundingBox(); ClipperLib::Polygons clipper_subject = tgPolygon::ToClipper( subject ); ClipperLib::Polygons clipper_result; @@ -2541,7 +2467,7 @@ tgPolygon tgAccumulator::Diff( const tgPolygon& subject ) // clip result against all polygons in the accum that intersect our bb for (unsigned int i=0; i < accum.size(); i++) { - tg::Rectangle box2 = BoundingBox_FromClipper( accum[i] ); + tgRectangle box2 = BoundingBox_FromClipper( accum[i] ); if ( box2.intersects(box1) ) { @@ -2674,14 +2600,12 @@ void tgChopper::Clip( const tgPolygon& subject, const std::string& type, SGBucket& b ) { - Point3D p; + // 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) ) { @@ -2730,30 +2654,6 @@ void tgChopper::Clip( const tgPolygon& subject, } 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(); @@ -2762,7 +2662,7 @@ void tgChopper::Clip( const tgPolygon& subject, void tgChopper::Add( const tgPolygon& subject, const std::string& type ) { - tg::Rectangle bb; + tgRectangle bb; SGGeod p; // bail out immediately if polygon is empty diff --git a/src/Lib/Polygon/polygon.hxx b/src/Lib/Polygon/polygon.hxx index a268e7a6..ad2226d0 100644 --- a/src/Lib/Polygon/polygon.hxx +++ b/src/Lib/Polygon/polygon.hxx @@ -139,7 +139,7 @@ public: } } - tg::Rectangle GetBoundingBox( void ) const; + tgRectangle GetBoundingBox( void ) const; double GetMinimumAngle( void ) const; double GetArea( void ) const; @@ -150,6 +150,7 @@ public: static tgContour SplitLongEdges( const tgContour& subject, double dist ); static tgContour RemoveSpikes( const tgContour& subject ); + static tgPolygon Union( const tgContour& subject, tgPolygon& clip ); static tgPolygon Diff( const tgContour& subject, tgPolygon& clip ); static tgContour AddColinearNodes( const tgContour& subject, UniqueSGGeodSet& nodes ); @@ -350,7 +351,7 @@ public: contours[c].AddNode( n ); } - tg::Rectangle GetBoundingBox( void ) const; + tgRectangle GetBoundingBox( void ) const; void InheritElevations( const tgPolygon& source ); diff --git a/src/Lib/Polygon/rectangle.cxx b/src/Lib/Polygon/rectangle.cxx index af3ce1bf..b3d135e1 100644 --- a/src/Lib/Polygon/rectangle.cxx +++ b/src/Lib/Polygon/rectangle.cxx @@ -6,83 +6,64 @@ #include "rectangle.hxx" -namespace tg { - -Rectangle::Rectangle () +tgRectangle::tgRectangle() { } -Rectangle::Rectangle (const Rectangle &r) +tgRectangle::tgRectangle (const tgRectangle &r) : _min(r.getMin()), _max(r.getMax()) { } -Rectangle::Rectangle (const SGGeod &min, const SGGeod &max) +tgRectangle::tgRectangle (const SGGeod &min, const SGGeod &max) : _min(min), _max(max) { } -Rectangle::~Rectangle () +tgRectangle::~tgRectangle () { } -void -Rectangle::setMin (const SGGeod &p) +void tgRectangle::setMin (const SGGeod &p) { - _min = p; + _min = p; } -void -Rectangle::setMax (const SGGeod &p) +void tgRectangle::setMax (const SGGeod &p) { - _max = p; + _max = p; } -void -Rectangle::sanify () +void tgRectangle::sanify () { - double tmp; - if (_min.getLongitudeDeg() > _max.getLongitudeDeg()) { - tmp = _min.getLongitudeDeg(); - _min.setLongitudeDeg(_max.getLongitudeDeg()); - _max.setLongitudeDeg(tmp); - } - if (_min.getLatitudeDeg() > _max.getLatitudeDeg()) { - tmp = _min.getLatitudeDeg(); - _min.setLatitudeDeg(_max.getLatitudeDeg()); - _max.setLatitudeDeg(tmp); - } + double tmp; + if (_min.getLongitudeDeg() > _max.getLongitudeDeg()) { + tmp = _min.getLongitudeDeg(); + _min.setLongitudeDeg(_max.getLongitudeDeg()); + _max.setLongitudeDeg(tmp); + } + if (_min.getLatitudeDeg() > _max.getLatitudeDeg()) { + tmp = _min.getLatitudeDeg(); + _min.setLatitudeDeg(_max.getLatitudeDeg()); + _max.setLatitudeDeg(tmp); + } } -bool -Rectangle::isInside (const SGGeod &p) const +bool tgRectangle::isInside (const SGGeod &p) const { return ((p.getLongitudeDeg() >= _min.getLongitudeDeg() && p.getLongitudeDeg() <= _max.getLongitudeDeg()) && - (p.getLatitudeDeg() >= _min.getLatitudeDeg() && p.getLatitudeDeg() <= _max.getLatitudeDeg())); + (p.getLatitudeDeg() >= _min.getLatitudeDeg() && p.getLatitudeDeg() <= _max.getLatitudeDeg())); } -bool -Rectangle::intersects (const Rectangle &r) const +bool tgRectangle::intersects (const tgRectangle &r) const { - const SGGeod &min = r.getMin(); - const SGGeod &max = r.getMax(); - return ((max.getLongitudeDeg() >= _min.getLongitudeDeg()) && (min.getLongitudeDeg() <= _max.getLongitudeDeg()) && - (max.getLatitudeDeg() >= _min.getLatitudeDeg()) && (min.getLatitudeDeg() <= _max.getLatitudeDeg())); + const SGGeod &min = r.getMin(); + const SGGeod &max = r.getMax(); + + return ((max.getLongitudeDeg() >= _min.getLongitudeDeg()) && (min.getLongitudeDeg() <= _max.getLongitudeDeg()) && + (max.getLatitudeDeg() >= _min.getLatitudeDeg()) && (min.getLatitudeDeg() <= _max.getLatitudeDeg())); } -/*const TGPolygon -Rectangle::toPoly () const -{ - TGPolygon poly; - poly.add_node(0, _min); - poly.add_node(0, Point3D(_max.x(), _min.y(), 0)); - poly.add_node(0, _max); - poly.add_node(0, Point3D(_min.x(), _max.y(), 0)); - return poly; -}*/ - -}; - // end of rectangle.cxx diff --git a/src/Lib/Polygon/rectangle.hxx b/src/Lib/Polygon/rectangle.hxx index 64e19bea..ab5fd1cd 100644 --- a/src/Lib/Polygon/rectangle.hxx +++ b/src/Lib/Polygon/rectangle.hxx @@ -14,9 +14,6 @@ #include #include - -namespace tg { - /** * A simple rectangle class for bounding rectangles. * @@ -25,113 +22,104 @@ namespace tg { * sanify the rectangle (to make certain that each point is correct) * and to test whether another point lies inside it. */ -class Rectangle +class tgRectangle { public: - /** - * Create a new empty rectangle with both points at 0,0. - */ - Rectangle (); + /** + * Create a new empty rectangle with both points at 0,0. + */ + tgRectangle(); - /** - * Copy an existing rectangle. - * - * @param r The rectangle to copy. - */ - Rectangle (const Rectangle &r); + /** + * Copy an existing rectangle. + * + * @param r The rectangle to copy. + */ + tgRectangle (const tgRectangle &r); - /** - * Convenience constructor. - */ - Rectangle (const SGGeod& min, const SGGeod& max); + /** + * Convenience constructor. + */ + tgRectangle (const SGGeod& min, const SGGeod& max); - /** - * Destructor. - */ - virtual ~Rectangle (); + /** + * Destructor. + */ + virtual ~tgRectangle(); - /** - * Get the minimum (top left) corner of the rectangle. - * - * @return The top-left vertex. - */ - virtual const SGGeod &getMin () const { return _min; } + /** + * Get the minimum (top left) corner of the rectangle. + * + * @return The top-left vertex. + */ + virtual const SGGeod &getMin () const { return _min; } - /** - * Get the maximum (bottom right) corner of the rectangle. - * - * @return The bottom-right vertex. - */ - virtual const SGGeod &getMax () const { return _max; } + /** + * Get the maximum (bottom right) corner of the rectangle. + * + * @return The bottom-right vertex. + */ + virtual const SGGeod &getMax () const { return _max; } - /** - * Get the minimum (top left) corner of the rectangle. - * - * @return The top-left vertex. - */ - virtual SGGeod &getMin () { return _min; } + /** + * Get the minimum (top left) corner of the rectangle. + * + * @return The top-left vertex. + */ + virtual SGGeod &getMin () { return _min; } - /** - * Get the maximum (bottom right) corner of the rectangle. - * - * @return The bottom-right vertex. - */ - virtual SGGeod &getMax () { return _max; } + /** + * Get the maximum (bottom right) corner of the rectangle. + * + * @return The bottom-right vertex. + */ + virtual SGGeod &getMax () { return _max; } - /** - * Set the minimum (top-left) corner of the rectangle. - * - * @param p The top-left vertex. - */ - virtual void setMin (const SGGeod& p); + /** + * Set the minimum (top-left) corner of the rectangle. + * + * @param p The top-left vertex. + */ + virtual void setMin (const SGGeod& p); - /** - * Set the maximum (bottom-right) corner of the rectangle. - * - * @param p The bottom-right vertex. - */ - virtual void setMax (const SGGeod& p); + /** + * Set the maximum (bottom-right) corner of the rectangle. + * + * @param p The bottom-right vertex. + */ + virtual void setMax (const SGGeod& p); - /** - * Make the rectangle sane. - * - * Ensure that the min vertex is less than the max vertex. - */ - virtual void sanify (); + /** + * Make the rectangle sane. + * + * Ensure that the min vertex is less than the max vertex. + */ + virtual void sanify (); - /** - * Test whether a point lies inside the rectangle. - * - * The z-coordinates are ignored. - * - * @param p The point to test. - * @return true if the point is inside or on the boundary of the - * rectangle, false if it is outside. - */ - virtual bool isInside (const SGGeod& p) const; + /** + * Test whether a point lies inside the rectangle. + * + * The z-coordinates are ignored. + * + * @param p The point to test. + * @return true if the point is inside or on the boundary of the + * rectangle, false if it is outside. + */ + virtual bool isInside (const SGGeod& p) const; - /** - * Test whether this rectangle overlaps with another one. - * - * @param r The rectangle to test. - * @return true if the rectangle is touching or overlapping, false - * otherwise. - */ - virtual bool intersects (const Rectangle &r) const; - - /** - * Create a polygon representation of this rectangle. - * - * @return A four-vertex polygon representing this rectangle. - */ -// virtual const TGPolygon toPoly () const; + /** + * Test whether this rectangle overlaps with another one. + * + * @param r The rectangle to test. + * @return true if the rectangle is touching or overlapping, false + * otherwise. + */ + virtual bool intersects (const tgRectangle &r) const; private: - SGGeod _min; - SGGeod _max; + SGGeod _min; + SGGeod _max; }; -}; - -#endif // __RECTANGLE_HXX +#endif // __TGRECTANGLE_HXX \ No newline at end of file diff --git a/src/Airports/GenAirports850/apt_surface.cxx b/src/Lib/Polygon/surface.cxx similarity index 55% rename from src/Airports/GenAirports850/apt_surface.cxx rename to src/Lib/Polygon/surface.cxx index 766ae4b1..15b6bbad 100644 --- a/src/Airports/GenAirports850/apt_surface.cxx +++ b/src/Lib/Polygon/surface.cxx @@ -24,21 +24,25 @@ # include #endif -#include - #include #include #include #include -#include "elevations.hxx" -#include "global.hxx" -#include "apt_surface.hxx" -#include "debug.hxx" +#include +#include "TNT/jama_qr.h" +#include "surface.hxx" -static bool limit_slope( SimpleMatrix *Pts, int i1, int j1, int i2, int j2, - double average_elev_m ) +// Final grid size for surface (in meters) +const double coarse_grid = 300.0; + +// compared to the average surface elevation, clamp all values within +// this many meters of the average +const double max_clamp = 100.0; + +static bool limit_slope( tgMatrix* Pts, int i1, int j1, int i2, int j2, + double average_elev_m, double slope_max, double slope_eps ) { bool slope_error = false; @@ -55,7 +59,7 @@ static bool limit_slope( SimpleMatrix *Pts, int i1, int j1, int i2, int j2, slope_error = true; - GENAPT_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of " << slope << " dist = " << dist ); + SG_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of " << slope << " dist = " << dist ); double e1 = fabs(average_elev_m - p1.getElevationM()); double e2 = fabs(average_elev_m - p2.getElevationM()); @@ -63,17 +67,17 @@ static bool limit_slope( SimpleMatrix *Pts, int i1, int j1, int i2, int j2, if ( e1 > e2 ) { // p1 error larger if ( slope > 0 ) { - p1.setElevationM( p2.getElevationM() - (dist * slope_max) ); + p1.setElevationM( p2.getElevationM() - (dist * slope_max) ); } else { - p1.setElevationM( p2.getElevationM() + (dist * slope_max) ); + p1.setElevationM( p2.getElevationM() + (dist * slope_max) ); } Pts->set(i1, j1, p1); } else { // p2 error larger if ( slope > 0 ) { - p2.setElevationM( p1.getElevationM() + (dist * slope_max) ); + p2.setElevationM( p1.getElevationM() + (dist * slope_max) ); } else { - p2.setElevationM( p1.getElevationM() - (dist * slope_max) ); + p2.setElevationM( p1.getElevationM() - (dist * slope_max) ); } Pts->set(i2, j2, p2); } @@ -83,12 +87,146 @@ static bool limit_slope( SimpleMatrix *Pts, int i1, int j1, int i2, int j2, } +// lookup node elevations for each point in the specified simple +// matrix. Returns average of all points. +static void tgCalcElevations( const std::string &root, const string_list elev_src, + tgMatrix &Pts, const double average ) +{ + bool done = false; + int i, j; + TGArray array; + + // just bail if no work to do + if ( Pts.rows() == 0 || Pts.cols() == 0 ) { + return; + } + + // set all elevations to -9999 + for ( j = 0; j < Pts.rows(); ++j ) { + for ( i = 0; i < Pts.cols(); ++i ) { + SGGeod p = Pts.element(i, j); + p.setElevationM(-9999.0); + Pts.set(i, j, p); + } + } + + while ( !done ) { + // find first node with -9999 elevation + SGGeod first = SGGeod(); + bool found_one = false; + for ( j = 0; j < Pts.rows(); ++j ) { + for ( i = 0; i < Pts.cols(); ++i ) { + SGGeod p = Pts.element(i,j); + if ( p.getElevationM() < -9000.0 && !found_one ) { + first = p; + found_one = true; + } + } + } + + if ( found_one ) { + SGBucket b( first ); + std::string base = b.gen_base_path(); + + // try the various elevation sources + j = 0; + bool found_file = false; + while ( !found_file && j < (int)elev_src.size() ) { + std::string array_path = root + "/" + elev_src[j] + "/" + base + "/" + b.gen_index_str(); + + if ( array.open(array_path) ) { + found_file = true; + SG_LOG( SG_GENERAL, SG_DEBUG, "Using array_path = " << array_path ); + } + j++; + } + + // this will fill in a zero structure if no array data + // found/opened + array.parse( b ); + + // this will do a hasty job of removing voids by inserting + // data from the nearest neighbor (sort of) + array.remove_voids(); + + // update all the non-updated elevations that are inside + // this array file + double elev; + done = true; + for ( j = 0; j < Pts.rows(); ++j ) { + for ( i = 0; i < Pts.cols(); ++i ) { + SGGeod p = Pts.element(i,j); + if ( p.getElevationM() < -9000.0 ) { + done = false; + elev = array.altitude_from_grid( p.getLongitudeDeg() * 3600.0, + p.getLatitudeDeg() * 3600.0 ); + if ( elev > -9000 ) { + p.setElevationM( elev ); + Pts.set(i, j, p); + } + } + } + } + + array.close(); + + } else { + done = true; + } + } + +#ifdef DEBUG + // do some post processing for sanity's sake + // find the average height of the queried points + double total = 0.0; + int count = 0; + for ( j = 0; j < Pts.rows(); ++j ) { + for ( i = 0; i < Pts.cols(); ++i ) { + SGGeod p = Pts.element(i,j); + total += p.getElevationM(); + count++; + } + } + double grid_average = total / (double) count; + SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of matrix = " << grid_average); +#endif +} + +// clamp all elevations to the specified range +static void tgClampElevations( tgMatrix& Pts, + double center_m, + double max_clamp_m ) +{ + int i, j; + + // go through the elevations and clamp all elevations to within + // +/-max_m of the center_m elevation. + for ( j = 0; j < Pts.rows(); ++j ) { + for ( i = 0; i < Pts.cols(); ++i ) { + SGGeod p = Pts.element(i,j); + if ( p.getElevationM() < center_m - max_clamp_m ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.getElevationM() << " to " << center_m - max_clamp_m ); + p.setElevationM( center_m - max_clamp_m ); + Pts.set(i, j, p); + } + if ( p.getElevationM() > center_m + max_clamp_m ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.getElevationM() << " to " << center_m + max_clamp_m ); + p.setElevationM( center_m + max_clamp_m ); + Pts.set(i, j, p); + } + } + } +} + // Constructor, specify min and max coordinates of desired area in // lon/lat degrees -TGAptSurface::TGAptSurface( const std::string& path, - const string_list& elev_src, - tg::Rectangle aptBounds, - double average_elev_m ) +tgSurface::tgSurface( const std::string& path, + const string_list& elev_src, + tgRectangle aptBounds, + double average_elev_m, + double slope_max, + double slope_eps + ) { // Calculate desired size of grid _aptBounds = aptBounds; @@ -110,8 +248,7 @@ TGAptSurface::TGAptSurface( const std::string& path, double x_nm = x_rad * SG_RAD_TO_NM * xfact; double x_m = x_nm * SG_NM_TO_METER; - GENAPT_LOG( SG_GENERAL, SG_DEBUG, - "Area size = " << y_m << " x " << x_m << " (m)" ); + SG_LOG( SG_GENERAL, SG_DEBUG, "Area size = " << y_m << " x " << x_m << " (m)" ); int xdivs = (int)(x_m / coarse_grid) + 1; int ydivs = (int)(y_m / coarse_grid) + 1; @@ -120,7 +257,8 @@ TGAptSurface::TGAptSurface( const std::string& path, // interesting if ( xdivs < 8 ) { xdivs = 8; } if ( ydivs < 8 ) { ydivs = 8; } - GENAPT_LOG(SG_GENERAL, SG_DEBUG, " M(" << ydivs << "," << xdivs << ")"); + + SG_LOG(SG_GENERAL, SG_DEBUG, " M(" << ydivs << "," << xdivs << ")"); double dlon = x_deg / xdivs; double dlat = y_deg / ydivs; @@ -131,15 +269,12 @@ TGAptSurface::TGAptSurface( const std::string& path, // Build the extra res input grid (shifted SW by half (dlon,dlat) // with an added major row column on the NE sides.) int mult = 10; - SimpleMatrix dPts( (xdivs + 1) * mult + 1, (ydivs + 1) * mult + 1 ); + tgMatrix dPts( (xdivs + 1) * mult + 1, (ydivs + 1) * mult + 1 ); for ( int j = 0; j < dPts.rows(); ++j ) { - for ( int i = 0; i < dPts.cols(); ++i ) { - dPts.set(i, j, SGGeod::fromDegM( _min_deg.getLongitudeDeg() - dlon_h - + i * (dlon / (double)mult), - _min_deg.getLatitudeDeg() - dlat_h - + j * (dlat / (double)mult), - -9999 ) - ); + for ( int i = 0; i < dPts.cols(); ++i ) { + dPts.set(i, j, SGGeod::fromDegM( _min_deg.getLongitudeDeg() - dlon_h + i * (dlon / (double)mult), + _min_deg.getLatitudeDeg() - dlat_h + j * (dlat / (double)mult), + -9999 ) ); } } @@ -151,18 +286,18 @@ TGAptSurface::TGAptSurface( const std::string& path, tgClampElevations( dPts, _average_elev_m, max_clamp ); // Build the normal res input grid from the double res version - Pts = new SimpleMatrix(xdivs + 1, ydivs + 1 ); + Pts = new tgMatrix(xdivs + 1, ydivs + 1 ); double ave_divider = (mult+1) * (mult+1); for ( int j = 0; j < Pts->rows(); ++j ) { for ( int i = 0; i < Pts->cols(); ++i ) { - GENAPT_LOG(SG_GENERAL, SG_DEBUG, i << "," << j); + SG_LOG(SG_GENERAL, SG_DEBUG, i << "," << j); double accum = 0.0; double lon_accum = 0.0; double lat_accum = 0.0; for ( int jj = 0; jj <= mult; ++jj ) { for ( int ii = 0; ii <= mult; ++ii ) { double value = dPts.element(mult*i + ii, mult*j + jj).getElevationM(); - GENAPT_LOG( SG_GENERAL, SG_DEBUG, "value = " << value ); + SG_LOG( SG_GENERAL, SG_DEBUG, "value = " << value ); accum += value; lon_accum += dPts.element(mult*i + ii, mult*j + jj).getLongitudeDeg(); lat_accum += dPts.element(mult*i + ii, mult*j + jj).getLatitudeDeg(); @@ -170,7 +305,7 @@ TGAptSurface::TGAptSurface( const std::string& path, } double val_ave = accum / ave_divider; - GENAPT_LOG( SG_GENERAL, SG_DEBUG, " val_ave = " << val_ave ); + SG_LOG( SG_GENERAL, SG_DEBUG, " val_ave = " << val_ave ); Pts->set(i, j, SGGeod::fromDegM( _min_deg.getLongitudeDeg() + i * dlon, _min_deg.getLatitudeDeg() + j * dlat, val_ave ) @@ -180,18 +315,18 @@ TGAptSurface::TGAptSurface( const std::string& path, bool slope_error = true; while ( slope_error ) { - GENAPT_LOG( SG_GENERAL, SG_DEBUG, "start of slope processing pass" ); + SG_LOG( SG_GENERAL, SG_DEBUG, "start of slope processing pass" ); slope_error = false; // Add some "slope" sanity to the resulting surface grid points for ( int j = 0; j < Pts->rows() - 1; ++j ) { for ( int i = 0; i < Pts->cols() - 1; ++i ) { - if ( limit_slope( Pts, i, j, i+1, j, _average_elev_m ) ) { + if ( limit_slope( Pts, i, j, i+1, j, _average_elev_m, slope_max, slope_eps ) ) { slope_error = true; } - if ( limit_slope( Pts, i, j, i, j+1, _average_elev_m ) ) { + if ( limit_slope( Pts, i, j, i, j+1, _average_elev_m, slope_max, slope_eps ) ) { slope_error = true; } - if ( limit_slope( Pts, i, j, i+1, j+1, _average_elev_m ) ) { + if ( limit_slope( Pts, i, j, i+1, j+1, _average_elev_m, slope_max, slope_eps ) ) { slope_error = true; } } @@ -202,23 +337,23 @@ TGAptSurface::TGAptSurface( const std::string& path, double clon = (_min_deg.getLongitudeDeg() + _max_deg.getLongitudeDeg()) / 2.0; double clat = (_min_deg.getLatitudeDeg() + _max_deg.getLatitudeDeg()) / 2.0; area_center = SGGeod::fromDegM( clon, clat, _average_elev_m ); - GENAPT_LOG(SG_GENERAL, SG_DEBUG, "Central offset point = " << area_center); + SG_LOG(SG_GENERAL, SG_DEBUG, "Central offset point = " << area_center); // Create the fitted surface - GENAPT_LOG(SG_GENERAL, SG_DEBUG, "ready to create fitted surface"); + SG_LOG(SG_GENERAL, SG_DEBUG, "ready to create fitted surface"); fit(); - GENAPT_LOG(SG_GENERAL, SG_DEBUG, " fit process successful."); + SG_LOG(SG_GENERAL, SG_DEBUG, " fit process successful."); } -TGAptSurface::~TGAptSurface() { +tgSurface::~tgSurface() { delete Pts; } // Use a linear least squares method to fit a 3d polynomial to the // sampled surface data -void TGAptSurface::fit() { +void tgSurface::fit() { // the fit function is: // f(x,y) = A1*x + A2*x*y + A3*y + @@ -228,7 +363,7 @@ void TGAptSurface::fit() { int nobs = Pts->cols() * Pts->rows(); // number of observations - GENAPT_LOG(SG_GENERAL, SG_DEBUG, "QR triangularisation" ); + SG_LOG(SG_GENERAL, SG_DEBUG, "QR triangularisation" ); // Create an array (matrix) with 16 columns (predictor values) A[n] TNT::Array2D mat(nobs,16); @@ -274,13 +409,12 @@ void TGAptSurface::fit() { // Query the elevation of a point, return -9999 if out of range -double TGAptSurface::query( SGGeod query ) { +double tgSurface::query( SGGeod query ) const { // sanity check if ( !_aptBounds.isInside(query) ) { - GENAPT_LOG(SG_GENERAL, SG_WARN, - "Warning: query out of bounds for fitted surface!"); + SG_LOG(SG_GENERAL, SG_WARN, "Warning: query out of bounds for fitted surface!"); return -9999.0; } @@ -303,4 +437,4 @@ double TGAptSurface::query( SGGeod query ) { result += area_center.getElevationM(); return result; -} +} \ No newline at end of file diff --git a/src/Airports/GenAirports850/apt_surface.hxx b/src/Lib/Polygon/surface.hxx similarity index 54% rename from src/Airports/GenAirports850/apt_surface.hxx rename to src/Lib/Polygon/surface.hxx index ff281cf1..f3a7cead 100644 --- a/src/Airports/GenAirports850/apt_surface.hxx +++ b/src/Lib/Polygon/surface.hxx @@ -21,71 +21,72 @@ // -#ifndef _APT_SURFACE_HXX -#define _APT_SURFACE_HXX - +#ifndef _SURFACE_HXX +#define _SURFACE_HXX #include -#include - #include -#include -#include "debug.hxx" +#include "TNT/tnt_array2d.h" +#include "polygon.hxx" /*** * A dirt simple matrix class for our convenience based on top of SGGeod */ +class tgMatrix { -class SimpleMatrix { - -private: - - int _rows; - int _cols; - tgContour m; +typedef std::vector GeodRow; +typedef std::vector GeodMatrix; public: + inline tgMatrix( unsigned int columns, unsigned int rows ) { + _cols = columns; + _rows = rows; - inline SimpleMatrix( int columns, int rows ) { - _cols = columns; - _rows = rows; - m.Resize( _cols * _rows ); - } - - inline SGGeod element( int col, int row ) { - int index = ( row * _cols ) + col; - if ( col < 0 || col >= _cols ) { - GENAPT_LOG(SG_GENERAL, SG_WARN, "column out of bounds on read (" << col << " >= " << _cols << ")"); - int *p = 0; *p = 1; // force crash - } else if ( row < 0 || row >= _rows ) { - GENAPT_LOG(SG_GENERAL, SG_WARN, "row out of bounds on read (" << row << " >= " << _rows << ")"); - int *p = 0; *p = 1; // force crash + m.resize( rows ); + for (unsigned int i=0; i= _cols ) { - GENAPT_LOG(SG_GENERAL, SG_WARN,"column out of bounds on set (" << col << " >= " << _cols << ")"); - int *p = 0; *p = 1; // force crash - } else if ( row < 0 || row >= _rows ) { - GENAPT_LOG(SG_GENERAL, SG_WARN,"row out of bounds on set (" << row << " >= " << _rows << ")"); - int *p = 0; *p = 1; // force crash + inline SGGeod const& element( unsigned int col, unsigned int row ) const { + //int index = ( row * _cols ) + col; + if ( col < 0 || col >= _cols ) { + SG_LOG(SG_GENERAL, SG_WARN, "column out of bounds on read (" << col << " >= " << _cols << ")"); + int *p = 0; *p = 1; // force crash + } else if ( row < 0 || row >= _rows ) { + SG_LOG(SG_GENERAL, SG_WARN, "row out of bounds on read (" << row << " >= " << _rows << ")"); + int *p = 0; *p = 1; // force crash + } + + return m[row][col]; } - m.SetNode(index, p); - } - inline int cols() const { return _cols; } - inline int rows() const { return _rows; } + inline void set( unsigned int col, unsigned int row, const SGGeod& p ) { + //int index = ( row * _cols ) + col; + if ( col < 0 || col >= _cols ) { + SG_LOG(SG_GENERAL, SG_WARN,"column out of bounds on set (" << col << " >= " << _cols << ")"); + int *p = 0; *p = 1; // force crash + } else if ( row < 0 || row >= _rows ) { + SG_LOG(SG_GENERAL, SG_WARN,"row out of bounds on set (" << row << " >= " << _rows << ")"); + int *p = 0; *p = 1; // force crash + } + m[row][col] = p; + } + + inline int cols() const { return _cols; } + inline int rows() const { return _rows; } + +private: + unsigned int _rows; + unsigned int _cols; + GeodMatrix m; }; - /*** - * Note of explanation. When a TGAptSurface instance is created, you - * must specify a min and max lon/lat containing the entire airport - * area. The class will divide up that area into a reasonably sized + * Note of explanation. When a tgSurface instance is created, you + * must specify a min and max lon/lat containing the entire area. + * The class will divide up that area into a reasonably sized * regular grid. It will then look up the elevation of each point on * the grid from the DEM/Array data. Finally it will fit do a linear * least squares polygonal surface approximation from this grid. Each @@ -94,36 +95,20 @@ public: * provides a) smoothing of noisy terrain data and b) natural rises * and dips in the airport surface. */ - -class TGAptSurface { - -private: - - // The actual nurbs surface approximation for the airport - SimpleMatrix *Pts; - TNT::Array1D surface_coefficients; - - tg::Rectangle _aptBounds; - - SGGeod _min_deg, _max_deg; - - // A central point in the airport area - SGGeod area_center; - - // externally seeded average airport elevation - double _average_elev_m; - +class tgSurface { public: // Constructor, specify min and max coordinates of desired area in // lon/lat degrees, also please specify an "average" airport // elevations in meters. - TGAptSurface( const std::string &path, const string_list& elev_src, - tg::Rectangle aptBounds, double average_elev_m ); + tgSurface( const std::string &path, const string_list& elev_src, + tgRectangle aptBounds, double average_elev_m, + double slope_max, double slope_eps + ); // Destructor - ~TGAptSurface(); + ~tgSurface(); // Use a linear least squares method to fit a 3d polynomial to the // sampled surface data @@ -132,9 +117,24 @@ public: // Query the elevation of a point, return -9999 if out of range. // This routine makes a simplistic assumption that X,Y space is // proportional to u,v space on the nurbs surface which it isn't. - double query( SGGeod query ); + double query( SGGeod query ) const; +private: + // The actual nurbs surface approximation for the airport + tgMatrix* Pts; + TNT::Array1D surface_coefficients; + + tgRectangle _aptBounds; + SGGeod _min_deg, _max_deg; + + // A central point in the airport area + SGGeod area_center; + + // externally seeded average airport elevation + double _average_elev_m; + +// double slope_max = 0.02; +// double slope_eps = 0.00001; }; - -#endif // _APT_SURFACE_HXX +#endif // _SURFACE_HXX diff --git a/src/Utils/poly2ogr/poly2ogr.cxx b/src/Utils/poly2ogr/poly2ogr.cxx index c8713fbb..fec6fb98 100644 --- a/src/Utils/poly2ogr/poly2ogr.cxx +++ b/src/Utils/poly2ogr/poly2ogr.cxx @@ -36,6 +36,7 @@ #endif #include +#include #include #include @@ -45,9 +46,6 @@ #include #include -#include - -#include using std::string; @@ -147,7 +145,7 @@ OGRLayer* get_layer_for_material(const std::string& material) { return layer; } -OGRLinearRing* make_ring_from_fan(const int_list& fan, const std::vector& nodes) { +OGRLinearRing* make_ring_from_fan(const int_list& fan, const std::vector& nodes) { OGRLinearRing* ring = new OGRLinearRing(); int_list::const_iterator vertex = fan.begin(); if (fan[1]==fan[fan.size()-1]) { @@ -156,10 +154,10 @@ OGRLinearRing* make_ring_from_fan(const int_list& fan, const std::vectorsetX(node.x()); - point->setY(node.y()); - point->setZ(node.z()); + const SGGeod& node = nodes[*vertex]; + point->setX(node.getLongitudeDeg()); + point->setY(node.getLatitudeDeg()); + point->setZ(node.getElevationM()); ring->addPoint(point); } @@ -169,25 +167,25 @@ OGRLinearRing* make_ring_from_fan(const int_list& fan, const std::vector& nodes) { +OGRLinearRing* make_ring_from_strip(const int_list& strip, const std::vector& nodes) { OGRLinearRing* ring = new OGRLinearRing(); const size_t vertex_count = strip.size(); unsigned int i; for (i=0;isetX(node.x()); - point->setY(node.y()); - point->setZ(node.z()); + const SGGeod& node = nodes[strip[i]]; + point->setX(node.getLongitudeDeg()); + point->setY(node.getLatitudeDeg()); + point->setZ(node.getElevationM()); ring->addPoint(point); } for (i--;i>0;i-=2) { OGRPoint *point=new OGRPoint(); - const Point3D& node = nodes[strip[i]]; - point->setX(node.x()); - point->setY(node.y()); - point->setZ(node.z()); + const SGGeod& node = nodes[strip[i]]; + point->setX(node.getLongitudeDeg()); + point->setY(node.getLatitudeDeg()); + point->setZ(node.getElevationM()); ring->addPoint(point); } @@ -221,7 +219,7 @@ void make_feature_from_ring(OGRLinearRing* ring, const std::string& material, co make_feature_from_polygon(polygon, material, path); } -void convert_triangles(const std::string& path, const group_list& verts, const string_list& materials, const std::vector& wgs84_nodes) { +void convert_triangles(const std::string& path, const group_list& verts, const string_list& materials, const std::vector& wgs84_nodes) { const size_t groups_count = verts.size(); for (unsigned int i=0;isetX(node.x()); - point->setY(node.y()); - point->setZ(node.z()); + const SGGeod& node = wgs84_nodes[tri_verts[j+k]]; + point->setX(node.getLongitudeDeg()); + point->setY(node.getLatitudeDeg()); + point->setZ(node.getElevationM()); ring->addPoint(point); } @@ -245,7 +243,7 @@ void convert_triangles(const std::string& path, const group_list& verts, const s } } -void convert_triangle_fans(const std::string& path, const group_list& verts, const string_list& materials, const std::vector& wgs84_nodes) { +void convert_triangle_fans(const std::string& path, const group_list& verts, const string_list& materials, const std::vector& wgs84_nodes) { const size_t groups_count = verts.size(); for (unsigned int i=0;i& wgs84_nodes) { +void convert_triangle_strips(const std::string& path, const group_list& verts, const string_list& materials, const std::vector& wgs84_nodes) { const size_t groups_count = verts.size(); for (unsigned int i=0;i& wgs84_nodes = binObject.get_wgs84_nodes(); - std::vector geod_nodes; + std::vector geod_nodes; const size_t node_count = wgs84_nodes.size(); for (unsigned int i=0;i