fly/terragear
1
0
Fork 0
Code Activity

- Point3D removed

Browse source 
- Moved AptSurface to the Polygon library - and made it more generic for
  possible future use in construct
This commit is contained in:
Peter Sadrozinski 2012-12-15 21:14:39 -05:00
parent 91b28879d4
commit b1012cd5e9
25 changed files with 506 additions and 1237 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/Airports/GenAirports850/CMakeLists.txt
View file

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

21
src/Airports/GenAirports850/airport.cxx
View file

@ -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<SGGeod> calc_elevations( TGAptSurface &surf, const std::vector<SGGeod>& geod_nodes, double offset )
static std::vector<SGGeod> calc_elevations( const tgSurface& surf, const std::vector<SGGeod>& geod_nodes, double offset )
{
std::vector<SGGeod> result = geod_nodes;
for ( unsigned int i = 0; i < result.size(); ++i ) {
@ -225,9 +225,7 @@ static std::vector<SGGeod> 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

1
src/Airports/GenAirports850/debug.hxx
View file

@ -4,6 +4,7 @@
#include <vector>
#include <simgear/threads/SGThread.hxx>
#include <simgear/debug/logstream.hxx>
#ifndef __DEBUG_HXX__
#define __DEBUG_HXX__

138
src/Airports/GenAirports850/elevations.cxx
View file

@ -31,7 +31,7 @@
#include <Array/array.hxx>
#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);
}
}
}
}
}

7
src/Airports/GenAirports850/elevations.hxx
View file

@ -20,8 +20,7 @@
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
#include "apt_surface.hxx"
#include <Polygon/surface.hxx>
// 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 );

12
src/Airports/GenAirports850/global.hxx
View file

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

5
src/Airports/GenAirports850/main.cxx
View file

@ -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 == "" )

2
src/Airports/GenAirports850/scheduler.cxx
View file

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

2
src/Airports/GenAirports850/scheduler.hxx
View file

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

8
src/BuildTiles/Main/tgconstruct.hxx
View file

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

8
src/BuildTiles/Main/tgconstruct_cleanup.cxx
View file

@ -25,19 +25,15 @@
# include <config.h>
#endif
//#include <iostream>
//#include <simgear/compiler.h>
#include <simgear/compiler.h>
#include <simgear/debug/logstream.hxx>
#include "tgconstruct.hxx"
using std::string;
void TGConstruct::FixTJunctions( void ) {
int before, after;
std::vector<SGGeod> points;
tg::Rectangle bb;
tgRectangle bb;
// traverse each poly, and add intermediate nodes
for ( unsigned int i = 0; i < TG_MAX_AREA_TYPES; ++i ) {

8
src/BuildTiles/Main/tgconstruct_elevation.cxx
View file

@ -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<SGGeod> const& corner_list = array.get_corner_list();
for (unsigned int i=0; i<corner_list.size(); i++) {
nodes.unique_add( corner_list[i].toSGGeod() );
nodes.unique_add( corner_list[i] );
}
point_list fit_list = array.get_fitted_list();
std::vector<SGGeod> const& fit_list = array.get_fitted_list();
for (unsigned int i=0; i<fit_list.size(); i++) {
nodes.unique_add( fit_list[i].toSGGeod() );
nodes.unique_add( fit_list[i] );
}
}
}

75
src/BuildTiles/Main/tgconstruct_landclass.cxx
View file

@ -26,6 +26,7 @@
#endif
#include <simgear/debug/logstream.hxx>
#include <simgear/math/SGMath.hxx>
#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)
{

3
src/BuildTiles/Main/tgconstruct_tesselate.cxx
View file

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

12
src/Lib/Array/array.cxx
View file

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

11
src/Lib/Array/array.hxx
View file

@ -28,8 +28,6 @@
#include <simgear/math/sg_types.hxx>
#include <simgear/misc/sgstream.hxx>
#include <Lib/Polygon/point3d.hxx>
class TGArray {
private:
@ -51,8 +49,8 @@ private:
short *in_data;
// output nodes
point_list corner_list;
point_list fitted_list;
std::vector<SGGeod> corner_list;
std::vector<SGGeod> 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<SGGeod> const& get_corner_list() const { return corner_list; }
inline std::vector<SGGeod> 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

3
src/Lib/Polygon/CMakeLists.txt
View file

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

569
src/Lib/Polygon/point3d.hxx
View file

@ -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 <simgear/compiler.h>
#include <stdio.h>
#include <cassert>
#include <cmath>
#include <istream>
#include <ostream>
#include <vector>
#include <simgear/math/SGMath.hxx>
#include <simgear/io/lowlevel.hxx>
//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<double>& cart);
static Point3D fromSGVec3(const SGVec3<float>& cart);
static Point3D fromSGVec2(const SGVec2<double>& cart);
static Point3D fromSGVec2(const SGVec2<float>& 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<double>& cart)
{
Point3D pt;
pt.setx(cart.x());
pt.sety(cart.y());
pt.setz(cart.z());
return pt;
}
inline Point3D Point3D::fromSGVec3(const SGVec3<float>& cart)
{
Point3D pt;
pt.setx(cart.x());
pt.sety(cart.y());
pt.setz(cart.z());
return pt;
}
inline Point3D Point3D::fromSGVec2(const SGVec2<double>& cart)
{
Point3D pt;
pt.setx(cart.x());
pt.sety(cart.y());
pt.setz(0);
return pt;
}
inline Point3D Point3D::fromSGVec2(const SGVec2<float>& 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<double>::round( n[PX]/grid );
n[PY] = grid * SGMisc<double>::round( n[PY]/grid );
n[PZ] = grid * SGMisc<double>::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<double> cart;
SGGeodesy::SGGeodToCart(SGGeod::fromRadM(geod.lon(), geod.lat(), geod.elev()), cart);
return Point3D::fromSGVec3(cart);
}
#endif // _POINT3D_HXX

178
src/Lib/Polygon/polygon.cxx
View file

@ -30,14 +30,13 @@
#include <simgear/threads/SGThread.hxx>
#include <simgear/threads/SGGuard.hxx>
#include <simgear/math/sg_geodesy.hxx>
#include <simgear/io/lowlevel.hxx>
#include <simgear/misc/texcoord.hxx>
#include <simgear/structure/exception.hxx>
#include <simgear/debug/logstream.hxx>
#include <simgear/bucket/newbucket.hxx>
#include <simgear/misc/sg_path.hxx>
#include <Polygon/point3d.hxx>
#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<SGGeod>& 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<SGGeod>& 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

5
src/Lib/Polygon/polygon.hxx
View file

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

75
src/Lib/Polygon/rectangle.cxx
View file

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

176
src/Lib/Polygon/rectangle.hxx
View file

@ -14,9 +14,6 @@
#include <simgear/compiler.h>
#include <simgear/math/SGMath.hxx>
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

226
src/Airports/GenAirports850/apt_surface.cxx → src/Lib/Polygon/surface.cxx
View file

@ -24,21 +24,25 @@
# include <config.h>
#endif
#include <Lib/Polygon/TNT/jama_qr.h>
#include <simgear/compiler.h>
#include <simgear/constants.h>
#include <simgear/math/SGMath.hxx>
#include <simgear/debug/logstream.hxx>
#include "elevations.hxx"
#include "global.hxx"
#include "apt_surface.hxx"
#include "debug.hxx"
#include <Array/array.hxx>
#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<double> 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;
}
}

142
src/Airports/GenAirports850/apt_surface.hxx → src/Lib/Polygon/surface.hxx
View file

@ -21,71 +21,72 @@
//
#ifndef _APT_SURFACE_HXX
#define _APT_SURFACE_HXX
#ifndef _SURFACE_HXX
#define _SURFACE_HXX
#include <string>
#include <Lib/Polygon/TNT/tnt_array2d.h>
#include <simgear/debug/logstream.hxx>
#include <Lib/Polygon/polygon.hxx>
#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<SGGeod> GeodRow;
typedef std::vector<GeodRow> 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<rows; i++ ) {
m[i].resize( columns );
}
}
return m.GetNode(index);
}
inline void set( int col, int row, SGGeod p ) {
int index = ( row * _cols ) + col;
if ( col < 0 || col >= _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<double> 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<double> 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

55
src/Utils/poly2ogr/poly2ogr.cxx
View file

@ -36,6 +36,7 @@
#endif
#include <boost/foreach.hpp>
#include <ogrsf_frmts.h>
#include <simgear/debug/logstream.hxx>
#include <simgear/io/sg_binobj.hxx>
@ -45,9 +46,6 @@
#include <simgear/misc/sg_dir.hxx>
#include <Polygon/polygon.hxx>
#include <Polygon/point3d.hxx>
#include <ogrsf_frmts.h>
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<Point3D>& nodes) {
OGRLinearRing* make_ring_from_fan(const int_list& fan, const std::vector<SGGeod>& 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::vector<Point3D
}
for (;vertex!=fan.end();++vertex) {
OGRPoint *point=new OGRPoint();
const Point3D& node = nodes[*vertex];
point->setX(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<Point3D
return ring;
}
OGRLinearRing* make_ring_from_strip(const int_list& strip, const std::vector<Point3D>& nodes) {
OGRLinearRing* make_ring_from_strip(const int_list& strip, const std::vector<SGGeod>& nodes) {
OGRLinearRing* ring = new OGRLinearRing();
const size_t vertex_count = strip.size();
unsigned int i;
for (i=0;i<vertex_count;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);
}
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<Point3D>& wgs84_nodes) {
void convert_triangles(const std::string& path, const group_list& verts, const string_list& materials, const std::vector<SGGeod>& wgs84_nodes) {
const size_t groups_count = verts.size();
for (unsigned int i=0;i<groups_count;i++) {
@ -232,10 +230,10 @@ void convert_triangles(const std::string& path, const group_list& verts, const s
OGRLinearRing* ring = new OGRLinearRing();
for (int k=0;k<3;k++) {
OGRPoint *point=new OGRPoint();
const Point3D& node = wgs84_nodes[tri_verts[j+k]];
point->setX(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<Point3D>& wgs84_nodes) {
void convert_triangle_fans(const std::string& path, const group_list& verts, const string_list& materials, const std::vector<SGGeod>& wgs84_nodes) {
const size_t groups_count = verts.size();
for (unsigned int i=0;i<groups_count;i++) {
@ -255,7 +253,7 @@ void convert_triangle_fans(const std::string& path, const group_list& verts, con
}
}
void convert_triangle_strips(const std::string& path, const group_list& verts, const string_list& materials, const std::vector<Point3D>& wgs84_nodes) {
void convert_triangle_strips(const std::string& path, const group_list& verts, const string_list& materials, const std::vector<SGGeod>& wgs84_nodes) {
const size_t groups_count = verts.size();
for (unsigned int i=0;i<groups_count;i++) {
@ -273,17 +271,14 @@ void process_scenery_file(const std::string& path) {
SGVec3d gbs_center = binObject.get_gbs_center();
const std::vector<SGVec3d>& wgs84_nodes = binObject.get_wgs84_nodes();
std::vector<Point3D> geod_nodes;
std::vector<SGGeod> geod_nodes;
const size_t node_count = wgs84_nodes.size();
for (unsigned int i=0;i<node_count;i++) {
SGVec3d wgs84 = wgs84_nodes[i];
Point3D raw = Point3D( gbs_center.x() + wgs84.x(),
SGVec3d raw = SGVec3d( gbs_center.x() + wgs84.x(),
gbs_center.y() + wgs84.y(),
gbs_center.z() + wgs84.z() );
Point3D radians = sgCartToGeod(raw);
Point3D geod = Point3D( radians.x() / SGD_DEGREES_TO_RADIANS,
radians.y() / SGD_DEGREES_TO_RADIANS,
radians.z() );
SGGeod geod = SGGeod::fromCart( raw );
geod_nodes.push_back(geod);
}