961 lines
24 KiB
C++
961 lines
24 KiB
C++
// -*- Mode: C++ -*-
|
|
//
|
|
// dem.c -- DEM management class
|
|
//
|
|
// Written by Curtis Olson, started March 1998.
|
|
//
|
|
// Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu
|
|
//
|
|
// This program is free software; you can redistribute it and/or
|
|
// modify it under the terms of the GNU General Public License as
|
|
// published by the Free Software Foundation; either version 2 of the
|
|
// License, or (at your option) any later version.
|
|
//
|
|
// This program 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
|
|
// 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., 675 Mass Ave, Cambridge, MA 02139, USA.
|
|
//
|
|
// $Id$
|
|
// (Log is kept at end of this file)
|
|
|
|
|
|
#ifdef HAVE_CONFIG_H
|
|
# include <config.h>
|
|
#endif
|
|
|
|
#include <Include/compiler.h>
|
|
|
|
#include <ctype.h> // isspace()
|
|
#include <stdlib.h> // atoi()
|
|
#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 "dem.hxx"
|
|
#include "leastsqs.hxx"
|
|
|
|
#include <Include/fg_constants.h>
|
|
|
|
|
|
#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 ) {
|
|
// printf("class fgDEM CONstructor called.\n");
|
|
dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
|
|
output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
|
|
}
|
|
|
|
|
|
#if 0
|
|
#ifdef WIN32
|
|
|
|
// return the file path name ( foo/bar/file.ext = foo/bar )
|
|
static void extract_path ( const char *in, char *base) {
|
|
int len, i;
|
|
|
|
len = strlen (in);
|
|
strcpy (base, in);
|
|
|
|
i = len - 1;
|
|
while ( (i >= 0) && (in[i] != '/') ) {
|
|
i--;
|
|
}
|
|
|
|
base[i] = '\0';
|
|
}
|
|
|
|
|
|
// Make a subdirectory
|
|
static int my_mkdir (const char *dir) {
|
|
struct stat stat_buf;
|
|
int result;
|
|
|
|
printf ("mk_dir() ");
|
|
|
|
result = stat (dir, &stat_buf);
|
|
|
|
if (result != 0) {
|
|
MKDIR (dir);
|
|
result = stat (dir, &stat_buf);
|
|
if (result != 0) {
|
|
printf ("problem creating %s\n", dir);
|
|
} else {
|
|
printf ("%s created\n", dir);
|
|
}
|
|
} else {
|
|
printf ("%s already exists\n", dir);
|
|
}
|
|
|
|
return (result);
|
|
}
|
|
|
|
#endif // WIN32
|
|
#endif //0
|
|
|
|
|
|
// open a DEM file
|
|
int fgDEM::open ( const string& file ) {
|
|
// open input file (or read from stdin)
|
|
if ( file == "-" ) {
|
|
printf("Loading DEM data file: stdin\n");
|
|
// fd = stdin;
|
|
// fd = gzdopen(STDIN_FILENO, "r");
|
|
printf("Not yet ported ...\n");
|
|
return 0;
|
|
} else {
|
|
in = new fg_gzifstream( file );
|
|
if ( !(*in) ) {
|
|
cout << "Cannot open " << file << endl;
|
|
return 0;
|
|
}
|
|
cout << "Loading DEM data file: " << file << endl;
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
// close a DEM file
|
|
int fgDEM::close () {
|
|
// the fg_gzifstream doesn't seem to have a close()
|
|
|
|
delete in;
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
// return next token from input stream
|
|
string fgDEM::next_token() {
|
|
string token;
|
|
|
|
*in >> token;
|
|
|
|
// cout << " returning " + token + "\n";
|
|
|
|
return token;
|
|
}
|
|
|
|
|
|
// return next integer from input stream
|
|
int fgDEM::next_int() {
|
|
int result;
|
|
|
|
*in >> result;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
// return next double from input stream
|
|
double fgDEM::next_double() {
|
|
double result;
|
|
|
|
*in >> result;
|
|
|
|
return result;
|
|
}
|
|
|
|
|
|
// return next exponential num from input stream
|
|
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;
|
|
|
|
for ( ; *p != 0; ++p )
|
|
{
|
|
if ( *p == 'D' )
|
|
*bp++ = 'E';
|
|
else
|
|
*bp++ = *p;
|
|
}
|
|
*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 i, inum;
|
|
double dnum;
|
|
string name, token;
|
|
char c;
|
|
char *ptr;
|
|
|
|
// get the name field (144 characters)
|
|
for ( i = 0; i < 144; i++ ) {
|
|
in->get(c);
|
|
name += c;
|
|
}
|
|
|
|
// clean off the trailing whitespace
|
|
name = trim(name);
|
|
cout << " Quad name field: " << name << endl;
|
|
|
|
// DEM level code, 3 reflects processing by DMA
|
|
inum = next_int();
|
|
cout << " DEM level code = " << inum << "\n";
|
|
|
|
if ( inum > 3 ) {
|
|
return 0;
|
|
}
|
|
|
|
// Pattern code, 1 indicates a regular elevation pattern
|
|
inum = next_int();
|
|
cout << " Pattern code = " << inum << "\n";
|
|
|
|
// Planimetric reference system code, 0 indicates geographic
|
|
// coordinate system.
|
|
inum = next_int();
|
|
cout << " Planimetric reference code = " << inum << "\n";
|
|
|
|
// Zone code
|
|
inum = next_int();
|
|
cout << " Zone code = " << inum << "\n";
|
|
|
|
// Map projection parameters (ignored)
|
|
for ( i = 0; i < 15; i++ ) {
|
|
dnum = next_exp();
|
|
// printf("%d: %f\n",i,dnum);
|
|
}
|
|
|
|
// Units code, 3 represents arc-seconds as the unit of measure for
|
|
// ground planimetric coordinates throughout the file.
|
|
inum = next_int();
|
|
if ( inum != 3 ) {
|
|
cout << " Unknown (X,Y) units code = " << inum << "!\n";
|
|
exit(-1);
|
|
}
|
|
|
|
// Units code; 2 represents meters as the unit of measure for
|
|
// elevation coordinates throughout the file.
|
|
inum = next_int();
|
|
if ( inum != 2 ) {
|
|
cout << " Unknown (Z) units code = " << inum << "!\n";
|
|
exit(-1);
|
|
}
|
|
|
|
// Number (n) of sides in the polygon which defines the coverage of
|
|
// the DEM file (usually equal to 4).
|
|
inum = next_int();
|
|
if ( inum != 4 ) {
|
|
cout << " Unknown polygon dimension = " << inum << "!\n";
|
|
exit(-1);
|
|
}
|
|
|
|
// Ground coordinates of bounding box in arc-seconds
|
|
dem_x1 = originx = next_exp();
|
|
dem_y1 = originy = next_exp();
|
|
cout << " Origin = (" << originx << "," << originy << ")\n";
|
|
|
|
dem_x2 = next_exp();
|
|
dem_y2 = next_exp();
|
|
|
|
dem_x3 = next_exp();
|
|
dem_y3 = next_exp();
|
|
|
|
dem_x4 = next_exp();
|
|
dem_y4 = next_exp();
|
|
|
|
// Minimum/maximum elevations in meters
|
|
dem_z1 = next_exp();
|
|
dem_z2 = next_exp();
|
|
cout << " Elevation range " << dem_z1 << " to " << dem_z2 << "\n";
|
|
|
|
// Counterclockwise angle from the primary axis of ground
|
|
// planimetric referenced to the primary axis of the DEM local
|
|
// reference system.
|
|
token = next_token();
|
|
|
|
// Accuracy code; 0 indicates that a record of accuracy does not
|
|
// exist and that no record type C will follow.
|
|
|
|
// DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1)
|
|
// depending on latitude
|
|
|
|
// I will eventually have to do something with this for data at
|
|
// higher latitudes */
|
|
token = next_token();
|
|
cout << " accuracy & spacial resolution string = " << token << endl;
|
|
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
|
|
cout << " Accuracy code = " << inum << "\n";
|
|
|
|
cout << " column step = " << col_step <<
|
|
" row step = " << row_step << "\n";
|
|
|
|
// dimension of arrays to follow (1)
|
|
token = next_token();
|
|
|
|
// number of profiles
|
|
dem_num_profiles = cols = next_int();
|
|
cout << " Expecting " << dem_num_profiles << " profiles\n";
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
// read and parse DEM "B" record
|
|
void fgDEM::read_b_record( ) {
|
|
string token;
|
|
int i;
|
|
|
|
// row / column id of this profile
|
|
prof_row = next_int();
|
|
prof_col = next_int();
|
|
// printf("col id = %d row id = %d\n", prof_col, prof_row);
|
|
|
|
// Number of columns and rows (elevations) in this profile
|
|
prof_num_rows = rows = next_int();
|
|
prof_num_cols = next_int();
|
|
// printf(" profile num rows = %d\n", prof_num_rows);
|
|
|
|
// Ground planimetric coordinates (arc-seconds) of the first
|
|
// elevation in the profile
|
|
prof_x1 = next_exp();
|
|
prof_y1 = next_exp();
|
|
// printf(" Starting at %.2f %.2f\n", prof_x1, prof_y1);
|
|
|
|
// Elevation of local datum for the profile. Always zero for
|
|
// 1-degree DEM, the reference is mean sea level.
|
|
token = next_token();
|
|
|
|
// Minimum and maximum elevations for the profile.
|
|
token = next_token();
|
|
token = next_token();
|
|
|
|
// One (usually) dimensional array (prof_num_cols,1) of elevations
|
|
for ( i = 0; i < prof_num_rows; i++ ) {
|
|
prof_data = next_int();
|
|
dem_data[cur_col][i] = (float)prof_data;
|
|
}
|
|
}
|
|
|
|
|
|
// parse dem file
|
|
int fgDEM::parse( ) {
|
|
int i;
|
|
|
|
cur_col = 0;
|
|
|
|
if ( !read_a_record() ) {
|
|
return(0);
|
|
}
|
|
|
|
for ( i = 0; i < dem_num_profiles; i++ ) {
|
|
// printf("Ready to read next b record\n");
|
|
read_b_record();
|
|
cur_col++;
|
|
|
|
if ( cur_col % 100 == 0 ) {
|
|
cout << " loaded " << cur_col << " profiles of data\n";
|
|
}
|
|
}
|
|
|
|
cout << " Done parsing\n";
|
|
|
|
return 1;
|
|
}
|
|
|
|
|
|
// return the current altitude based on mesh 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 > FG_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 > FG_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, fgBUCKET *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->x * ( (cols - 1) / 8);
|
|
colmax = colmin + ( (cols - 1) / 8);
|
|
rowmin = p->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, fgBUCKET *p ) {
|
|
double exnodes[MAX_EX_NODES][3];
|
|
struct stat stat_buf;
|
|
string dir;
|
|
char base_path[256], 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->x * ( (cols - 1) / 8);
|
|
colmax = colmin + ( (cols - 1) / 8);
|
|
rowmin = p->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);
|
|
cout << "fg_root = " << fg_root << " Base Path = " << base_path << endl;
|
|
dir = fg_root + "/Scenery/" + 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";
|
|
|
|
// #ifndef WIN32
|
|
|
|
command = "mkdir -p " + dir + "\n";
|
|
system( command.c_str() );
|
|
|
|
#if 0
|
|
// #else // WIN32
|
|
|
|
// Cygwin crashes when trying to output to node file
|
|
// explicitly making directory structure seems OK on Win95
|
|
|
|
extract_path (base_path, tmp_path);
|
|
|
|
dir = fg_root + "/Scenery";
|
|
if (my_mkdir ( dir.c_str() )) { exit (-1); }
|
|
|
|
dir = fg_root + "/Scenery/" + tmp_path;
|
|
if (my_mkdir ( dir.c_str() )) { exit (-1); }
|
|
|
|
dir = fg_root + "/Scenery/" + base_path;
|
|
if (my_mkdir ( dir.c_str() )) { exit (-1); }
|
|
|
|
// #endif // WIN32
|
|
#endif //0
|
|
|
|
} else {
|
|
// assume directory exists
|
|
}
|
|
|
|
// get index and generate output file name
|
|
index = fgBucketGenIndex(p);
|
|
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);
|
|
}
|
|
|
|
|
|
fgDEM::~fgDEM( void ) {
|
|
// printf("class fgDEM DEstructor called.\n");
|
|
delete [] dem_data;
|
|
delete [] output_data;
|
|
}
|
|
|
|
|
|
// $Log$
|
|
// Revision 1.22 1999/01/19 20:56:56 curt
|
|
// MacOS portability changes contributed by "Robert Puyol" <puyol@abvent.fr>
|
|
//
|
|
// Revision 1.21 1998/11/06 14:04:32 curt
|
|
// Changes due to updates in fgstream.
|
|
//
|
|
// Revision 1.20 1998/10/28 19:38:20 curt
|
|
// Elliminate some unnecessary win32 specific stuff (by Norman Vine)
|
|
//
|
|
// Revision 1.19 1998/10/22 21:59:19 curt
|
|
// Fixed a couple subtle bugs that resulted from some of my c++ conversions.
|
|
// One bug could cause a segfault on certain input, and the other bug could
|
|
// cause the whole procedure to go balistic and generate huge files (also only
|
|
// on rare input combinations.)
|
|
//
|
|
// Revision 1.18 1998/10/18 01:17:09 curt
|
|
// Point3D tweaks.
|
|
//
|
|
// Revision 1.17 1998/10/16 19:08:12 curt
|
|
// Portability updates from Bernie Bright.
|
|
//
|
|
// Revision 1.16 1998/10/02 21:41:39 curt
|
|
// Fixes for win32.
|
|
//
|
|
// Revision 1.15 1998/09/21 20:53:59 curt
|
|
// minor tweaks to clean a few additional things up after the rewrite.
|
|
//
|
|
// Revision 1.14 1998/09/19 17:59:45 curt
|
|
// Use c++ streams (fg_gzifstream). Also converted many character arrays to
|
|
// the string class.
|
|
//
|
|
// Revision 1.13 1998/09/09 16:24:04 curt
|
|
// Fixed a bug in the handling of exclude files which was causing
|
|
// a crash by calling fclose() on an invalid file handle.
|
|
//
|
|
// Revision 1.12 1998/08/24 20:03:31 curt
|
|
// Eliminated a possible memory overrun error.
|
|
// Use the proper free() rather than the incorrect delete().
|
|
//
|
|
// Revision 1.11 1998/07/20 12:46:11 curt
|
|
// When outputing to a .node file, first 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. This is my first pass at adding an area "cutout"
|
|
// feature to the terrain generation pipeline.
|
|
//
|
|
// Revision 1.10 1998/07/13 20:58:02 curt
|
|
// .
|
|
//
|
|
// Revision 1.9 1998/07/13 15:29:49 curt
|
|
// Added #ifdef HAVE_CONFIG_H
|
|
//
|
|
// Revision 1.8 1998/07/04 00:47:18 curt
|
|
// typedef'd struct fgBUCKET.
|
|
//
|
|
// Revision 1.7 1998/06/05 18:14:39 curt
|
|
// Abort out early when reading the "A" record if it doesn't look like
|
|
// a proper DEM file.
|
|
//
|
|
// Revision 1.6 1998/05/02 01:49:21 curt
|
|
// Fixed a bug where the wrong variable was being initialized.
|
|
//
|
|
// Revision 1.5 1998/04/25 15:00:32 curt
|
|
// Changed "r" to "rb" in gzopen() options. This fixes bad behavior in win32.
|
|
//
|
|
// Revision 1.4 1998/04/22 13:14:46 curt
|
|
// Fixed a bug in zlib usage.
|
|
//
|
|
// Revision 1.3 1998/04/18 03:53:05 curt
|
|
// Added zlib support.
|
|
//
|
|
// Revision 1.2 1998/04/14 02:43:27 curt
|
|
// Used "new" to auto-allocate large DEM parsing arrays in class constructor.
|
|
//
|
|
// Revision 1.1 1998/04/08 22:57:22 curt
|
|
// Adopted Gnu automake/autoconf system.
|
|
//
|
|
// Revision 1.3 1998/04/06 21:09:41 curt
|
|
// Additional win32 support.
|
|
// Fixed a bad bug in dem file parsing that was causing the output to be
|
|
// flipped about x = y.
|
|
//
|
|
// Revision 1.2 1998/03/23 20:35:41 curt
|
|
// Updated to use FG_EPSILON
|
|
//
|
|
// Revision 1.1 1998/03/19 02:54:47 curt
|
|
// Reorganized into a class lib called fgDEM.
|
|
//
|
|
// Revision 1.1 1998/03/19 01:46:28 curt
|
|
// Initial revision.
|
|
//
|