DEM/dem.cxx

   1 // -*- Mode: C++ -*-
   2 //
   3 // dem.c -- DEM management class
   4 //
   5 // Written by Curtis Olson, started March 1998.
   6 //
   7 // Copyright (C) 1998  Curtis L. Olson  - curt@me.umn.edu
   8 //
   9 // This program is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU General Public License as
  11 // published by the Free Software Foundation; either version 2 of the
  12 // License, or (at your option) any later version.
  13 //
  14 // This program is distributed in the hope that it will be useful, but
  15 // WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17 // General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU General Public License
  20 // along with this program; if not, write to the Free Software
  21 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  22 //
  23 // $Id$
  24 // (Log is kept at end of this file)
  25
  26
  27 #ifdef HAVE_CONFIG_H
  28 #  include <config.h>
  29 #endif
  30
  31 #include <Include/compiler.h>
  32
  33 #include <ctype.h>    // isspace()
  34 #include <stdlib.h>   // atoi()
  35 #include <math.h>     // rint()
  36 #include <stdio.h>
  37 #include <string.h>
  38 #ifdef HAVE_SYS_STAT_H
  39 #  include <sys/stat.h> // stat()
  40 #endif
  41 #ifdef FG_HAVE_STD_INCLUDES
  42 #  include <cerrno>
  43 #else
  44 #  include <errno.h>
  45 #endif
  46 #ifdef HAVE_UNISTD_H
  47 # include <unistd.h>   // stat()
  48 #endif
  49 #include <string>
  50
  51 #include STL_IOSTREAM
  52
  53 // #include <zlib/zlib.h>
  54 #include <Misc/fgstream.hxx>
  55 #include <Misc/strutils.hxx>
  56 #include <Include/compiler.h>
  57
  58 FG_USING_NAMESPACE(std);
  59
  60 #include "dem.hxx"
  61 // #include "leastsqs.hxx"
  62
  63 #include <Include/fg_constants.h>
  64
  65
  66 #define MAX_EX_NODES 10000
  67
  68 #if 0
  69 #ifdef WIN32
  70 # ifdef __BORLANDC__
  71 #  include <dir.h>
  72 #  define MKDIR(a) mkdir(a)
  73 # else
  74 #  define MKDIR(a) mkdir(a,S_IRWXU)  // I am just guessing at this flag (NHV)
  75 # endif // __BORLANDC__
  76 #endif // WIN32
  77 #endif //0
  78
  79
  80 FGDem::FGDem( void ) {
  81     // cout << "class FGDem CONstructor called." << endl;
  82     dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
  83     output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
  84 }
  85
  86
  87 FGDem::FGDem( const string &file ) {
  88     // cout << "class FGDem CONstructor called." << endl;
  89     dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
  90     output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
  91
  92     FGDem::open(file);
  93 }
  94
  95
  96 // open a DEM file
  97 int
  98 FGDem::open ( const string& file ) {
  99     // open input file (or read from stdin)
 100     if ( file ==  "-" ) {
 101         printf("Loading DEM data file: stdin\n");
 102         // fd = stdin;
 103         // fd = gzdopen(STDIN_FILENO, "r");
 104         printf("Not yet ported ...\n");
 105         return 0;
 106     } else {
 107         in = new fg_gzifstream( file );
 108         if ( !(*in) ) {
 109             cout << "Cannot open " << file << endl;
 110             return 0;
 111         }
 112         cout << "Loading DEM data file: " << file << endl;
 113     }
 114
 115     return 1;
 116 }
 117
 118
 119 // close a DEM file
 120 int
 121 FGDem::close () {
 122     // the fg_gzifstream doesn't seem to have a close()
 123
 124     delete in;
 125
 126     return 1;
 127 }
 128
 129
 130 // return next token from input stream
 131 string
 132 FGDem::next_token() {
 133     string token;
 134
 135     *in >> token;
 136
 137     // cout << "    returning " + token + "\n";
 138
 139     return token;
 140 }
 141
 142
 143 // return next integer from input stream
 144 int
 145 FGDem::next_int() {
 146     int result;
 147
 148     *in >> result;
 149
 150     return result;
 151 }
 152
 153
 154 // return next double from input stream
 155 double
 156 FGDem::next_double() {
 157     double result;
 158
 159     *in >> result;
 160
 161     return result;
 162 }
 163
 164
 165 // return next exponential num from input stream
 166 double
 167 FGDem::next_exp() {
 168     string token;
 169     double mantissa;
 170     int exp, acc;
 171     int i;
 172
 173     token = next_token();
 174
 175 #if 1
 176     const char* p = token.c_str();
 177     char buf[64];
 178     char* bp = buf;
 179
 180     for ( ; *p != 0; ++p )
 181     {
 182         if ( *p == 'D' )
 183             *bp++ = 'E';
 184         else
 185             *bp++ = *p;
 186     }
 187     *bp = 0;
 188     return ::atof( buf );
 189 #else
 190     sscanf(token.c_str(), "%lfD%d", &mantissa, &exp);
 191
 192     // cout << "    Mantissa = " << mantissa << "  Exp = " << exp << "\n";
 193
 194     acc = 1;
 195     if ( exp > 0 ) {
 196         for ( i = 1; i <= exp; i++ ) {
 197             acc *= 10;
 198         }
 199     } else if ( exp < 0 ) {
 200         for ( i = -1; i >= exp; i-- ) {
 201             acc /= 10;
 202         }
 203     }
 204
 205     return( (int)rint(mantissa * (double)acc) );
 206 #endif
 207 }
 208
 209
 210 // read and parse DEM "A" record
 211 int
 212 FGDem::read_a_record() {
 213     int i, inum;
 214     double dnum;
 215     string name, token;
 216     char c;
 217     char *ptr;
 218
 219     // get the name field (144 characters)
 220     for ( i = 0; i < 144; i++ ) {
 221         in->get(c);
 222         name += c;
 223     }
 224
 225     // clean off the trailing whitespace
 226     name = trim(name);
 227     cout << "    Quad name field: " << name << endl;
 228
 229     // DEM level code, 3 reflects processing by DMA
 230     inum = next_int();
 231     cout << "    DEM level code = " << inum << "\n";
 232
 233     if ( inum > 3 ) {
 234         return 0;
 235     }
 236
 237     // Pattern code, 1 indicates a regular elevation pattern
 238     inum = next_int();
 239     cout << "    Pattern code = " << inum << "\n";
 240
 241     // Planimetric reference system code, 0 indicates geographic
 242     // coordinate system.
 243     inum = next_int();
 244     cout << "    Planimetric reference code = " << inum << "\n";
 245
 246     // Zone code
 247     inum = next_int();
 248     cout << "    Zone code = " << inum << "\n";
 249
 250     // Map projection parameters (ignored)
 251     for ( i = 0; i < 15; i++ ) {
 252         dnum = next_exp();
 253         // printf("%d: %f\n",i,dnum);
 254     }
 255
 256     // Units code, 3 represents arc-seconds as the unit of measure for
 257     // ground planimetric coordinates throughout the file.
 258     inum = next_int();
 259     if ( inum != 3 ) {
 260         cout << "    Unknown (X,Y) units code = " << inum << "!\n";
 261         exit(-1);
 262     }
 263
 264     // Units code; 2 represents meters as the unit of measure for
 265     // elevation coordinates throughout the file.
 266     inum = next_int();
 267     if ( inum != 2 ) {
 268         cout << "    Unknown (Z) units code = " << inum << "!\n";
 269         exit(-1);
 270     }
 271
 272     // Number (n) of sides in the polygon which defines the coverage of
 273     // the DEM file (usually equal to 4).
 274     inum = next_int();
 275     if ( inum != 4 ) {
 276         cout << "    Unknown polygon dimension = " << inum << "!\n";
 277         exit(-1);
 278     }
 279
 280     // Ground coordinates of bounding box in arc-seconds
 281     dem_x1 = originx = next_exp();
 282     dem_y1 = originy = next_exp();
 283     cout << "    Origin = (" << originx << "," << originy << ")\n";
 284
 285     dem_x2 = next_exp();
 286     dem_y2 = next_exp();
 287
 288     dem_x3 = next_exp();
 289     dem_y3 = next_exp();
 290
 291     dem_x4 = next_exp();
 292     dem_y4 = next_exp();
 293
 294     // Minimum/maximum elevations in meters
 295     dem_z1 = next_exp();
 296     dem_z2 = next_exp();
 297     cout << "    Elevation range " << dem_z1 << " to " << dem_z2 << "\n";
 298
 299     // Counterclockwise angle from the primary axis of ground
 300     // planimetric referenced to the primary axis of the DEM local
 301     // reference system.
 302     token = next_token();
 303
 304     // Accuracy code; 0 indicates that a record of accuracy does not
 305     // exist and that no record type C will follow.
 306
 307     // DEM spacial resolution.  Usually (3,3,1) (3,6,1) or (3,9,1)
 308     // depending on latitude
 309
 310     // I will eventually have to do something with this for data at
 311     // higher latitudes */
 312     token = next_token();
 313     cout << "    accuracy & spacial resolution string = " << token << endl;
 314     i = token.length();
 315     cout << "    length = " << i << "\n";
 316
 317     inum = atoi( token.substr( 0, i - 36 ) );
 318     row_step = atof( token.substr( i - 24, 12 ) );
 319     col_step = atof( token.substr( i - 36, 12 ) );
 320     cout << "    Accuracy code = " << inum << "\n";
 321     cout << "    column step = " << col_step <<
 322         "  row step = " << row_step << "\n";
 323
 324     // dimension of arrays to follow (1)
 325     token = next_token();
 326
 327     // number of profiles
 328     dem_num_profiles = cols = next_int();
 329     cout << "    Expecting " << dem_num_profiles << " profiles\n";
 330
 331     return 1;
 332 }
 333
 334
 335 // read and parse DEM "B" record
 336 void
 337 FGDem::read_b_record( ) {
 338     string token;
 339     int i;
 340
 341     // row / column id of this profile
 342     prof_row = next_int();
 343     prof_col = next_int();
 344     // printf("col id = %d  row id = %d\n", prof_col, prof_row);
 345
 346     // Number of columns and rows (elevations) in this profile
 347     prof_num_rows = rows = next_int();
 348     prof_num_cols = next_int();
 349     // printf("    profile num rows = %d\n", prof_num_rows);
 350
 351     // Ground planimetric coordinates (arc-seconds) of the first
 352     // elevation in the profile
 353     prof_x1 = next_exp();
 354     prof_y1 = next_exp();
 355     // printf("    Starting at %.2f %.2f\n", prof_x1, prof_y1);
 356
 357     // Elevation of local datum for the profile.  Always zero for
 358     // 1-degree DEM, the reference is mean sea level.
 359     token = next_token();
 360
 361     // Minimum and maximum elevations for the profile.
 362     token = next_token();
 363     token = next_token();
 364
 365     // One (usually) dimensional array (prof_num_cols,1) of elevations
 366     for ( i = 0; i < prof_num_rows; i++ ) {
 367         prof_data = next_int();
 368         dem_data[cur_col][i] = (float)prof_data;
 369     }
 370 }
 371
 372
 373 // parse dem file
 374 int
 375 FGDem::parse( ) {
 376     int i;
 377
 378     cur_col = 0;
 379
 380     if ( !read_a_record() ) {
 381         return(0);
 382     }
 383
 384     for ( i = 0; i < dem_num_profiles; i++ ) {
 385         // printf("Ready to read next b record\n");
 386         read_b_record();
 387         cur_col++;
 388
 389         if ( cur_col % 100 == 0 ) {
 390             cout << "    loaded " << cur_col << " profiles of data\n";
 391         }
 392     }
 393
 394     cout << "    Done parsing\n";
 395
 396     return 1;
 397 }
 398
 399
 400 // write out the area of data covered by the specified bucket.  Data
 401 // is written out column by column starting at the lower left hand
 402 // corner.
 403 int
 404 FGDem::write_area( const string& root, FGBucket& b, bool compress ) {
 405     char tile_name[256];
 406
 407     // calculate some boundaries
 408     double min_x = ( b.get_center_lon() - 0.5 * b.get_width() ) * 3600.0;
 409     double max_x = ( b.get_center_lon() + 0.5 * b.get_width() ) * 3600.0;
 410
 411     double min_y = ( b.get_center_lat() - 0.5 * b.get_height() ) * 3600.0;
 412     double max_y = ( b.get_center_lat() + 0.5 * b.get_height() ) * 3600.0;
 413
 414     cout << b << endl;
 415     cout << "width = " << b.get_width() << " height = " << b.get_height()
 416          << endl;
 417
 418     int start_x = (int)((min_x - originx) / col_step);
 419     int span_x = (int)(b.get_width() * 3600.0 / col_step);
 420
 421     int start_y = (int)((min_y - originy) / row_step);
 422     int span_y = (int)(b.get_height() * 3600.0 / row_step);
 423
 424     cout << "start_x = " << start_x << "  span_x = " << span_x << endl;
 425     cout << "start_y = " << start_y << "  span_y = " << span_y << endl;
 426
 427     // Do a simple sanity checking.  But, please, please be nice to
 428     // this write_area() routine and feed it buckets that coincide
 429     // well with the underlying grid structure and spacing.
 430
 431     if ( ( min_x < originx )
 432          || ( max_x > originx + cols * col_step )
 433          || ( min_y < originy )
 434          || ( max_y > originy + rows * row_step ) ) {
 435         cout << "  ERROR: bucket at least partially outside DEM data range!" <<
 436             endl;
 437         return -1;
 438     }
 439
 440     // generate output file name
 441     string base = b.gen_base_path();
 442     string path = root + "/Scenery/" + base;
 443     string command = "mkdir -p " + path;
 444     system( command.c_str() );
 445
 446     long int b_index = b.gen_index();
 447     sprintf(tile_name, "%ld", b_index);
 448     string demfile = path + "/" + tile_name + ".dem";
 449     cout << "demfile = " << demfile << endl;
 450
 451     // write the file
 452     FILE *fp;
 453     if ( (fp = fopen(demfile.c_str(), "w")) == NULL ) {
 454         cout << "cannot open " << demfile << " for writing!" << endl;
 455         exit(-1);
 456     }
 457
 458     fprintf( fp, "%d %d\n", (int)min_x, (int)min_y );
 459     fprintf( fp, "%d %d %d %d\n", span_x + 1, (int)col_step,
 460              span_y + 1, (int)row_step );
 461     for ( int i = start_x; i <= start_x + span_x; ++i ) {
 462         for ( int j = start_y; j <= start_y + span_y; ++j ) {
 463             fprintf( fp, "%d ", (int)dem_data[i][j] );
 464         }
 465         fprintf( fp, "\n" );
 466     }
 467     fclose(fp);
 468
 469     if ( compress ) {
 470         string command = "gzip --best -f " + demfile;
 471         system( command.c_str() );
 472     }
 473 }
 474
 475
 476 #if 0
 477
 478 // return the current altitude based on grid data.  We should rewrite
 479 // this to interpolate exact values, but for now this is good enough
 480 double FGDem::interpolate_altitude( double lon, double lat ) {
 481     // we expect incoming (lon,lat) to be in arcsec for now
 482
 483     double xlocal, ylocal, dx, dy, zA, zB, elev;
 484     int x1, x2, x3, y1, y2, y3;
 485     float z1, z2, z3;
 486     int xindex, yindex;
 487
 488     /* determine if we are in the lower triangle or the upper triangle
 489        ______
 490        |   /|
 491        |  / |
 492        | /  |
 493        |/   |
 494        ------
 495
 496        then calculate our end points
 497      */
 498
 499     xlocal = (lon - originx) / col_step;
 500     ylocal = (lat - originy) / row_step;
 501
 502     xindex = (int)(xlocal);
 503     yindex = (int)(ylocal);
 504
 505     // printf("xindex = %d  yindex = %d\n", xindex, yindex);
 506
 507     if ( xindex + 1 == cols ) {
 508         xindex--;
 509     }
 510
 511     if ( yindex + 1 == rows ) {
 512         yindex--;
 513     }
 514
 515     if ( (xindex < 0) || (xindex + 1 >= cols) ||
 516          (yindex < 0) || (yindex + 1 >= rows) ) {
 517         return(-9999);
 518     }
 519
 520     dx = xlocal - xindex;
 521     dy = ylocal - yindex;
 522
 523     if ( dx > dy ) {
 524         // lower triangle
 525         // printf("  Lower triangle\n");
 526
 527         x1 = xindex;
 528         y1 = yindex;
 529         z1 = dem_data[x1][y1];
 530
 531         x2 = xindex + 1;
 532         y2 = yindex;
 533         z2 = dem_data[x2][y2];
 534
 535         x3 = xindex + 1;
 536         y3 = yindex + 1;
 537         z3 = dem_data[x3][y3];
 538
 539         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
 540         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
 541         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
 542         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
 543
 544         zA = dx * (z2 - z1) + z1;
 545         zB = dx * (z3 - z1) + z1;
 546
 547         // printf("  zA = %.2f  zB = %.2f\n", zA, zB);
 548
 549         if ( dx > FG_EPSILON ) {
 550             elev = dy * (zB - zA) / dx + zA;
 551         } else {
 552             elev = zA;
 553         }
 554     } else {
 555         // upper triangle
 556         // printf("  Upper triangle\n");
 557
 558         x1 = xindex;
 559         y1 = yindex;
 560         z1 = dem_data[x1][y1];
 561
 562         x2 = xindex;
 563         y2 = yindex + 1;
 564         z2 = dem_data[x2][y2];
 565
 566         x3 = xindex + 1;
 567         y3 = yindex + 1;
 568         z3 = dem_data[x3][y3];
 569
 570         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
 571         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
 572         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
 573         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
 574
 575         zA = dy * (z2 - z1) + z1;
 576         zB = dy * (z3 - z1) + z1;
 577
 578         // printf("  zA = %.2f  zB = %.2f\n", zA, zB );
 579         // printf("  xB - xA = %.2f\n", col_step * dy / row_step);
 580
 581         if ( dy > FG_EPSILON ) {
 582             elev = dx * (zB - zA) / dy    + zA;
 583         } else {
 584             elev = zA;
 585         }
 586     }
 587
 588     return(elev);
 589 }
 590
 591
 592 // Use least squares to fit a simpler data set to dem data
 593 void FGDem::fit( double error, FGBucket& p ) {
 594     double x[DEM_SIZE_1], y[DEM_SIZE_1];
 595     double m, b, ave_error, max_error;
 596     double cury, lasty;
 597     int n, row, start, end;
 598     int colmin, colmax, rowmin, rowmax;
 599     bool good_fit;
 600     // FILE *dem, *fit, *fit1;
 601
 602     printf("Initializing output mesh structure\n");
 603     outputmesh_init();
 604
 605     // determine dimensions
 606     colmin = p.get_x() * ( (cols - 1) / 8);
 607     colmax = colmin + ( (cols - 1) / 8);
 608     rowmin = p.get_y() * ( (rows - 1) / 8);
 609     rowmax = rowmin + ( (rows - 1) / 8);
 610     printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
 611
 612     // include the corners explicitly
 613     outputmesh_set_pt(colmin, rowmin, dem_data[colmin][rowmin]);
 614     outputmesh_set_pt(colmin, rowmax, dem_data[colmin][rowmax]);
 615     outputmesh_set_pt(colmax, rowmax, dem_data[colmax][rowmax]);
 616     outputmesh_set_pt(colmax, rowmin, dem_data[colmax][rowmin]);
 617
 618     printf("Beginning best fit procedure\n");
 619
 620     for ( row = rowmin; row <= rowmax; row++ ) {
 621         // fit  = fopen("fit.dat",  "w");
 622         // fit1 = fopen("fit1.dat", "w");
 623
 624         start = colmin;
 625
 626         // printf("    fitting row = %d\n", row);
 627
 628         while ( start < colmax ) {
 629             end = start + 1;
 630             good_fit = true;
 631
 632             x[(end - start) - 1] = 0.0 + ( start * col_step );
 633             y[(end - start) - 1] = dem_data[start][row];
 634
 635             while ( (end <= colmax) && good_fit ) {
 636                 n = (end - start) + 1;
 637                 // printf("Least square of first %d points\n", n);
 638                 x[end - start] = 0.0 + ( end * col_step );
 639                 y[end - start] = dem_data[end][row];
 640                 least_squares(x, y, n, &m, &b);
 641                 ave_error = least_squares_error(x, y, n, m, b);
 642                 max_error = least_squares_max_error(x, y, n, m, b);
 643
 644                 /*
 645                 printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n",
 646                 start, end, ave_error, max_error, m, b);
 647
 648                 f = fopen("gnuplot.dat", "w");
 649                 for ( j = 0; j <= end; j++) {
 650                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 651                             dem_data[row][j]);
 652                 }
 653                 for ( j = start; j <= end; j++) {
 654                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 655                             dem_data[row][j]);
 656                 }
 657                 fclose(f);
 658
 659                 printf("Please hit return: "); gets(junk);
 660                 */
 661
 662                 if ( max_error > error ) {
 663                     good_fit = false;
 664                 }
 665
 666                 end++;
 667             }
 668
 669             if ( !good_fit ) {
 670                 // error exceeded the threshold, back up
 671                 end -= 2;  // back "end" up to the last good enough fit
 672                 n--;       // back "n" up appropriately too
 673             } else {
 674                 // we popped out of the above loop while still within
 675                 // the error threshold, so we must be at the end of
 676                 // the data set
 677                 end--;
 678             }
 679
 680             least_squares(x, y, n, &m, &b);
 681             ave_error = least_squares_error(x, y, n, m, b);
 682             max_error = least_squares_max_error(x, y, n, m, b);
 683
 684             /*
 685             printf("\n");
 686             printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n",
 687                    start, end, ave_error, max_error, m, b);
 688             printf("\n");
 689
 690             fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
 691             fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
 692             */
 693
 694             if ( start > colmin ) {
 695                 // skip this for the first line segment
 696                 cury = m * x[0] + b;
 697                 outputmesh_set_pt(start, row, (lasty + cury) / 2);
 698                 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
 699             }
 700
 701             lasty = m * x[end-start] + b;
 702             start = end;
 703         }
 704
 705         /*
 706         fclose(fit);
 707         fclose(fit1);
 708
 709         dem = fopen("gnuplot.dat", "w");
 710         for ( j = 0; j < DEM_SIZE_1; j++) {
 711             fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 712                     dem_data[j][row]);
 713         }
 714         fclose(dem);
 715         */
 716
 717         // NOTICE, this is for testing only.  This instance of
 718         // output_nodes should be removed.  It should be called only
 719         // once at the end once all the nodes have been generated.
 720         // newmesh_output_nodes(&nm, "mesh.node");
 721         // printf("Please hit return: "); gets(junk);
 722     }
 723
 724     // outputmesh_output_nodes(fg_root, p);
 725 }
 726
 727
 728 // Initialize output mesh structure
 729 void FGDem::outputmesh_init( void ) {
 730     int i, j;
 731
 732     for ( j = 0; j < DEM_SIZE_1; j++ ) {
 733         for ( i = 0; i < DEM_SIZE_1; i++ ) {
 734             output_data[i][j] = -9999.0;
 735         }
 736     }
 737 }
 738
 739
 740 // Get the value of a mesh node
 741 double FGDem::outputmesh_get_pt( int i, int j ) {
 742     return ( output_data[i][j] );
 743 }
 744
 745
 746 // Set the value of a mesh node
 747 void FGDem::outputmesh_set_pt( int i, int j, double value ) {
 748     // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
 749    output_data[i][j] = value;
 750 }
 751
 752
 753 // Write out a node file that can be used by the "triangle" program.
 754 // Check for an optional "index.node.ex" file in case there is a .poly
 755 // file to go along with this node file.  Include these nodes first
 756 // since they are referenced by position from the .poly file.
 757 void FGDem::outputmesh_output_nodes( const string& fg_root, FGBucket& p )
 758 {
 759     double exnodes[MAX_EX_NODES][3];
 760     struct stat stat_buf;
 761     string dir;
 762     char file[256], exfile[256];
 763 #ifdef WIN32
 764     char tmp_path[256];
 765 #endif
 766     string command;
 767     FILE *fd;
 768     long int index;
 769     int colmin, colmax, rowmin, rowmax;
 770     int i, j, count, excount, result;
 771
 772     // determine dimensions
 773     colmin = p.get_x() * ( (cols - 1) / 8);
 774     colmax = colmin + ( (cols - 1) / 8);
 775     rowmin = p.get_y() * ( (rows - 1) / 8);
 776     rowmax = rowmin + ( (rows - 1) / 8);
 777     cout << "  dumping region = " << colmin << "," << rowmin << " to " <<
 778         colmax << "," << rowmax << "\n";
 779
 780     // generate the base directory
 781     string base_path = p.gen_base_path();
 782     cout << "fg_root = " << fg_root << "  Base Path = " << base_path << endl;
 783     dir = fg_root + "/Scenery/" + base_path;
 784     cout << "Dir = " << dir << endl;
 785
 786     // stat() directory and create if needed
 787     errno = 0;
 788     result = stat(dir.c_str(), &stat_buf);
 789     if ( result != 0 && errno == ENOENT ) {
 790         cout << "Creating directory\n";
 791
 792 // #ifndef WIN32
 793
 794         command = "mkdir -p " + dir + "\n";
 795         system( command.c_str() );
 796
 797 #if 0
 798 // #else // WIN32
 799
 800         // Cygwin crashes when trying to output to node file
 801         // explicitly making directory structure seems OK on Win95
 802
 803         extract_path (base_path, tmp_path);
 804
 805         dir = fg_root + "/Scenery";
 806         if (my_mkdir ( dir.c_str() )) { exit (-1); }
 807
 808         dir = fg_root + "/Scenery/" + tmp_path;
 809         if (my_mkdir ( dir.c_str() )) { exit (-1); }
 810
 811         dir = fg_root + "/Scenery/" + base_path;
 812         if (my_mkdir ( dir.c_str() )) { exit (-1); }
 813
 814 // #endif // WIN32
 815 #endif //0
 816
 817     } else {
 818         // assume directory exists
 819     }
 820
 821     // get index and generate output file name
 822     index = p.gen_index();
 823     sprintf(file, "%s/%ld.node", dir.c_str(), index);
 824
 825     // get (optional) extra node file name (in case there is matching
 826     // .poly file.
 827     strcpy(exfile, file);
 828     strcat(exfile, ".ex");
 829
 830     // load extra nodes if they exist
 831     excount = 0;
 832     if ( (fd = fopen(exfile, "r")) != NULL ) {
 833         int junki;
 834         fscanf(fd, "%d %d %d %d", &excount, &junki, &junki, &junki);
 835
 836         if ( excount > MAX_EX_NODES - 1 ) {
 837             printf("Error, too many 'extra' nodes, increase array size\n");
 838             exit(-1);
 839         } else {
 840             printf("    Expecting %d 'extra' nodes\n", excount);
 841         }
 842
 843         for ( i = 1; i <= excount; i++ ) {
 844             fscanf(fd, "%d %lf %lf %lf\n", &junki,
 845                    &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]);
 846             printf("(extra) %d %.2f %.2f %.2f\n",
 847                     i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
 848         }
 849         fclose(fd);
 850     }
 851
 852     printf("Creating node file:  %s\n", file);
 853     fd = fopen(file, "w");
 854
 855     // first count regular nodes to generate header
 856     count = 0;
 857     for ( j = rowmin; j <= rowmax; j++ ) {
 858         for ( i = colmin; i <= colmax; i++ ) {
 859             if ( output_data[i][j] > -9000.0 ) {
 860                 count++;
 861             }
 862         }
 863         // printf("    count = %d\n", count);
 864     }
 865     fprintf(fd, "%d 2 1 0\n", count + excount);
 866
 867     // now write out extra node data
 868     for ( i = 1; i <= excount; i++ ) {
 869         fprintf(fd, "%d %.2f %.2f %.2f\n",
 870                 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
 871     }
 872
 873     // write out actual node data
 874     count = excount + 1;
 875     for ( j = rowmin; j <= rowmax; j++ ) {
 876         for ( i = colmin; i <= colmax; i++ ) {
 877             if ( output_data[i][j] > -9000.0 ) {
 878                 fprintf(fd, "%d %.2f %.2f %.2f\n",
 879                         count++,
 880                         originx + (double)i * col_step,
 881                         originy + (double)j * row_step,
 882                         output_data[i][j]);
 883             }
 884         }
 885         // printf("    count = %d\n", count);
 886     }
 887
 888     fclose(fd);
 889 }
 890 #endif
 891
 892
 893 FGDem::~FGDem( void ) {
 894     // printf("class FGDem DEstructor called.\n");
 895     delete [] dem_data;
 896     delete [] output_data;
 897 }
 898
 899
 900 // $Log$
 901 // Revision 1.25  1999/03/12 22:53:07  curt
 902 // Added a routine to dump out the portion of the dem data covered by a
 903 // specified bucket.  Other changes related to needs of scenery tools overhaul.
 904 //
 905 // Revision 1.24  1999/03/11 23:31:56  curt
 906 // Tweaks to use newbucket.hxx
 907 //
 908 // Revision 1.23  1999/03/10 01:09:12  curt
 909 // Tweaks to go along with scenery tools overhaul.
 910 // Added a new constructor that accepts the file name.
 911 //
 912 // Revision 1.22  1999/01/19 20:56:56  curt
 913 // MacOS portability changes contributed by "Robert Puyol" <puyol@abvent.fr>
 914 //
 915 // Revision 1.21  1998/11/06 14:04:32  curt
 916 // Changes due to updates in fgstream.
 917 //
 918 // Revision 1.20  1998/10/28 19:38:20  curt
 919 // Elliminate some unnecessary win32 specific stuff (by Norman Vine)
 920 //
 921 // Revision 1.19  1998/10/22 21:59:19  curt
 922 // Fixed a couple subtle bugs that resulted from some of my c++ conversions.
 923 // One bug could cause a segfault on certain input, and the other bug could
 924 // cause the whole procedure to go balistic and generate huge files (also only
 925 // on rare input combinations.)
 926 //
 927 // Revision 1.18  1998/10/18 01:17:09  curt
 928 // Point3D tweaks.
 929 //
 930 // Revision 1.17  1998/10/16 19:08:12  curt
 931 // Portability updates from Bernie Bright.
 932 //
 933 // Revision 1.16  1998/10/02 21:41:39  curt
 934 // Fixes for win32.
 935 //
 936 // Revision 1.15  1998/09/21 20:53:59  curt
 937 // minor tweaks to clean a few additional things up after the rewrite.
 938 //
 939 // Revision 1.14  1998/09/19 17:59:45  curt
 940 // Use c++ streams (fg_gzifstream).  Also converted many character arrays to
 941 // the string class.
 942 //
 943 // Revision 1.13  1998/09/09 16:24:04  curt
 944 // Fixed a bug in the handling of exclude files which was causing
 945 // a crash by calling fclose() on an invalid file handle.
 946 //
 947 // Revision 1.12  1998/08/24 20:03:31  curt
 948 // Eliminated a possible memory overrun error.
 949 // Use the proper free() rather than the incorrect delete().
 950 //
 951 // Revision 1.11  1998/07/20 12:46:11  curt
 952 // When outputing to a .node file, first check for an optional
 953 // "index.node.ex" file in case there is a .poly file to go along with this
 954 // node file.  Include these nodes first since they are referenced by position
 955 // from the .poly file.  This is my first pass at adding an area "cutout"
 956 // feature to the terrain generation pipeline.
 957 //
 958 // Revision 1.10  1998/07/13 20:58:02  curt
 959 // .
 960 //
 961 // Revision 1.9  1998/07/13 15:29:49  curt
 962 // Added #ifdef HAVE_CONFIG_H
 963 //
 964 // Revision 1.8  1998/07/04 00:47:18  curt
 965 // typedef'd struct fgBUCKET.
 966 //
 967 // Revision 1.7  1998/06/05 18:14:39  curt
 968 // Abort out early when reading the "A" record if it doesn't look like
 969 // a proper DEM file.
 970 //
 971 // Revision 1.6  1998/05/02 01:49:21  curt
 972 // Fixed a bug where the wrong variable was being initialized.
 973 //
 974 // Revision 1.5  1998/04/25 15:00:32  curt
 975 // Changed "r" to "rb" in gzopen() options.  This fixes bad behavior in win32.
 976 //
 977 // Revision 1.4  1998/04/22 13:14:46  curt
 978 // Fixed a bug in zlib usage.
 979 //
 980 // Revision 1.3  1998/04/18 03:53:05  curt
 981 // Added zlib support.
 982 //
 983 // Revision 1.2  1998/04/14 02:43:27  curt
 984 // Used "new" to auto-allocate large DEM parsing arrays in class constructor.
 985 //
 986 // Revision 1.1  1998/04/08 22:57:22  curt
 987 // Adopted Gnu automake/autoconf system.
 988 //
 989 // Revision 1.3  1998/04/06 21:09:41  curt
 990 // Additional win32 support.
 991 // Fixed a bad bug in dem file parsing that was causing the output to be
 992 // flipped about x = y.
 993 //
 994 // Revision 1.2  1998/03/23 20:35:41  curt
 995 // Updated to use FG_EPSILON
 996 //
 997 // Revision 1.1  1998/03/19 02:54:47  curt
 998 // Reorganized into a class lib called fgDEM.
 999 //
1000 // Revision 1.1  1998/03/19 01:46:28  curt
1001 // Initial revision.
1002 //