diff --git a/src/Airports/GenAirports/apt_surface.cxx b/src/Airports/GenAirports/apt_surface.cxx index 9e66a5f4..63ac256f 100644 --- a/src/Airports/GenAirports/apt_surface.cxx +++ b/src/Airports/GenAirports/apt_surface.cxx @@ -40,133 +40,58 @@ SG_USING_NAMESPACE( PLib ); -#if 0 -// fix node elevations. Offset is added to the final elevation, -// returns average of all points. -static double calc_elevations( const string &root, const string_list elev_src, - Matrix_Point3Dd &Pts ) { - bool done = false; - int i, j; - TGArray array; +static bool limit_slope( Matrix_Point3Dd &Pts, int i1, int j1, int i2, int j2, + double average_elev_m ) +{ + bool slope_error = false; - // just bail if no work to do - if ( Pts.rows() == 0 || Pts.cols() == 0 ) { - return 0.0; - } + Point3Dd p1, p2; + p1 = Pts(j1,i1); + p2 = Pts(j2,i2); - // set all elevations to -9999 - for ( i = 0; i < Pts.cols(); ++i ) { - for ( j = 0; j < Pts.rows(); ++j ) { - Point3Dd p = Pts(j,i); - p.z() = -9999.0; + double az1, az2, dist; + double slope; + + geo_inverse_wgs_84( 0, p1.y(), p1.x(), p2.y(), p2.x(), + &az1, &az2, &dist ); + slope = (p2.z() - p1.z()) / dist; + + if ( fabs(slope) > (slope_max + slope_eps) ) { + // need to throttle slope, let's move the point + // furthest away from average towards the center. + + slope_error = true; + + SG_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of " + << slope << " dist = " << dist ); + + double e1 = fabs(average_elev_m - p1.z()); + double e2 = fabs(average_elev_m - p2.z()); + // cout << " z1 = " << p1.z() << " z2 = " << p2.z() << endl; + // cout << " e1 = " << e1 << " e2 = " << e2 << endl; + + if ( e1 > e2 ) { + // cout << " p1 error larger" << endl; + if ( slope > 0 ) { + p1.z() = p2.z() - (dist * slope_max); + } else { + p1.z() = p2.z() + (dist * slope_max); + } + Pts(j1,i1) = p1; + } else { + // cout << " p2 error larger" << endl; + if ( slope > 0 ) { + p2.z() = p1.z() + (dist * slope_max); + } else { + p2.z() = p1.z() - (dist * slope_max); + } + Pts(j2,i2) = p2; } + // cout << " z1 = " << p1.z() << " z2 = " << p2.z() << endl; } - while ( !done ) { - // find first node with -9999 elevation - Point3Dd first; - bool found_one = false; - for ( i = 0; i < Pts.cols(); ++i ) { - for ( j = 0; j < Pts.rows(); ++j ) { - Point3Dd p = Pts(j,i); - if ( p.z() < -9000.0 && !found_one ) { - first = p; - found_one = true; - } - } - } - - if ( found_one ) { - SGBucket b( first.x(), first.y() ); - 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() ) { - 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 ( i = 0; i < Pts.cols(); ++i ) { - for ( j = 0; j < Pts.rows(); ++j ) { - Point3Dd p = Pts(j,i); - if ( p.z() < -9000.0 ) { - done = false; - elev = array.altitude_from_grid( p.x() * 3600.0, - p.y() * 3600.0 ); - if ( elev > -9000 ) { - p.z() = elev; - Pts(j,i) = p; - // cout << "interpolating for " << p << endl; - // cout << p.x() << " " << p.y() << " " << p.z() - // << endl; - } - } - } - } - - array.close(); - } else { - done = true; - } - } - - // do some post processing for sanity's sake - - // find the average height of the queried points - double total = 0.0; - int count = 0; - for ( i = 0; i < Pts.cols(); ++i ) { - for ( j = 0; j < Pts.rows(); ++j ) { - Point3Dd p = Pts(j,i); - total += p.z(); - count++; - } - } - double average = total / (double) count; - SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height = " << average); - - // now go through the elevations and clamp them all to within - // +/-50m (164') of the average. - for ( i = 0; i < Pts.cols(); ++i ) { - for ( j = 0; j < Pts.rows(); ++j ) { - Point3Dd p = Pts(j,i); - if ( p.z() < average - max_clamp ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.z() - << " to " << average - max_clamp ); - p.z() = average - max_clamp; - Pts(j,i) = p; - } - if ( p.z() > average + max_clamp ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.z() - << " to " << average + max_clamp ); - p.z() = average + max_clamp; - Pts(j,i) = p; - } - } - } - - return average; + return slope_error; } -#endif // Constructor, specify min and max coordinates of desired area in @@ -226,7 +151,7 @@ TGAptSurface::TGAptSurface( const string& path, } // Lookup the elevations of all the grid points - tgCalcElevations( path, elev_src, dPts ); + tgCalcElevations( path, elev_src, dPts, 0.0 ); // Clamp the elevations against the externally provided average // elevation. @@ -273,83 +198,14 @@ TGAptSurface::TGAptSurface( const string& path, // Add some "slope" sanity to the resulting surface grid points for ( int i = 0; i < Pts.cols() - 1; ++i ) { for ( int j = 0; j < Pts.rows() - 1; ++j ) { - Point3Dd p1, p2; - double az1, az2, dist; - double slope; - - p1 = Pts(j,i); - p2 = Pts(j,i+1); - geo_inverse_wgs_84( 0, p1.y(), p1.x(), p2.y(), p2.x(), - &az1, &az2, &dist ); - slope = (p2.z() - p1.z()) / dist; - - if ( fabs(slope) > (slope_max + slope_eps) ) { - // need to throttle slope, let's move the point - // furthest away from average towards the center. - + if ( limit_slope( Pts, i, j, i+1, j, average_elev_m ) ) { slope_error = true; - - SG_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of " - << slope << " dist = " << dist ); - - double e1 = average_elev_m - p1.z(); - double e2 = average_elev_m - p2.z(); - // cout << " z1 = " << p1.z() << " z2 = " << p2.z() << endl; - // cout << " e1 = " << e1 << " e2 = " << e2 << endl; - - if ( fabs(e1) > fabs(e2) ) { - // cout << " p1 error larger" << endl; - if ( slope > 0 ) { - p1.z() = p2.z() - (dist * slope_max); - } else { - p1.z() = p2.z() + (dist * slope_max); - } - Pts(j,i) = p1; - } else { - // cout << " p2 error larger" << endl; - if ( slope > 0 ) { - p2.z() = p1.z() + (dist * slope_max); - } else { - p2.z() = p1.z() - (dist * slope_max); - } - Pts(j,i+1) = p2; - } - // cout << " z1 = " << p1.z() << " z2 = " << p2.z() << endl; } - - p1 = Pts(j,i); - p2 = Pts(j+1,i); - geo_inverse_wgs_84( 0, p1.y(), p1.x(), p2.y(), p2.x(), - &az1, &az2, &dist ); - slope = ( p2.z() - p1.z() ) / dist; - - if ( fabs(slope) > (slope_max+slope_eps) ) { - // need to throttle slope, let's move the point - // furthest away from average towards the center. - + if ( limit_slope( Pts, i, j, i, j+1, average_elev_m ) ) { + slope_error = true; + } + if ( limit_slope( Pts, i, j, i+1, j+1, average_elev_m ) ) { slope_error = true; - - SG_LOG( SG_GENERAL, SG_DEBUG, " (b) detected slope of " - << slope << " dist = " << dist ); - - double e1 = average_elev_m - p1.z(); - double e2 = average_elev_m - p2.z(); - - if ( fabs(e1) > fabs(e2) ) { - if ( slope > 0 ) { - p1.z() = p2.z() - (dist * slope_max); - } else { - p1.z() = p2.z() + (dist * slope_max); - } - Pts(j,i) = p1; - } else { - if ( slope > 0 ) { - p2.z() = p1.z() + (dist * slope_max); - } else { - p2.z() = p1.z() - (dist * slope_max); - } - Pts(j+1,i) = p2; - } } } } diff --git a/src/Airports/GenAirports/elevations.cxx b/src/Airports/GenAirports/elevations.cxx index 36990dfa..7ecd497a 100644 --- a/src/Airports/GenAirports/elevations.cxx +++ b/src/Airports/GenAirports/elevations.cxx @@ -141,15 +141,15 @@ double tgAverageElevation( const string &root, const string_list elev_src, // lookup node elevations for each point in the specified nurbs++ // matrix. Returns average of all points. -double tgCalcElevations( const string &root, const string_list elev_src, - Matrix_Point3Dd &Pts ) { +void tgCalcElevations( const string &root, const string_list elev_src, + Matrix_Point3Dd &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 0.0; + return; } // set all elevations to -9999 @@ -240,31 +240,9 @@ double tgCalcElevations( const string &root, const string_list elev_src, count++; } } - double average = total / (double) count; + double grid_average = total / (double) count; SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of nurbs matrix = " - << average); - - // now go through the elevations and clamp them all to within - // +/-50m (164') of the average. - for ( i = 0; i < Pts.cols(); ++i ) { - for ( j = 0; j < Pts.rows(); ++j ) { - Point3Dd p = Pts(j,i); - if ( p.z() < average - max_clamp ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.z() - << " to " << average - max_clamp ); - p.z() = average - max_clamp; - Pts(j,i) = p; - } - if ( p.z() > average + max_clamp ) { - SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.z() - << " to " << average + max_clamp ); - p.z() = average + max_clamp; - Pts(j,i) = p; - } - } - } - - return average; + << grid_average); } diff --git a/src/Airports/GenAirports/elevations.hxx b/src/Airports/GenAirports/elevations.hxx index 3b85d124..85546190 100644 --- a/src/Airports/GenAirports/elevations.hxx +++ b/src/Airports/GenAirports/elevations.hxx @@ -45,9 +45,9 @@ double tgAverageElevation( const string &root, const string_list elev_src, const point_list points_source ); // lookup node elevations for each point in the specified nurbs++ -// matrix. Returns average of all points. -double tgCalcElevations( const string &root, const string_list elev_src, - Matrix_Point3Dd &Pts ); +// matrix. +void tgCalcElevations( const string &root, const string_list elev_src, + Matrix_Point3Dd &Pts, double average ); // clamp all elevations to the specified range void tgClampElevations( Matrix_Point3Dd &Pts, diff --git a/src/Airports/GenAirports/global.hxx b/src/Airports/GenAirports/global.hxx index ac5b70d1..f431e13c 100644 --- a/src/Airports/GenAirports/global.hxx +++ b/src/Airports/GenAirports/global.hxx @@ -36,7 +36,7 @@ const double coarse_grid = 700.0; const double max_clamp = 100.0; // maximum slope (rise/run) allowed on an airport surface -const double slope_max = 0.01; +const double slope_max = 0.009; const double slope_eps = 0.00001; #endif // _GEN_AIRPORT_GLOBAL_HXX