]> git.mxchange.org Git - flightgear.git/blobdiff - DEM/dem.cxx
Dynamically update "error" until the resulting tile data scales within
[flightgear.git] / DEM / dem.cxx
index f59d4831e09309e0be8583b64138c1be0c94f360..e4661a3c70b7896816700732797e08c243a752da 100644 (file)
@@ -1,10 +1,8 @@
-// -*- Mode: C++ -*-
-//
-// dem.c -- DEM management class
+// dem.cxx -- DEM management class
 //
 // Written by Curtis Olson, started March 1998.
 //
-// Copyright (C) 1998  Curtis L. Olson  - curt@me.umn.edu
+// Copyright (C) 1998  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
 #include <math.h>     // rint()
 #include <stdio.h>
 #include <string.h>
+
 #ifdef HAVE_SYS_STAT_H
 #  include <sys/stat.h> // stat()
 #endif
+
 #ifdef FG_HAVE_STD_INCLUDES
 #  include <cerrno>
 #else
 #  include <errno.h>
 #endif
+
 #ifdef HAVE_UNISTD_H
 # include <unistd.h>   // stat()
 #endif
-#include <string>
 
-#include STL_IOSTREAM
-
-// #include <zlib/zlib.h>
 #include <Misc/fgstream.hxx>
 #include <Misc/strutils.hxx>
-#include <Include/compiler.h>
-FG_USING_NAMESPACE(std);
+#include <Include/fg_constants.h>
 
 #include "dem.hxx"
-#include "leastsqs.hxx"
-
-#include <Include/fg_constants.h>
 
 
 #define MAX_EX_NODES 10000
@@ -76,24 +69,25 @@ FG_USING_NAMESPACE(std);
 #endif //0
 
 
-fgDEM::fgDEM( void ) {
-    // cout << "class fgDEM CONstructor called." << endl;
+FGDem::FGDem( void ) {
+    // cout << "class FGDem CONstructor called." << endl;
     dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
     output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
 }
 
 
-fgDEM::fgDEM( const string &file ) {
-    // cout << "class fgDEM CONstructor called." << endl;
+FGDem::FGDem( const string &file ) {
+    // cout << "class FGDem CONstructor called." << endl;
     dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
     output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
 
-    fgDEM::open(file);
+    FGDem::open(file);
 }
 
 
 // open a DEM file
-int fgDEM::open ( const string& file ) {
+int
+FGDem::open ( const string& file ) {
     // open input file (or read from stdin)
     if ( file ==  "-" ) {
        printf("Loading DEM data file: stdin\n");
@@ -115,7 +109,8 @@ int fgDEM::open ( const string& file ) {
 
 
 // close a DEM file
-int fgDEM::close () {
+int
+FGDem::close () {
     // the fg_gzifstream doesn't seem to have a close()
 
     delete in;
@@ -125,7 +120,8 @@ int fgDEM::close () {
 
 
 // return next token from input stream
-string fgDEM::next_token() {
+string
+FGDem::next_token() {
     string token;
 
     *in >> token;
@@ -137,7 +133,8 @@ string fgDEM::next_token() {
 
 
 // return next integer from input stream
-int fgDEM::next_int() {
+int
+FGDem::next_int() {
     int result;
     
     *in >> result;
@@ -147,7 +144,8 @@ int fgDEM::next_int() {
 
 
 // return next double from input stream
-double fgDEM::next_double() {
+double
+FGDem::next_double() {
     double result;
 
     *in >> result;
@@ -157,15 +155,12 @@ double fgDEM::next_double() {
 
 
 // return next exponential num from input stream
-double fgDEM::next_exp() {
+double
+FGDem::next_exp() {
     string token;
-    double mantissa;
-    int exp, acc;
-    int i;
 
     token = next_token();
 
-#if 1
     const char* p = token.c_str();
     char buf[64];
     char* bp = buf;
@@ -179,34 +174,16 @@ double fgDEM::next_exp() {
     }
     *bp = 0;
     return ::atof( buf );
-#else
-    sscanf(token.c_str(), "%lfD%d", &mantissa, &exp);
-
-    // cout << "    Mantissa = " << mantissa << "  Exp = " << exp << "\n";
-
-    acc = 1;
-    if ( exp > 0 ) {
-       for ( i = 1; i <= exp; i++ ) {
-           acc *= 10;
-       }
-    } else if ( exp < 0 ) {
-       for ( i = -1; i >= exp; i-- ) {
-           acc /= 10;
-       }
-    }
-
-    return( (int)rint(mantissa * (double)acc) );
-#endif
 }
 
 
 // read and parse DEM "A" record
-int fgDEM::read_a_record() {
+int
+FGDem::read_a_record() {
     int i, inum;
     double dnum;
     string name, token;
     char c;
-    char *ptr;
 
     // get the name field (144 characters)
     for ( i = 0; i < 144; i++ ) {
@@ -306,32 +283,10 @@ int fgDEM::read_a_record() {
     i = token.length();
     cout << "    length = " << i << "\n";
 
-#if 1
     inum = atoi( token.substr( 0, i - 36 ) );
-    row_step = atof( token.substr( i - 36, 12 ) );
-    col_step = atof( token.substr( i - 24, 12 ) );
-    //token.substr( 25, 12 )
-#else
-    ptr = token.c_str() + i - 12;
-    cout << "    last field = " << ptr << " = " << atof(ptr) << "\n";
-    ptr[0] = '\0';
-
-    ptr = ptr - 12;
-    col_step = atof(ptr);
-    cout << "    last field = " << ptr << " = " << col_step << "\n";
-    ptr[0] = '\0';
-
-    ptr = ptr - 12;
-    row_step = atof(ptr);
-    cout << "    last field = " << ptr << " = " << row_step << "\n";
-    ptr[0] = '\0';
-
-    // accuracy code = atod(token)
-    ptr = ptr - 12;
-    inum = atoi(ptr);
-#endif
+    row_step = atof( token.substr( i - 24, 12 ) );
+    col_step = atof( token.substr( i - 36, 12 ) );
     cout << "    Accuracy code = " << inum << "\n";
-
     cout << "    column step = " << col_step << 
        "  row step = " << row_step << "\n";
 
@@ -347,7 +302,8 @@ int fgDEM::read_a_record() {
 
 
 // read and parse DEM "B" record
-void fgDEM::read_b_record( ) {
+void
+FGDem::read_b_record( ) {
     string token;
     int i;
 
@@ -384,7 +340,8 @@ void fgDEM::read_b_record( ) {
 
 
 // parse dem file
-int fgDEM::parse( ) {
+int
+FGDem::parse( ) {
     int i;
 
     cur_col = 0;
@@ -409,9 +366,85 @@ int fgDEM::parse( ) {
 }
 
 
-// return the current altitude based on mesh data.  We should rewrite
+// write out the area of data covered by the specified bucket.  Data
+// is written out column by column starting at the lower left hand
+// corner.
+int
+FGDem::write_area( const string& root, FGBucket& b, bool compress ) {
+    // calculate some boundaries
+    double min_x = ( b.get_center_lon() - 0.5 * b.get_width() ) * 3600.0;
+    double max_x = ( b.get_center_lon() + 0.5 * b.get_width() ) * 3600.0;
+
+    double min_y = ( b.get_center_lat() - 0.5 * b.get_height() ) * 3600.0;
+    double max_y = ( b.get_center_lat() + 0.5 * b.get_height() ) * 3600.0;
+
+    cout << b << endl;
+    cout << "width = " << b.get_width() << " height = " << b.get_height() 
+        << endl;
+
+    int start_x = (int)((min_x - originx) / col_step);
+    int span_x = (int)(b.get_width() * 3600.0 / col_step);
+
+    int start_y = (int)((min_y - originy) / row_step);
+    int span_y = (int)(b.get_height() * 3600.0 / row_step);
+
+    cout << "start_x = " << start_x << "  span_x = " << span_x << endl;
+    cout << "start_y = " << start_y << "  span_y = " << span_y << endl;
+
+    // Do a simple sanity checking.  But, please, please be nice to
+    // this write_area() routine and feed it buckets that coincide
+    // well with the underlying grid structure and spacing.
+
+    if ( ( min_x < originx )
+        || ( max_x > originx + cols * col_step )
+        || ( min_y < originy )
+        || ( max_y > originy + rows * row_step ) ) {
+       cout << "  ERROR: bucket at least partially outside DEM data range!" <<
+           endl;
+       return 0;
+    }
+
+    // generate output file name
+    string base = b.gen_base_path();
+    string path = root + "/Scenery/" + base;
+    string command = "mkdir -p " + path;
+    system( command.c_str() );
+
+    string demfile = path + "/" + b.gen_index_str() + ".dem";
+    cout << "demfile = " << demfile << endl;
+
+    // write the file
+    FILE *fp;
+    if ( (fp = fopen(demfile.c_str(), "w")) == NULL ) {
+       cout << "cannot open " << demfile << " for writing!" << endl;
+       exit(-1);
+    }
+
+    fprintf( fp, "%d %d\n", (int)min_x, (int)min_y );
+    fprintf( fp, "%d %d %d %d\n", span_x + 1, (int)col_step, 
+            span_y + 1, (int)row_step );
+    for ( int i = start_x; i <= start_x + span_x; ++i ) {
+       for ( int j = start_y; j <= start_y + span_y; ++j ) {
+           fprintf( fp, "%d ", (int)dem_data[i][j] );
+       }
+       fprintf( fp, "\n" );
+    }
+    fclose(fp);
+
+    if ( compress ) {
+       string command = "gzip --best -f " + demfile;
+       system( command.c_str() );
+    }
+
+    return 1;
+}
+
+
+#if 0
+
+// return the current altitude based on grid data.  We should rewrite
 // this to interpolate exact values, but for now this is good enough
-double fgDEM::interpolate_altitude( double lon, double lat ) {
+double FGDem::interpolate_altitude( double lon, double lat ) {
     // we expect incoming (lon,lat) to be in arcsec for now
 
     double xlocal, ylocal, dx, dy, zA, zB, elev;
@@ -524,7 +557,7 @@ double fgDEM::interpolate_altitude( double lon, double lat ) {
 
 
 // Use least squares to fit a simpler data set to dem data
-void fgDEM::fit( double error, fgBUCKET *p ) {
+void FGDem::fit( double error, FGBucket& p ) {
     double x[DEM_SIZE_1], y[DEM_SIZE_1];
     double m, b, ave_error, max_error;
     double cury, lasty;
@@ -537,9 +570,9 @@ void fgDEM::fit( double error, fgBUCKET *p ) {
     outputmesh_init();
 
     // determine dimensions
-    colmin = p->x * ( (cols - 1) / 8);
+    colmin = p.get_x() * ( (cols - 1) / 8);
     colmax = colmin + ( (cols - 1) / 8);
-    rowmin = p->y * ( (rows - 1) / 8);
+    rowmin = p.get_y() * ( (rows - 1) / 8);
     rowmax = rowmin + ( (rows - 1) / 8);
     printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
     
@@ -660,7 +693,7 @@ void fgDEM::fit( double error, fgBUCKET *p ) {
 
 
 // Initialize output mesh structure
-void fgDEM::outputmesh_init( void ) {
+void FGDem::outputmesh_init( void ) {
     int i, j;
     
     for ( j = 0; j < DEM_SIZE_1; j++ ) {
@@ -672,13 +705,13 @@ void fgDEM::outputmesh_init( void ) {
 
 
 // Get the value of a mesh node
-double fgDEM::outputmesh_get_pt( int i, int j ) {
+double FGDem::outputmesh_get_pt( int i, int j ) {
     return ( output_data[i][j] );
 }
 
 
 // Set the value of a mesh node
-void fgDEM::outputmesh_set_pt( int i, int j, double value ) {
+void FGDem::outputmesh_set_pt( int i, int j, double value ) {
     // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
    output_data[i][j] = value;
 }
@@ -688,11 +721,12 @@ void fgDEM::outputmesh_set_pt( int i, int j, double value ) {
 // Check for an optional "index.node.ex" file in case there is a .poly
 // file to go along with this node file.  Include these nodes first
 // since they are referenced by position from the .poly file.
-void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
+void FGDem::outputmesh_output_nodes( const string& fg_root, FGBucket& p )
+{
     double exnodes[MAX_EX_NODES][3];
     struct stat stat_buf;
     string dir;
-    char base_path[256], file[256], exfile[256];
+    char file[256], exfile[256];
 #ifdef WIN32
     char tmp_path[256];
 #endif
@@ -703,15 +737,15 @@ void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
     int i, j, count, excount, result;
 
     // determine dimensions
-    colmin = p->x * ( (cols - 1) / 8);
+    colmin = p.get_x() * ( (cols - 1) / 8);
     colmax = colmin + ( (cols - 1) / 8);
-    rowmin = p->y * ( (rows - 1) / 8);
+    rowmin = p.get_y() * ( (rows - 1) / 8);
     rowmax = rowmin + ( (rows - 1) / 8);
     cout << "  dumping region = " << colmin << "," << rowmin << " to " <<
        colmax << "," << rowmax << "\n";
 
     // generate the base directory
-    fgBucketGenBasePath(p, base_path);
+    string base_path = p.gen_base_path();
     cout << "fg_root = " << fg_root << "  Base Path = " << base_path << endl;
     dir = fg_root + "/Scenery/" + base_path;
     cout << "Dir = " << dir << endl;
@@ -752,7 +786,7 @@ void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
     }
 
     // get index and generate output file name
-    index = fgBucketGenIndex(p);
+    index = p.gen_index();
     sprintf(file, "%s/%ld.node", dir.c_str(), index);
 
     // get (optional) extra node file name (in case there is matching
@@ -820,16 +854,32 @@ void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
 
     fclose(fd);
 }
+#endif
 
 
-fgDEM::~fgDEM( void ) {
-    // printf("class fgDEM DEstructor called.\n");
+FGDem::~FGDem( void ) {
+    // printf("class FGDem DEstructor called.\n");
     delete [] dem_data;
     delete [] output_data;
 }
 
 
 // $Log$
+// Revision 1.27  1999/03/25 19:04:36  curt
+// Minor tweaks related to FGBucket usage.
+//
+// Revision 1.26  1999/03/13 17:40:37  curt
+// Moved point interpolation and least squares fitting to contruction program
+// area.
+// Moved leastsqs.* to Lib/Math/
+//
+// Revision 1.25  1999/03/12 22:53:07  curt
+// Added a routine to dump out the portion of the dem data covered by a
+// specified bucket.  Other changes related to needs of scenery tools overhaul.
+//
+// Revision 1.24  1999/03/11 23:31:56  curt
+// Tweaks to use newbucket.hxx
+//
 // Revision 1.23  1999/03/10 01:09:12  curt
 // Tweaks to go along with scenery tools overhaul.
 // Added a new constructor that accepts the file name.