1
0
Fork 0

First stab at code to build a nurbs surface that approximates the surface

of an airport.
This commit is contained in:
curt 2003-03-09 18:11:29 +00:00
parent b628d7cdd8
commit 15ae058d4a
2 changed files with 280 additions and 0 deletions

View file

@ -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();
}

View file

@ -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