diff --git a/src/Lib/DEM/dem.cxx b/src/Lib/DEM/dem.cxx index 67b328b9..afee8cd7 100644 --- a/src/Lib/DEM/dem.cxx +++ b/src/Lib/DEM/dem.cxx @@ -61,21 +61,7 @@ SG_USING_STD(cout); SG_USING_STD(endl); -#define MAX_EX_NODES 10000 - -#if 0 -#ifdef WIN32 -# ifdef __BORLANDC__ -# include <dir.h> -# define MKDIR(a) mkdir(a) -# else -# define MKDIR(a) mkdir(a,S_IRWXU) // I am just guessing at this flag (NHV) -# endif // __BORLANDC__ -#endif // WIN32 -#endif //0 - - -FGDem::FGDem( void ) : +FGDem::FGDem() : z_units(2) // meters { // cout << "class FGDem CONstructor called." << endl; @@ -476,406 +462,7 @@ FGDem::write_area( const string& root, SGBucket& b, bool compress ) { } -#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 ) { - // we expect incoming (lon,lat) to be in arcsec for now - - double xlocal, ylocal, dx, dy, zA, zB, elev; - int x1, x2, x3, y1, y2, y3; - float z1, z2, z3; - int xindex, yindex; - - /* determine if we are in the lower triangle or the upper triangle - ______ - | /| - | / | - | / | - |/ | - ------ - - then calculate our end points - */ - - xlocal = (lon - originx) / col_step; - ylocal = (lat - originy) / row_step; - - xindex = (int)(xlocal); - yindex = (int)(ylocal); - - // printf("xindex = %d yindex = %d\n", xindex, yindex); - - if ( xindex + 1 == cols ) { - xindex--; - } - - if ( yindex + 1 == rows ) { - yindex--; - } - - if ( (xindex < 0) || (xindex + 1 >= cols) || - (yindex < 0) || (yindex + 1 >= rows) ) { - return(-9999); - } - - dx = xlocal - xindex; - dy = ylocal - yindex; - - if ( dx > dy ) { - // lower triangle - // printf(" Lower triangle\n"); - - x1 = xindex; - y1 = yindex; - z1 = dem_data[x1][y1]; - - x2 = xindex + 1; - y2 = yindex; - z2 = dem_data[x2][y2]; - - x3 = xindex + 1; - y3 = yindex + 1; - z3 = dem_data[x3][y3]; - - // printf(" dx = %.2f dy = %.2f\n", dx, dy); - // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1); - // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); - // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); - - zA = dx * (z2 - z1) + z1; - zB = dx * (z3 - z1) + z1; - - // printf(" zA = %.2f zB = %.2f\n", zA, zB); - - if ( dx > SG_EPSILON ) { - elev = dy * (zB - zA) / dx + zA; - } else { - elev = zA; - } - } else { - // upper triangle - // printf(" Upper triangle\n"); - - x1 = xindex; - y1 = yindex; - z1 = dem_data[x1][y1]; - - x2 = xindex; - y2 = yindex + 1; - z2 = dem_data[x2][y2]; - - x3 = xindex + 1; - y3 = yindex + 1; - z3 = dem_data[x3][y3]; - - // printf(" dx = %.2f dy = %.2f\n", dx, dy); - // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1); - // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2); - // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3); - - zA = dy * (z2 - z1) + z1; - zB = dy * (z3 - z1) + z1; - - // printf(" zA = %.2f zB = %.2f\n", zA, zB ); - // printf(" xB - xA = %.2f\n", col_step * dy / row_step); - - if ( dy > SG_EPSILON ) { - elev = dx * (zB - zA) / dy + zA; - } else { - elev = zA; - } - } - - return(elev); -} - - -// Use least squares to fit a simpler data set to dem data -void FGDem::fit( double error, SGBucket& p ) { - double x[DEM_SIZE_1], y[DEM_SIZE_1]; - double m, b, ave_error, max_error; - double cury, lasty; - int n, row, start, end; - int colmin, colmax, rowmin, rowmax; - bool good_fit; - // FILE *dem, *fit, *fit1; - - printf("Initializing output mesh structure\n"); - outputmesh_init(); - - // determine dimensions - colmin = p.get_x() * ( (cols - 1) / 8); - colmax = colmin + ( (cols - 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); - - // include the corners explicitly - outputmesh_set_pt(colmin, rowmin, dem_data[colmin][rowmin]); - outputmesh_set_pt(colmin, rowmax, dem_data[colmin][rowmax]); - outputmesh_set_pt(colmax, rowmax, dem_data[colmax][rowmax]); - outputmesh_set_pt(colmax, rowmin, dem_data[colmax][rowmin]); - - printf("Beginning best fit procedure\n"); - - for ( row = rowmin; row <= rowmax; row++ ) { - // fit = fopen("fit.dat", "w"); - // fit1 = fopen("fit1.dat", "w"); - - start = colmin; - - // printf(" fitting row = %d\n", row); - - while ( start < colmax ) { - end = start + 1; - good_fit = true; - - x[(end - start) - 1] = 0.0 + ( start * col_step ); - y[(end - start) - 1] = dem_data[start][row]; - - while ( (end <= colmax) && good_fit ) { - n = (end - start) + 1; - // printf("Least square of first %d points\n", n); - x[end - start] = 0.0 + ( end * col_step ); - y[end - start] = dem_data[end][row]; - least_squares(x, y, n, &m, &b); - ave_error = least_squares_error(x, y, n, m, b); - max_error = least_squares_max_error(x, y, n, m, b); - - /* - printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n", - start, end, ave_error, max_error, m, b); - - f = fopen("gnuplot.dat", "w"); - for ( j = 0; j <= end; j++) { - fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ), - dem_data[row][j]); - } - for ( j = start; j <= end; j++) { - fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ), - dem_data[row][j]); - } - fclose(f); - - printf("Please hit return: "); gets(junk); - */ - - if ( max_error > error ) { - good_fit = false; - } - - end++; - } - - if ( !good_fit ) { - // error exceeded the threshold, back up - end -= 2; // back "end" up to the last good enough fit - n--; // back "n" up appropriately too - } else { - // we popped out of the above loop while still within - // the error threshold, so we must be at the end of - // the data set - end--; - } - - least_squares(x, y, n, &m, &b); - ave_error = least_squares_error(x, y, n, m, b); - max_error = least_squares_max_error(x, y, n, m, b); - - /* - printf("\n"); - printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n", - start, end, ave_error, max_error, m, b); - printf("\n"); - - fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b); - fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b); - */ - - if ( start > colmin ) { - // skip this for the first line segment - cury = m * x[0] + b; - outputmesh_set_pt(start, row, (lasty + cury) / 2); - // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2); - } - - lasty = m * x[end-start] + b; - start = end; - } - - /* - fclose(fit); - fclose(fit1); - - dem = fopen("gnuplot.dat", "w"); - for ( j = 0; j < DEM_SIZE_1; j++) { - fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ), - dem_data[j][row]); - } - fclose(dem); - */ - - // NOTICE, this is for testing only. This instance of - // output_nodes should be removed. It should be called only - // once at the end once all the nodes have been generated. - // newmesh_output_nodes(&nm, "mesh.node"); - // printf("Please hit return: "); gets(junk); - } - - // outputmesh_output_nodes(fg_root, p); -} - - -// Initialize output mesh structure -void FGDem::outputmesh_init( void ) { - int i, j; - - for ( j = 0; j < DEM_SIZE_1; j++ ) { - for ( i = 0; i < DEM_SIZE_1; i++ ) { - output_data[i][j] = -9999.0; - } - } -} - - -// Get the value of a mesh node -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 ) { - // printf("Setting data[%d][%d] = %.2f\n", i, j, value); - output_data[i][j] = value; -} - - -// Write out a node file that can be used by the "triangle" program. -// 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, SGBucket& p ) -{ - double exnodes[MAX_EX_NODES][3]; - struct stat stat_buf; - string dir; - char file[256], exfile[256]; -#ifdef WIN32 - char tmp_path[256]; -#endif - string command; - FILE *fd; - long int index; - int colmin, colmax, rowmin, rowmax; - int i, j, count, excount, result; - - // determine dimensions - colmin = p.get_x() * ( (cols - 1) / 8); - colmax = colmin + ( (cols - 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 - string base_path = p.gen_base_path(); - cout << "fg_root = " << fg_root << " Base Path = " << base_path << endl; - dir = fg_root + "/" + base_path; - cout << "Dir = " << dir << endl; - - // stat() directory and create if needed - errno = 0; - result = stat(dir.c_str(), &stat_buf); - if ( result != 0 && errno == ENOENT ) { - cout << "Creating directory\n"; - -#ifdef _MSC_VER - fg_mkdir( dir.c_str() ); -#else - command = "mkdir -p " + dir + "\n"; - system( command.c_str() ); -#endif - - } else { - // assume directory exists - } - - // get index and generate output file name - index = p.gen_index(); - sprintf(file, "%s/%ld.node", dir.c_str(), index); - - // get (optional) extra node file name (in case there is matching - // .poly file. - strcpy(exfile, file); - strcat(exfile, ".ex"); - - // load extra nodes if they exist - excount = 0; - if ( (fd = fopen(exfile, "r")) != NULL ) { - int junki; - fscanf(fd, "%d %d %d %d", &excount, &junki, &junki, &junki); - - if ( excount > MAX_EX_NODES - 1 ) { - printf("Error, too many 'extra' nodes, increase array size\n"); - exit(-1); - } else { - printf(" Expecting %d 'extra' nodes\n", excount); - } - - for ( i = 1; i <= excount; i++ ) { - fscanf(fd, "%d %lf %lf %lf\n", &junki, - &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]); - printf("(extra) %d %.2f %.2f %.2f\n", - i, exnodes[i][0], exnodes[i][1], exnodes[i][2]); - } - fclose(fd); - } - - printf("Creating node file: %s\n", file); - fd = fopen(file, "w"); - - // first count regular nodes to generate header - count = 0; - for ( j = rowmin; j <= rowmax; j++ ) { - for ( i = colmin; i <= colmax; i++ ) { - if ( output_data[i][j] > -9000.0 ) { - count++; - } - } - // printf(" count = %d\n", count); - } - fprintf(fd, "%d 2 1 0\n", count + excount); - - // now write out extra node data - for ( i = 1; i <= excount; i++ ) { - fprintf(fd, "%d %.2f %.2f %.2f\n", - i, exnodes[i][0], exnodes[i][1], exnodes[i][2]); - } - - // write out actual node data - count = excount + 1; - for ( j = rowmin; j <= rowmax; j++ ) { - for ( i = colmin; i <= colmax; i++ ) { - if ( output_data[i][j] > -9000.0 ) { - fprintf(fd, "%d %.2f %.2f %.2f\n", - count++, - originx + (double)i * col_step, - originy + (double)j * row_step, - output_data[i][j]); - } - } - // printf(" count = %d\n", count); - } - - fclose(fd); -} -#endif - - -FGDem::~FGDem( void ) { +FGDem::~FGDem() { // printf("class FGDem DEstructor called.\n"); delete [] dem_data; delete [] output_data; diff --git a/src/Lib/DEM/dem.hxx b/src/Lib/DEM/dem.hxx index 28d9cb86..7fe63c9b 100644 --- a/src/Lib/DEM/dem.hxx +++ b/src/Lib/DEM/dem.hxx @@ -96,11 +96,11 @@ private: public: // Constructor - FGDem( void ); + FGDem(); FGDem( const string& file ); // Destructor - ~FGDem( void ); + ~FGDem(); // open a DEM file (use "-" if input is coming from stdin) int open ( const string& file ); @@ -122,28 +122,6 @@ public: // hand corner. int write_area( const string& root, SGBucket& b, bool compress ); -#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 interpolate_altitude( double lon, double lat ); - - // Use least squares to fit a simpler data set to dem data - void fit( double error, SGBucket& p ); - - // Initialize output mesh structure - void outputmesh_init( void ); - - // Get the value of a mesh node - double outputmesh_get_pt( int i, int j ); - - // Set the value of a mesh node - void outputmesh_set_pt( int i, int j, double value ); - - // Write out a node file that can be used by the "triangle" program - void outputmesh_output_nodes( const string& fg_root, SGBucket& p ); -#endif - // Informational methods inline double get_originx() const { return originx; } inline double get_originy() const { return originy; }