fly/flightgear
1
0
Fork 0
Code Activity

Initial revision. (derived from libDEM.a code.)

Browse source
This commit is contained in:
curt 1999-03-13 18:45:02 +00:00
parent 5fe3a80ae3
commit 2e20f32267
4 changed files with 752 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

22
Array/Makefile.am Normal file
View file

@ -0,0 +1,22 @@
noinst_LIBRARIES = libArray.a
libArray_a_SOURCES = array.cxx array.hxx
bin_PROGRAMS = testarray
testarray_SOURCES = testarray.cxx
testarray_LDADD = \
$(top_builddir)/Tools/Construct/Array/libArray.a \
$(top_builddir)/Lib/Bucket/libBucket.a \
$(top_builddir)/Lib/Math/libMath.a \
$(top_builddir)/Lib/Misc/libMisc.a \
$(top_builddir)/Lib/zlib/libz.a
INCLUDES += -I$(top_builddir) -I$(top_builddir)/Lib
# We can't build this with "-O2" (optimization) since this causes a seg fault
# I haven't found a way to strip this out of the CXXFLAGS, so I'm just
# setting it to "-g"
CXXFLAGS = -g

557
Array/array.cxx Normal file
View file

@ -0,0 +1,557 @@
// array.cxx -- Array management class
//
// Written by Curtis Olson, started March 1998.
//
// Copyright (C) 1998 - 1999 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
// 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 STL_STRING
#include <Include/fg_constants.h>
#include <Misc/fgstream.hxx>
#include <Misc/strutils.hxx>
#include <Math/leastsqs.hxx>
#include "array.hxx"
FG_USING_STD(string);
FGArray::FGArray( void ) {
// cout << "class FGArray CONstructor called." << endl;
in_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
out_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
}
FGArray::FGArray( const string &file ) {
// cout << "class FGArray CONstructor called." << endl;
in_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
out_data = new float[ARRAY_SIZE_1][ARRAY_SIZE_1];
FGArray::open(file);
}
// open an Array file
int
FGArray::open( const string& file ) {
// open input file (or read from stdin)
if ( file == "-" ) {
cout << "Opening array data pipe from stdin" << endl;
// fd = stdin;
// fd = gzdopen(STDIN_FILENO, "r");
cout << "Not yet ported ..." << endl;
return 0;
} else {
in = new fg_gzifstream( file );
if ( !(*in) ) {
cout << "Cannot open " << file << endl;
return 0;
}
cout << "Opening array data file: " << file << endl;
}
return 1;
}
// close an Array file
int
FGArray::close() {
// the fg_gzifstream doesn't seem to have a close()
delete in;
return 1;
}
// parse Array file
int
FGArray::parse() {
int i;
// cur_col = 0;
*in >> originx >> originy;
*in >> cols >> col_step;
*in >> rows >> row_step;
cout << " origin = " << originx << " " << originy << endl;
cout << " cols = " << cols << " rows = " << rows << endl;
cout << " col_step = " << col_step << " row_step = " << row_step <<endl;
for ( int i = 0; i < cols; i++ ) {
for ( int j = 0; j < rows; j++ ) {
*in >> in_data[i][j];
}
}
cout << " Done parsing\n";
return 1;
}
// Initialize output mesh structure
void FGArray::outputmesh_init( void ) {
int i, j;
for ( j = 0; j < ARRAY_SIZE_1; j++ ) {
for ( i = 0; i < ARRAY_SIZE_1; i++ ) {
out_data[i][j] = -9999.0;
}
}
}
// Get the value of a mesh node
double FGArray::outputmesh_get_pt( int i, int j ) {
return ( out_data[i][j] );
}
// Set the value of a mesh node
void FGArray::outputmesh_set_pt( int i, int j, double value ) {
// cout << "Setting data[" << i << "][" << j << "] = " << value << endl;
out_data[i][j] = value;
}
// Use least squares to fit a simpler data set to dem data
void FGArray::fit( FGBucket& p, double error ) {
double x[ARRAY_SIZE_1], y[ARRAY_SIZE_1];
double m, b, max_error, error_sq;
double x1, y1;
// double ave_error;
double cury, lasty;
int n, row, start, end;
int colmin, colmax, rowmin, rowmax;
bool good_fit;
// FILE *dem, *fit, *fit1;
error_sq = error * error;
cout << "Initializing output mesh structure" << endl;
outputmesh_init();
// determine dimensions
colmin = 0;
colmax = cols;
rowmin = 0;
rowmax = rows;
cout << "Fitting region = " << colmin << "," << rowmin << " to "
<< colmax << "," << rowmax << endl;;
// include the corners explicitly
outputmesh_set_pt(colmin, rowmin, in_data[colmin][rowmin]);
outputmesh_set_pt(colmin, rowmax, in_data[colmin][rowmax]);
outputmesh_set_pt(colmax, rowmax, in_data[colmax][rowmax]);
outputmesh_set_pt(colmax, rowmin, in_data[colmax][rowmin]);
cout << "Beginning best fit procedure" << endl;
for ( row = rowmin; row < rowmax; row++ ) {
// fit = fopen("fit.dat", "w");
// fit1 = fopen("fit1.dat", "w");
start = colmin;
// cout << " fitting row = " << row << endl;
while ( start < colmax - 1 ) {
end = start + 1;
good_fit = true;
x[0] = start * col_step;
y[0] = in_data[start][row];
x[1] = end * col_step;
y[1] = in_data[end][row];
n = 2;
// cout << "Least square of first 2 points" << endl;
least_squares(x, y, n, &m, &b);
end++;
while ( (end < colmax) && good_fit ) {
++n;
// cout << "Least square of first " << n << " points" << endl;
x[n-1] = x1 = end * col_step;
y[n-1] = y1 = in_data[end][row];
least_squares_update(x1, y1, &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 ),
in_data[row][j]);
}
for ( j = start; j <= end; j++) {
fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
in_data[row][j]);
}
fclose(f);
printf("Please hit return: "); gets(junk);
*/
if ( max_error > error_sq ) {
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 < ARRAY_SIZE_1; j++) {
fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
in_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);
}
// return the current altitude based on grid data. We should rewrite
// this to interpolate exact values, but for now this is good enough
double FGArray::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) ) {
cout << "WARNING: Attempt to interpolate value outside of array!!!"
<< endl;
return(-9999);
}
dx = xlocal - xindex;
dy = ylocal - yindex;
if ( dx > dy ) {
// lower triangle
// printf(" Lower triangle\n");
x1 = xindex;
y1 = yindex;
z1 = in_data[x1][y1];
x2 = xindex + 1;
y2 = yindex;
z2 = in_data[x2][y2];
x3 = xindex + 1;
y3 = yindex + 1;
z3 = in_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 = in_data[x1][y1];
x2 = xindex;
y2 = yindex + 1;
z2 = in_data[x2][y2];
x3 = xindex + 1;
y3 = yindex + 1;
z3 = in_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);
}
#if 0
// 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 FGArray::outputmesh_output_nodes( const string& fg_root, FGBucket& 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 + "/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";
command = "mkdir -p " + dir + "\n";
system( command.c_str() );
} 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 ( out_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 ( out_data[i][j] > -9000.0 ) {
fprintf(fd, "%d %.2f %.2f %.2f\n",
count++,
originx + (double)i * col_step,
originy + (double)j * row_step,
out_data[i][j]);
}
}
// printf(" count = %d\n", count);
}
fclose(fd);
}
#endif
FGArray::~FGArray( void ) {
// printf("class FGArray DEstructor called.\n");
delete [] in_data;
delete [] out_data;
}
// $Log$
// Revision 1.1 1999/03/13 18:45:02 curt
// Initial revision. (derived from libDEM.a code.)
//

136
Array/array.hxx Normal file
View file

@ -0,0 +1,136 @@
// array.hxx -- Array management class
//
// Written by Curtis Olson, started March 1998.
//
// Copyright (C) 1998 - 1999 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
// 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)
#ifndef _ARRAY_HXX
#define _ARRAY_HXX
#ifndef __cplusplus
# error This library requires C++
#endif
#include <Bucket/newbucket.hxx>
#include <Misc/fgstream.hxx>
#define ARRAY_SIZE 1200
#define ARRAY_SIZE_1 1201
class FGArray {
private:
// file pointer for input
// gzFile fd;
fg_gzifstream *in;
// coordinates (in arc seconds) of south west corner
double originx, originy;
// number of columns and rows
int cols, rows;
// Distance between column and row data points (in arc seconds)
double col_step, row_step;
// pointers to the actual grid data allocated here
float (*in_data)[ARRAY_SIZE_1];
float (*out_data)[ARRAY_SIZE_1];
// Current "A" Record Information
// char dem_description[80], dem_quadrangle[80];
// double dem_x1, dem_y1, dem_x2, dem_y2, dem_x3, dem_y3, dem_x4, dem_y4;
// double dem_z1, dem_z2;
// int dem_resolution, dem_num_profiles;
// Current "B" Record Information
// int prof_col, prof_row;
// int prof_num_cols, prof_num_rows;
// double prof_x1, prof_y1;
// int prof_data;
// temporary values for the class to use
// char option_name[32];
// int do_data;
// int cur_col, cur_row;
// 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 );
#if 0
// Write out a node file that can be used by the "triangle" program
void outputmesh_output_nodes( const string& fg_root, FGBucket& p );
#endif
public:
// Constructor
FGArray( void );
FGArray( const string& file );
// Destructor
~FGArray( void );
// open an Array file (use "-" if input is coming from stdin)
int open ( const string& file );
// close a Array file
int close();
// parse a Array file
int parse();
// Use least squares to fit a simpler data set to dem data
void fit( FGBucket& p, double error );
// 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 );
// Informational methods
inline double get_originx() const { return originx; }
inline double get_originy() const { return originy; }
inline int get_cols() const { return cols; }
inline int get_rows() const { return rows; }
inline double get_col_step() const { return col_step; }
inline double get_row_step() const { return row_step; }
};
#endif // _ARRAY_HXX
// $Log$
// Revision 1.1 1999/03/13 18:45:02 curt
// Initial revision. (derived from libDEM.a code.)
//

37
Array/testarray.cxx Normal file
View file

@ -0,0 +1,37 @@
#include <Bucket/newbucket.hxx>
#include "array.hxx"
main(int argc, char **argv) {
char tile_name[256];
double lon, lat;
if ( argc != 2 ) {
cout << "Usage: " << argv[0] << " work_dir" << endl;
exit(-1);
}
string work_dir = argv[1];
lon = -146.248360; lat = 61.133950; // PAVD (Valdez, AK)
lon = -110.664244; lat = 33.352890; // P13
FGBucket b( lon, lat );
string base = b.gen_base_path();
string path = work_dir + "/Scenery/" + base;
long int b_index = b.gen_index();
sprintf(tile_name, "%ld", b_index);
string arrayfile = path + "/" + tile_name + ".dem";
cout << "arrayfile = " << arrayfile << endl;
FGArray a(arrayfile);
a.parse();
lon *= 3600;
lat *= 3600;
cout << a.interpolate_altitude(lon, lat) << endl;
a.fit(b, 100);
}