Array/array.cxx

   1 // array.cxx -- Array management class
   2 //
   3 // Written by Curtis Olson, started March 1998.
   4 //
   5 // Copyright (C) 1998 - 1999  Curtis L. Olson  - curt@flightgear.org
   6 //
   7 // This program is free software; you can redistribute it and/or
   8 // modify it under the terms of the GNU General Public License as
   9 // published by the Free Software Foundation; either version 2 of the
  10 // License, or (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful, but
  13 // WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 // General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20 //
  21 // $Id$
  22 // (Log is kept at end of this file)
  23
  24
  25 #ifdef HAVE_CONFIG_H
  26 #  include <config.h>
  27 #endif
  28
  29 #include <Include/compiler.h>
  30
  31 // #include <ctype.h>    // isspace()
  32 // #include <stdlib.h>   // atoi()
  33 // #include <math.h>     // rint()
  34 // #include <stdio.h>
  35 // #include <string.h>
  36 // #ifdef HAVE_SYS_STAT_H
  37 // #  include <sys/stat.h> // stat()
  38 // #endif
  39 // #ifdef FG_HAVE_STD_INCLUDES
  40 // #  include <cerrno>
  41 // #else
  42 // #  include <errno.h>
  43 // #endif
  44 // #ifdef HAVE_UNISTD_H
  45 // # include <unistd.h>   // stat()
  46 // #endif
  47
  48 #include STL_STRING
  49
  50 #include <Include/fg_constants.h>
  51 #include <Misc/fgstream.hxx>
  52 #include <Misc/strutils.hxx>
  53 #include <Math/leastsqs.hxx>
  54
  55 #include "array.hxx"
  56
  57 FG_USING_STD(string);
  58
  59
  60 FGArray::FGArray( void ) {
  61     // cout << "class FGArray CONstructor called." << endl;
  62     in_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
  63     // out_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
  64 }
  65
  66
  67 FGArray::FGArray( const string &file ) {
  68     // cout << "class FGArray CONstructor called." << endl;
  69     in_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
  70     // out_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
  71
  72     FGArray::open(file);
  73 }
  74
  75
  76 // open an Array file
  77 int
  78 FGArray::open( const string& file ) {
  79     // open input file (or read from stdin)
  80     if ( file ==  "-" ) {
  81         cout << "  Opening array data pipe from stdin" << endl;
  82         // fd = stdin;
  83         // fd = gzdopen(STDIN_FILENO, "r");
  84         cout << "  Not yet ported ..." << endl;
  85         return 0;
  86     } else {
  87         in = new fg_gzifstream( file );
  88         if ( !(*in) ) {
  89             cout << "  Cannot open " << file << endl;
  90             return 0;
  91         }
  92         cout << "  Opening array data file: " << file << endl;
  93     }
  94
  95     return 1;
  96 }
  97
  98
  99 // close an Array file
 100 int
 101 FGArray::close() {
 102     // the fg_gzifstream doesn't seem to have a close()
 103
 104     delete in;
 105
 106     return 1;
 107 }
 108
 109
 110 // parse Array file
 111 int
 112 FGArray::parse() {
 113     *in >> originx >> originy;
 114     *in >> cols >> col_step;
 115     *in >> rows >> row_step;
 116
 117     cout << "    origin  = " << originx << "  " << originy << endl;
 118     cout << "    cols = " << cols << "  rows = " << rows << endl;
 119     cout << "    col_step = " << col_step << "  row_step = " << row_step <<endl;
 120
 121     for ( int i = 0; i < cols; i++ ) {
 122         for ( int j = 0; j < rows; j++ ) {
 123             *in >> in_data[i][j];
 124         }
 125     }
 126
 127     cout << "    Done parsing\n";
 128
 129     return 1;
 130 }
 131
 132
 133 // add a node to the output (fitted) node list
 134 void FGArray::add_fit_node( int i, int j, double val ) {
 135     double x = (originx + i * col_step) / 3600.0;
 136     double y = (originy + j * row_step) / 3600.0;
 137     cout << Point3D(x, y, val) << endl;
 138     node_list.push_back( Point3D(x, y, val) );
 139 }
 140
 141
 142 #if 0
 143 // Initialize output mesh structure
 144 void FGArray::outputmesh_init( void ) {
 145     int i, j;
 146
 147     for ( j = 0; j < ARRAY_SIZE_1; j++ ) {
 148         for ( i = 0; i < ARRAY_SIZE_1; i++ ) {
 149             out_data[i][j] = -9999.0;
 150         }
 151     }
 152 }
 153
 154
 155 // Get the value of a mesh node
 156 double FGArray::outputmesh_get_pt( int i, int j ) {
 157     return ( out_data[i][j] );
 158 }
 159
 160
 161 // Set the value of a mesh node
 162 void FGArray::outputmesh_set_pt( int i, int j, double value ) {
 163     // cout << "Setting data[" << i << "][" << j << "] = " << value << endl;
 164    out_data[i][j] = value;
 165 }
 166 #endif
 167
 168
 169 // Use least squares to fit a simpler data set to dem data
 170 void FGArray::fit( double error ) {
 171     double x[ARRAY_SIZE_1], y[ARRAY_SIZE_1];
 172     double m, b, max_error, error_sq;
 173     double x1, y1;
 174     // double ave_error;
 175     double cury, lasty;
 176     int n, row, start, end;
 177     int colmin, colmax, rowmin, rowmax;
 178     bool good_fit;
 179     // FILE *dem, *fit, *fit1;
 180
 181     error_sq = error * error;
 182
 183     cout << "  Initializing output mesh structure" << endl;
 184     // outputmesh_init();
 185
 186     // determine dimensions
 187     colmin = 0;
 188     colmax = cols;
 189     rowmin = 0;
 190     rowmax = rows;
 191     cout << "  Fitting region = " << colmin << "," << rowmin << " to "
 192          << colmax << "," << rowmax << endl;;
 193
 194     // include the corners explicitly
 195     add_fit_node( colmin, rowmin, in_data[colmin][rowmin] );
 196     add_fit_node( colmin, rowmax-1, in_data[colmin][rowmax] );
 197     add_fit_node( colmax-1, rowmin, in_data[colmax][rowmin] );
 198     add_fit_node( colmax-1, rowmax-1, in_data[colmax][rowmax] );
 199
 200     cout << "  Beginning best fit procedure" << endl;
 201     lasty = 0;
 202
 203     for ( row = rowmin; row < rowmax; row++ ) {
 204         // fit  = fopen("fit.dat",  "w");
 205         // fit1 = fopen("fit1.dat", "w");
 206
 207         start = colmin;
 208
 209         // cout << "    fitting row = " << row << endl;
 210
 211         while ( start < colmax - 1 ) {
 212             end = start + 1;
 213             good_fit = true;
 214
 215             x[0] = start * col_step;
 216             y[0] = in_data[start][row];
 217
 218             x[1] = end * col_step;
 219             y[1] = in_data[end][row];
 220
 221             n = 2;
 222
 223             // cout << "Least square of first 2 points" << endl;
 224             least_squares(x, y, n, &m, &b);
 225
 226             end++;
 227
 228             while ( (end < colmax) && good_fit ) {
 229                 ++n;
 230                 // cout << "Least square of first " << n << " points" << endl;
 231                 x[n-1] = x1 = end * col_step;
 232                 y[n-1] = y1 = in_data[end][row];
 233                 least_squares_update(x1, y1, &m, &b);
 234                 // ave_error = least_squares_error(x, y, n, m, b);
 235                 max_error = least_squares_max_error(x, y, n, m, b);
 236
 237                 /*
 238                 printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n",
 239                 start, end, ave_error, max_error, m, b);
 240
 241                 f = fopen("gnuplot.dat", "w");
 242                 for ( j = 0; j <= end; j++) {
 243                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 244                             in_data[row][j]);
 245                 }
 246                 for ( j = start; j <= end; j++) {
 247                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 248                             in_data[row][j]);
 249                 }
 250                 fclose(f);
 251
 252                 printf("Please hit return: "); gets(junk);
 253                 */
 254
 255                 if ( max_error > error_sq ) {
 256                     good_fit = false;
 257                 }
 258
 259                 end++;
 260             }
 261
 262             if ( !good_fit ) {
 263                 // error exceeded the threshold, back up
 264                 end -= 2;  // back "end" up to the last good enough fit
 265                 n--;       // back "n" up appropriately too
 266             } else {
 267                 // we popped out of the above loop while still within
 268                 // the error threshold, so we must be at the end of
 269                 // the data set
 270                 end--;
 271             }
 272
 273             least_squares(x, y, n, &m, &b);
 274             // ave_error = least_squares_error(x, y, n, m, b);
 275             max_error = least_squares_max_error(x, y, n, m, b);
 276
 277             /*
 278             printf("\n");
 279             printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n",
 280                    start, end, ave_error, max_error, m, b);
 281             printf("\n");
 282
 283             fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
 284             fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
 285             */
 286
 287             if ( start > colmin ) {
 288                 // skip this for the first line segment
 289                 cury = m * x[0] + b;
 290                 add_fit_node( start, row, (lasty + cury) / 2 );
 291                 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
 292             }
 293
 294             lasty = m * x[end-start] + b;
 295             start = end;
 296         }
 297
 298         /*
 299         fclose(fit);
 300         fclose(fit1);
 301
 302         dem = fopen("gnuplot.dat", "w");
 303         for ( j = 0; j < ARRAY_SIZE_1; j++) {
 304             fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
 305                     in_data[j][row]);
 306         }
 307         fclose(dem);
 308         */
 309
 310         // NOTICE, this is for testing only.  This instance of
 311         // output_nodes should be removed.  It should be called only
 312         // once at the end once all the nodes have been generated.
 313         // newmesh_output_nodes(&nm, "mesh.node");
 314         // printf("Please hit return: "); gets(junk);
 315     }
 316
 317     // outputmesh_output_nodes(fg_root, p);
 318 }
 319
 320
 321 // return the current altitude based on grid data.  We should rewrite
 322 // this to interpolate exact values, but for now this is good enough
 323 double FGArray::interpolate_altitude( double lon, double lat ) {
 324     // we expect incoming (lon,lat) to be in arcsec for now
 325
 326     double xlocal, ylocal, dx, dy, zA, zB, elev;
 327     int x1, x2, x3, y1, y2, y3;
 328     float z1, z2, z3;
 329     int xindex, yindex;
 330
 331     /* determine if we are in the lower triangle or the upper triangle
 332        ______
 333        |   /|
 334        |  / |
 335        | /  |
 336        |/   |
 337        ------
 338
 339        then calculate our end points
 340      */
 341
 342     xlocal = (lon - originx) / col_step;
 343     ylocal = (lat - originy) / row_step;
 344
 345     xindex = (int)(xlocal);
 346     yindex = (int)(ylocal);
 347
 348     // printf("xindex = %d  yindex = %d\n", xindex, yindex);
 349
 350     if ( xindex + 1 == cols ) {
 351         xindex--;
 352     }
 353
 354     if ( yindex + 1 == rows ) {
 355         yindex--;
 356     }
 357
 358     if ( (xindex < 0) || (xindex + 1 >= cols) ||
 359          (yindex < 0) || (yindex + 1 >= rows) ) {
 360         cout << "WARNING: Attempt to interpolate value outside of array!!!"
 361              << endl;
 362         return(-9999);
 363     }
 364
 365     dx = xlocal - xindex;
 366     dy = ylocal - yindex;
 367
 368     if ( dx > dy ) {
 369         // lower triangle
 370         // printf("  Lower triangle\n");
 371
 372         x1 = xindex;
 373         y1 = yindex;
 374         z1 = in_data[x1][y1];
 375
 376         x2 = xindex + 1;
 377         y2 = yindex;
 378         z2 = in_data[x2][y2];
 379
 380         x3 = xindex + 1;
 381         y3 = yindex + 1;
 382         z3 = in_data[x3][y3];
 383
 384         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
 385         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
 386         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
 387         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
 388
 389         zA = dx * (z2 - z1) + z1;
 390         zB = dx * (z3 - z1) + z1;
 391
 392         // printf("  zA = %.2f  zB = %.2f\n", zA, zB);
 393
 394         if ( dx > FG_EPSILON ) {
 395             elev = dy * (zB - zA) / dx + zA;
 396         } else {
 397             elev = zA;
 398         }
 399     } else {
 400         // upper triangle
 401         // printf("  Upper triangle\n");
 402
 403         x1 = xindex;
 404         y1 = yindex;
 405         z1 = in_data[x1][y1];
 406
 407         x2 = xindex;
 408         y2 = yindex + 1;
 409         z2 = in_data[x2][y2];
 410
 411         x3 = xindex + 1;
 412         y3 = yindex + 1;
 413         z3 = in_data[x3][y3];
 414
 415         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
 416         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
 417         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
 418         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
 419
 420         zA = dy * (z2 - z1) + z1;
 421         zB = dy * (z3 - z1) + z1;
 422
 423         // printf("  zA = %.2f  zB = %.2f\n", zA, zB );
 424         // printf("  xB - xA = %.2f\n", col_step * dy / row_step);
 425
 426         if ( dy > FG_EPSILON ) {
 427             elev = dx * (zB - zA) / dy    + zA;
 428         } else {
 429             elev = zA;
 430         }
 431     }
 432
 433     return(elev);
 434 }
 435
 436
 437 #if 0
 438 // Write out a node file that can be used by the "triangle" program.
 439 // Check for an optional "index.node.ex" file in case there is a .poly
 440 // file to go along with this node file.  Include these nodes first
 441 // since they are referenced by position from the .poly file.
 442 void FGArray::outputmesh_output_nodes( const string& fg_root, FGBucket& p )
 443 {
 444     double exnodes[MAX_EX_NODES][3];
 445     struct stat stat_buf;
 446     string dir;
 447     char file[256], exfile[256];
 448 #ifdef WIN32
 449     char tmp_path[256];
 450 #endif
 451     string command;
 452     FILE *fd;
 453     long int index;
 454     int colmin, colmax, rowmin, rowmax;
 455     int i, j, count, excount, result;
 456
 457     // determine dimensions
 458     colmin = p.get_x() * ( (cols - 1) / 8);
 459     colmax = colmin + ( (cols - 1) / 8);
 460     rowmin = p.get_y() * ( (rows - 1) / 8);
 461     rowmax = rowmin + ( (rows - 1) / 8);
 462     cout << "  dumping region = " << colmin << "," << rowmin << " to " <<
 463         colmax << "," << rowmax << "\n";
 464
 465     // generate the base directory
 466     string base_path = p.gen_base_path();
 467     cout << "  fg_root = " << fg_root << "  Base Path = " << base_path << endl;
 468     dir = fg_root + "/Scenery/" + base_path;
 469     cout << "  Dir = " << dir << endl;
 470
 471     // stat() directory and create if needed
 472     errno = 0;
 473     result = stat(dir.c_str(), &stat_buf);
 474     if ( result != 0 && errno == ENOENT ) {
 475         cout << "  Creating directory\n";
 476
 477         command = "mkdir -p " + dir + "\n";
 478         system( command.c_str() );
 479     } else {
 480         // assume directory exists
 481     }
 482
 483     // get index and generate output file name
 484     index = p.gen_index();
 485     sprintf(file, "%s/%ld.node", dir.c_str(), index);
 486
 487     // get (optional) extra node file name (in case there is matching
 488     // .poly file.
 489     strcpy(exfile, file);
 490     strcat(exfile, ".ex");
 491
 492     // load extra nodes if they exist
 493     excount = 0;
 494     if ( (fd = fopen(exfile, "r")) != NULL ) {
 495         int junki;
 496         fscanf(fd, "%d %d %d %d", &excount, &junki, &junki, &junki);
 497
 498         if ( excount > MAX_EX_NODES - 1 ) {
 499             printf("Error, too many 'extra' nodes, increase array size\n");
 500             exit(-1);
 501         } else {
 502             printf("    Expecting %d 'extra' nodes\n", excount);
 503         }
 504
 505         for ( i = 1; i <= excount; i++ ) {
 506             fscanf(fd, "%d %lf %lf %lf\n", &junki,
 507                    &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]);
 508             printf("(extra) %d %.2f %.2f %.2f\n",
 509                     i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
 510         }
 511         fclose(fd);
 512     }
 513
 514     printf("Creating node file:  %s\n", file);
 515     fd = fopen(file, "w");
 516
 517     // first count regular nodes to generate header
 518     count = 0;
 519     for ( j = rowmin; j <= rowmax; j++ ) {
 520         for ( i = colmin; i <= colmax; i++ ) {
 521             if ( out_data[i][j] > -9000.0 ) {
 522                 count++;
 523             }
 524         }
 525         // printf("    count = %d\n", count);
 526     }
 527     fprintf(fd, "%d 2 1 0\n", count + excount);
 528
 529     // now write out extra node data
 530     for ( i = 1; i <= excount; i++ ) {
 531         fprintf(fd, "%d %.2f %.2f %.2f\n",
 532                 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
 533     }
 534
 535     // write out actual node data
 536     count = excount + 1;
 537     for ( j = rowmin; j <= rowmax; j++ ) {
 538         for ( i = colmin; i <= colmax; i++ ) {
 539             if ( out_data[i][j] > -9000.0 ) {
 540                 fprintf(fd, "%d %.2f %.2f %.2f\n",
 541                         count++,
 542                         originx + (double)i * col_step,
 543                         originy + (double)j * row_step,
 544                         out_data[i][j]);
 545             }
 546         }
 547         // printf("    count = %d\n", count);
 548     }
 549
 550     fclose(fd);
 551 }
 552 #endif
 553
 554
 555 FGArray::~FGArray( void ) {
 556     // printf("class FGArray DEstructor called.\n");
 557     delete [] in_data;
 558     // delete [] out_data;
 559 }
 560
 561
 562 // $Log$
 563 // Revision 1.4  1999/03/20 20:32:51  curt
 564 // First mostly successful tile triangulation works.  There's plenty of tweaking
 565 // to do, but we are marching in the right direction.
 566 //
 567 // Revision 1.3  1999/03/17 23:48:17  curt
 568 // Removed forced -g compile flag.
 569 // Fixed a couple compiler warnings.
 570 //
 571 // Revision 1.2  1999/03/13 23:50:26  curt
 572 // Tweaked output formatting a bit.
 573 //
 574 // Revision 1.1  1999/03/13 18:45:02  curt
 575 // Initial revision. (derived from libDEM.a code.)
 576 //