From 15ae058d4afb8e66e0be6a435e7baf395f7d6017 Mon Sep 17 00:00:00 2001
From: curt <curt>
Date: Sun, 9 Mar 2003 18:11:29 +0000
Subject: [PATCH] First stab at code to build a nurbs surface that approximates
the surface of an airport.
---
src/Airports/GenAirports/apt_surface.cxx | 206 +++++++++++++++++++++++
src/Airports/GenAirports/apt_surface.hxx | 74 ++++++++
2 files changed, 280 insertions(+)
create mode 100644 src/Airports/GenAirports/apt_surface.cxx
create mode 100644 src/Airports/GenAirports/apt_surface.hxx
diff --git a/src/Airports/GenAirports/apt_surface.cxx b/src/Airports/GenAirports/apt_surface.cxx
new file mode 100644
index 00000000..2842712c
--- /dev/null
+++ b/src/Airports/GenAirports/apt_surface.cxx
@@ -0,0 +1,206 @@
+// apt_surface.cxx -- class to manage airport terrain surface
+// approximation and smoothing
+//
+// Written by Curtis Olson, started March 2003.
+//
+// Copyright (C) 2003 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 <nurbs++/nurbsS.hh>
+#include <nurbs++/nurbsSub.hh>
+
+#include <simgear/constants.h>
+#include <simgear/math/sg_types.hxx>
+
+#include <Array/array.hxx>
+
+#include "apt_surface.hxx"
+
+SG_USING_NAMESPACE( PLib );
+
+
+// fix node elevations. Offset is added to the final elevation
+static void calc_elevations( const string &root, Matrix_Point3Df &Pts ) {
+ string_list elev_src;
+ elev_src.clear();
+ elev_src.push_back( "SRTM-1" );
+ elev_src.push_back( "SRTM-3" );
+ elev_src.push_back( "DEM-3" );
+ elev_src.push_back( "DEM-30" );
+
+ 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 ( i = 0; i < Pts.rows(); ++i ) {
+ for ( j = 0; j < Pts.cols(); ++j ) {
+ Point3Df p = Pts(i,j);
+ p.z() = -9999.0;
+ }
+ }
+
+ while ( !done ) {
+ // find first node with -9999 elevation
+ Point3Df first;
+ bool found_one = false;
+ for ( i = 0; i < Pts.rows(); ++i ) {
+ for ( j = 0; j < Pts.cols(); ++j ) {
+ Point3Df p = Pts(i,j);
+ 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() + ".arr";
+ if ( array.open(array_path) ) {
+ found_file = true;
+ cout << "Using array_path = " << array_path << endl;
+ }
+ j++;
+ }
+
+ // this will fill in a zero structure if no array data
+ // found/opened
+ array.parse( b );
+
+ // update all the non-updated elevations that are inside
+ // this array file
+ double elev;
+ done = true;
+ for ( i = 0; i < Pts.rows(); ++i ) {
+ for ( j = 0; j < Pts.cols(); ++j ) {
+ Point3Df p = Pts(i,j);
+ if ( p.z() < -9000.0 ) {
+ done = false;
+ elev = array.interpolate_altitude( p.x() * 3600.0,
+ p.y() * 3600.0 );
+ if ( elev > -9000 ) {
+ p.z() = elev;
+ Pts(i,j) = p;
+ cout << "interpolating for " << p << endl;
+ }
+ }
+ }
+ }
+
+ array.close();
+ } else {
+ done = true;
+ }
+ }
+}
+
+
+// Constructor, specify min and max coordinates of desired area in
+// lon/lat degrees
+TGAptSurface::TGAptSurface( const string& path,
+ Point3D _min_deg, Point3D _max_deg )
+{
+ // Calculate desired size of grid
+ min_deg = _min_deg;
+ max_deg = _max_deg;
+
+ // The following size calculations are for the purpose of
+ // determining grid divisions so it's not important that they be
+ // *exact*, just ball park.
+ double y_deg = max_deg.lat() - min_deg.lat();
+ double y_rad = y_deg * SG_DEGREES_TO_RADIANS;
+ double y_nm = y_rad * SG_RAD_TO_NM;
+ double y_m = y_nm * SG_NM_TO_METER;
+
+ double xfact = cos( min_deg.lat() * SG_DEGREES_TO_RADIANS );
+ double x_deg = max_deg.lon() - min_deg.lon();
+ double x_rad = x_deg * SG_DEGREES_TO_RADIANS;
+ double x_nm = x_rad * SG_RAD_TO_NM * xfact;
+ double x_m = x_nm * SG_NM_TO_METER;
+
+ cout << "Area size = " << x_m << " x " << y_m << " (m)" << endl;
+
+ int xdivs = (int)(x_m / 400.0) + 1;
+ int ydivs = (int)(y_m / 400.0) + 1;
+
+ // Build the input grid
+
+ cout << " M(" << xdivs << "," << ydivs << ")" << endl;
+ Matrix_Point3Df Pts(xdivs, ydivs);
+
+ double dlon = x_deg / xdivs;
+ double dlat = y_deg / ydivs;
+
+ for ( int i = 0; i < xdivs; ++i ) {
+ for ( int j = 0; j < ydivs; ++j ) {
+ Pts(i,j) = Point3Df( min_deg.lon() + i * dlon,
+ min_deg.lat() + j * dlat,
+ -9999 );
+ }
+ }
+
+ // Determine elevation of the grid points
+ calc_elevations( path, Pts );
+
+ // Create the nurbs surface
+
+ apt_surf = new PlNurbsSurfacef;
+ apt_surf->globalInterp( Pts, 3, 3);
+}
+
+
+TGAptSurface::~TGAptSurface() {
+ delete apt_surf;
+}
+
+
+// Query the elevation of a point, return -9999 if out of range
+double TGAptSurface::query( double lon_deg, double lat_deg ) {
+ // sanity check
+ if ( lon_deg < min_deg.lon() || lon_deg > max_deg.lon() ||
+ lat_deg < min_deg.lat() || lat_deg > max_deg.lat() )
+ {
+ cout << "Warning: query out of bounds for NURBS surface!" << endl;
+ return -9999.0;
+ }
+
+ // convert lon,lat to nurbs space
+ double x = (lon_deg - min_deg.lon()) / (max_deg.lon() - min_deg.lon());
+ double y = (lat_deg - min_deg.lat()) / (max_deg.lat() - min_deg.lat());
+
+ cout << " querying for " << x << ", " << y << " = ";
+ Point3Df p = apt_surf->pointAt( x, y );
+ cout << p.z() << endl;
+
+ return p.z();
+}
diff --git a/src/Airports/GenAirports/apt_surface.hxx b/src/Airports/GenAirports/apt_surface.hxx
new file mode 100644
index 00000000..993305bc
--- /dev/null
+++ b/src/Airports/GenAirports/apt_surface.hxx
@@ -0,0 +1,74 @@
+// apt_surface.hxx -- class to manage airport terrain surface
+// approximation and smoothing
+//
+// Written by Curtis Olson, started March 2003.
+//
+// Copyright (C) 2003 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$
+//
+
+
+#ifndef _APT_SURFACE_HXX
+#define _APT_SURFACE_HXX
+
+
+#include <string>
+
+#include <nurbs++/nurbsS.hh>
+
+#include <simgear/math/point3d.hxx>
+
+
+/***
+ * 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
+ * regular grid. It will then look up the elevation of each point on
+ * the grid from the DEM/Array data. Finally it will build a nurbs
+ * surface from this grid. Each vertex of the actual airport model is
+ * drapped over this nurbs surface rather than over the underlying
+ * terrain data. This provides a) smoothing of noisy terrain data, b)
+ * natural rises and dips in the airport surface, c) chance to impress
+ * colleages by using something or other nurbs surfaces -- or whatever
+ * they are called. :-)
+ */
+
+class TGAptSurface {
+
+private:
+
+ // The actual nurbs surface approximation for the airport
+ PlNurbsSurfacef *apt_surf;
+
+ Point3D min_deg, max_deg;
+
+public:
+
+ // Constructor, specify min and max coordinates of desired area in
+ // lon/lat degrees
+ TGAptSurface( const string &path, Point3D _min_deg, Point3D _max_deg );
+
+ // Destructor
+ ~TGAptSurface();
+
+ // Query the elevation of a point, return -9999 if out of range
+ double query( double lon_deg, double lat_deg );
+};
+
+
+#endif // _APT_SURFACE_HXX