From fd0a6e29fcff80914cec2c751da71fd3fb2e21cd Mon Sep 17 00:00:00 2001 From: curt Date: Wed, 21 Apr 2004 17:03:59 +0000 Subject: [PATCH] Ok, this change seems to help a lot. There are quite a few difficult cases where an airport is located on top of a hill, or in a bowl, or has a significant elevation change near by. I take the average elevation of the area and clamp the outlyers. However these difficult cases "bias" the average elevation because the airport surface may include much of the surrounding area. This change to the code computes the airport elevation *only* based on the actual airport geometry node and ignores all the surrounding nonsense that might exist. This doesn't make things perfect, but is a *big* step forward for airports in areas with significant elevation change nearby. --- src/Airports/GenAirports/Makefile.am | 1 + src/Airports/GenAirports/apt_surface.cxx | 33 ++- src/Airports/GenAirports/apt_surface.hxx | 8 +- src/Airports/GenAirports/build.cxx | 26 +- src/Airports/GenAirports/elevations.cxx | 297 +++++++++++++++++++++++ src/Airports/GenAirports/elevations.hxx | 54 +++++ src/Airports/GenAirports/main.cxx | 2 +- 7 files changed, 400 insertions(+), 21 deletions(-) create mode 100644 src/Airports/GenAirports/elevations.cxx create mode 100644 src/Airports/GenAirports/elevations.hxx diff --git a/src/Airports/GenAirports/Makefile.am b/src/Airports/GenAirports/Makefile.am index c0f7cdc1..3930364a 100644 --- a/src/Airports/GenAirports/Makefile.am +++ b/src/Airports/GenAirports/Makefile.am @@ -29,6 +29,7 @@ genapts_SOURCES = \ apt_surface.hxx apt_surface.cxx \ build.cxx build.hxx \ convex_hull.cxx convex_hull.hxx \ + elevations.cxx elevations.hxx \ global.hxx \ lights.hxx lights.cxx \ main.cxx \ diff --git a/src/Airports/GenAirports/apt_surface.cxx b/src/Airports/GenAirports/apt_surface.cxx index 269ec7f5..316d0d69 100644 --- a/src/Airports/GenAirports/apt_surface.cxx +++ b/src/Airports/GenAirports/apt_surface.cxx @@ -33,12 +33,14 @@ #include +#include "elevations.hxx" #include "global.hxx" #include "apt_surface.hxx" 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, @@ -164,17 +166,20 @@ static double calc_elevations( const string &root, const string_list elev_src, return average; } +#endif // Constructor, specify min and max coordinates of desired area in // lon/lat degrees TGAptSurface::TGAptSurface( const string& path, const string_list& elev_src, - Point3D _min_deg, Point3D _max_deg ) + Point3D _min_deg, Point3D _max_deg, + double _average_elev_m ) { // Calculate desired size of grid min_deg = _min_deg; max_deg = _max_deg; + average_elev_m = _average_elev_m; // The following size calculations are for the purpose of // determining grid divisions so it's not important that they be @@ -220,8 +225,12 @@ TGAptSurface::TGAptSurface( const string& path, } } - // Determine elevation of the grid points - double average = calc_elevations( path, elev_src, dPts ); + // Lookup the elevations of all the grid points + tgCalcElevations( path, elev_src, dPts ); + + // Clamp the elevations against the externally provided average + // elevation. + tgClampElevations( dPts, average_elev_m, max_clamp ); // Build the normal res input grid from the double res version Matrix_Point3Dd Pts(ydivs + 1, xdivs + 1); @@ -259,7 +268,7 @@ TGAptSurface::TGAptSurface( const string& path, bool slope_error = true; while ( slope_error ) { - // cout << "start of slope processing pass" << endl; + 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 i = 0; i < Pts.cols() - 1; ++i ) { @@ -280,11 +289,11 @@ TGAptSurface::TGAptSurface( const string& path, slope_error = true; - // cout << " (a) detected slope of " << slope - // << " dist = " << dist << endl; + SG_LOG( SG_GENERAL, SG_DEBUG, " (a) detected slope of " + << slope << " dist = " << dist ); - double e1 = average - p1.z(); - double e2 = average - p2.z(); + 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; @@ -320,11 +329,11 @@ TGAptSurface::TGAptSurface( const string& path, slope_error = true; - // cout << " (b) detected slope of " << slope - // << " dist = " << dist << endl; + SG_LOG( SG_GENERAL, SG_DEBUG, " (b) detected slope of " + << slope << " dist = " << dist ); - double e1 = average - p1.z(); - double e2 = average - p2.z(); + double e1 = average_elev_m - p1.z(); + double e2 = average_elev_m - p2.z(); if ( fabs(e1) > fabs(e2) ) { if ( slope > 0 ) { diff --git a/src/Airports/GenAirports/apt_surface.hxx b/src/Airports/GenAirports/apt_surface.hxx index dffd66fb..046b961c 100644 --- a/src/Airports/GenAirports/apt_surface.hxx +++ b/src/Airports/GenAirports/apt_surface.hxx @@ -57,12 +57,16 @@ private: Point3D min_deg, max_deg; + // externally seeded average airport elevation + double average_elev_m; + public: // Constructor, specify min and max coordinates of desired area in - // lon/lat degrees + // lon/lat degrees, also please specify an "average" airport + // elevations in meters. TGAptSurface( const string &path, const string_list& elev_src, - Point3D _min_deg, Point3D _max_deg ); + Point3D _min_deg, Point3D _max_deg, double _average_elev_m ); // Destructor ~TGAptSurface(); diff --git a/src/Airports/GenAirports/build.cxx b/src/Airports/GenAirports/build.cxx index 31880929..774fe212 100644 --- a/src/Airports/GenAirports/build.cxx +++ b/src/Airports/GenAirports/build.cxx @@ -58,6 +58,7 @@ #include "apt_surface.hxx" #include "convex_hull.hxx" +#include "elevations.hxx" #include "lights.hxx" #include "point2d.hxx" #include "poly_extra.hxx" @@ -961,6 +962,18 @@ void build_airport( string airport_id, float alt_m, tris_tc.push_back( base_tc ); } + // Now that we have assembled all the airport geometry nodes into + // a list, calculate an "average" airport elevation based on all + // the actual airport node points. This is more useful than + // calculating an average over the entire airport surface because + // it avoids biases introduced from the surrounding area if the + // airport is located in a bowl or on a hill. + + double average = tgAverageElevation( root, elev_src, + nodes.get_node_list() ); + + // Now build the airport nurbs surface ... + // calculation min/max coordinates of airport area Point3D min_deg(9999.0, 9999.0, 0), max_deg(-9999.0, -9999.0, 0); for ( j = 0; j < (int)nodes.get_node_list().size(); ++j ) { @@ -1005,16 +1018,17 @@ void build_airport( string airport_id, float alt_m, // 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(); - min_deg.setlon( min_deg.lon() - 0.2 * dlon ); - max_deg.setlon( max_deg.lon() + 0.2 * dlon ); - min_deg.setlat( min_deg.lat() - 0.2 * dlat ); - max_deg.setlat( max_deg.lat() + 0.2 * dlat ); + min_deg.setlon( min_deg.lon() - 0.3 * dlon ); + max_deg.setlon( max_deg.lon() + 0.3 * dlon ); + min_deg.setlat( min_deg.lat() - 0.3 * dlat ); + max_deg.setlat( max_deg.lat() + 0.3 * dlat ); - TGAptSurface apt_surf( root, elev_src, min_deg, max_deg ); + TGAptSurface apt_surf( root, elev_src, min_deg, max_deg, average ); SG_LOG(SG_GENERAL, SG_DEBUG, "Surface created"); // calculate node elevations - point_list geod_nodes = calc_elevations( apt_surf, nodes.get_node_list(), + point_list geod_nodes = calc_elevations( apt_surf, + nodes.get_node_list(), 0.0 ); divided_base = calc_elevations( apt_surf, divided_base, 0.0 ); SG_LOG(SG_GENERAL, SG_DEBUG, "DIVIDED"); diff --git a/src/Airports/GenAirports/elevations.cxx b/src/Airports/GenAirports/elevations.cxx new file mode 100644 index 00000000..36990dfa --- /dev/null +++ b/src/Airports/GenAirports/elevations.cxx @@ -0,0 +1,297 @@ +// elevations.cxx -- routines to help calculate DEM elevations for a +// set of points +// +// Written by Curtis Olson, started April 2004. +// +// Copyright (C) 2004 Curtis L. Olson - curt@flightgear.org +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ +// + + +#include +#include + +#include +#include +#include +#include + +#include + +#include "global.hxx" +#include "apt_surface.hxx" + +SG_USING_NAMESPACE( PLib ); + + +// lookup node elevations for each point in the point_list. Returns +// average of all points. Doesn't modify the original list. + +double tgAverageElevation( const string &root, const string_list elev_src, + const point_list points_source ) +{ + bool done = false; + unsigned int i; + TGArray array; + + // make a copy so our routine is non-destructive. + point_list points = points_source; + + // just bail if no work to do + if ( points.size() == 0 ) { + return 0.0; + } + + // set all elevations to -9999 + for ( i = 0; i < points.size(); ++i ) { + points[i].setz( -9999.0 ); + } + + while ( !done ) { + // find first node with -9999 elevation + Point3D first; + bool found_one = false; + for ( i = 0; i < points.size(); ++i ) { + if ( points[i].z() < -9000.0 && !found_one ) { + first = points[i]; + found_one = true; + } + } + + if ( found_one ) { + SGBucket b( first.x(), first.y() ); + string base = b.gen_base_path(); + + // try the various elevation sources + i = 0; + bool found_file = false; + while ( !found_file && i < elev_src.size() ) { + string array_path = root + "/" + elev_src[i] + "/" + base + + "/" + b.gen_index_str(); + if ( array.open(array_path) ) { + found_file = true; + SG_LOG( SG_GENERAL, SG_DEBUG, "Using array_path = " + << array_path ); + } + i++; + } + + // 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 < points.size(); ++i ) { + if ( points[i].z() < -9000.0 ) { + done = false; + elev = array.altitude_from_grid( points[i].lon() * 3600.0, + points[i].lat() * 3600.0 ); + if ( elev > -9000 ) { + points[i].setz( elev ); + // cout << "interpolating for " << p << endl; + // cout << p.x() << " " << p.y() << " " << p.z() + // << endl; + } + } + } + + array.close(); + } else { + done = true; + } + } + + // now find the average height of the queried points + double total = 0.0; + int count = 0; + for ( i = 0; i < points.size(); ++i ) { + total += points[i].z(); + count++; + } + double average = total / (double) count; + SG_LOG(SG_GENERAL, SG_DEBUG, "Average surface height of point list = " + << average); + + return average; +} + + +// 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 ) { + 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; + } + + // 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; + } + } + + 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 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; +} + + +// clamp all elevations to the specified range + +void tgClampElevations( Matrix_Point3Dd &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 ( i = 0; i < Pts.cols(); ++i ) { + for ( j = 0; j < Pts.rows(); ++j ) { + Point3Dd p = Pts(j,i); + if ( p.z() < center_m - max_clamp_m ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.z() + << " to " << center_m - max_clamp_m ); + p.z() = center_m - max_clamp_m; + Pts(j,i) = p; + } + if ( p.z() > center_m + max_clamp_m ) { + SG_LOG(SG_GENERAL, SG_DEBUG, " clamping " << p.z() + << " to " << center_m + max_clamp_m ); + p.z() = center_m + max_clamp_m; + Pts(j,i) = p; + } + } + } +} diff --git a/src/Airports/GenAirports/elevations.hxx b/src/Airports/GenAirports/elevations.hxx new file mode 100644 index 00000000..3b85d124 --- /dev/null +++ b/src/Airports/GenAirports/elevations.hxx @@ -0,0 +1,54 @@ +// elevations.hxx -- routines to help calculate DEM elevations for a +// set of points +// +// Written by Curtis Olson, started April 2004. +// +// Copyright (C) 2004 Curtis L. Olson - curt@flightgear.org +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. +// +// $Id$ +// + + +#include +#include + +#include +#include +#include +#include + +#include + +#include "global.hxx" +#include "apt_surface.hxx" + +SG_USING_NAMESPACE( PLib ); + + +// lookup node elevations for each point in the point_list. Returns +// average of all points. Doesn't modify the original list. +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 ); + +// clamp all elevations to the specified range +void tgClampElevations( Matrix_Point3Dd &Pts, + double center_m, double max_clamp_m ); diff --git a/src/Airports/GenAirports/main.cxx b/src/Airports/GenAirports/main.cxx index af92784f..bd0e1110 100644 --- a/src/Airports/GenAirports/main.cxx +++ b/src/Airports/GenAirports/main.cxx @@ -118,7 +118,7 @@ int main( int argc, char **argv ) { string_list elev_src; elev_src.clear(); - sglog().setLogLevels( SG_GENERAL, SG_DEBUG ); + sglog().setLogLevels( SG_GENERAL, SG_INFO ); // parse arguments string work_dir = "";