1 // array.cxx -- Array management class
3 // Written by Curtis Olson, started March 1998.
5 // Copyright (C) 1998 - 1999 Curtis L. Olson - curt@flightgear.org
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.
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.
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.
22 // (Log is kept at end of this file)
29 #include <Include/compiler.h>
31 // #include <ctype.h> // isspace()
32 // #include <stdlib.h> // atoi()
33 // #include <math.h> // rint()
35 // #include <string.h>
36 // #ifdef HAVE_SYS_STAT_H
37 // # include <sys/stat.h> // stat()
39 // #ifdef FG_HAVE_STD_INCLUDES
42 // # include <errno.h>
44 // #ifdef HAVE_UNISTD_H
45 // # include <unistd.h> // stat()
50 #include <Include/fg_constants.h>
51 #include <Misc/fgstream.hxx>
52 #include <Misc/strutils.hxx>
53 #include <Math/leastsqs.hxx>
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];
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];
78 FGArray::open( const string& file ) {
79 // open input file (or read from stdin)
81 cout << " Opening array data pipe from stdin" << endl;
83 // fd = gzdopen(STDIN_FILENO, "r");
84 cout << " Not yet ported ..." << endl;
87 in = new fg_gzifstream( file );
89 cout << " Cannot open " << file << endl;
92 cout << " Opening array data file: " << file << endl;
99 // close an Array file
102 // the fg_gzifstream doesn't seem to have a close()
113 *in >> originx >> originy;
114 *in >> cols >> col_step;
115 *in >> rows >> row_step;
117 cout << " origin = " << originx << " " << originy << endl;
118 cout << " cols = " << cols << " rows = " << rows << endl;
119 cout << " col_step = " << col_step << " row_step = " << row_step <<endl;
121 for ( int i = 0; i < cols; i++ ) {
122 for ( int j = 0; j < rows; j++ ) {
123 *in >> in_data[i][j];
127 cout << " Done parsing\n";
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) );
143 // Initialize output mesh structure
144 void FGArray::outputmesh_init( void ) {
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;
155 // Get the value of a mesh node
156 double FGArray::outputmesh_get_pt( int i, int j ) {
157 return ( out_data[i][j] );
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;
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;
176 int n, row, start, end;
177 int colmin, colmax, rowmin, rowmax;
179 // FILE *dem, *fit, *fit1;
181 error_sq = error * error;
183 cout << " Initializing output mesh structure" << endl;
184 // outputmesh_init();
186 // determine dimensions
191 cout << " Fitting region = " << colmin << "," << rowmin << " to "
192 << colmax << "," << rowmax << endl;;
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] );
200 cout << " Beginning best fit procedure" << endl;
203 for ( row = rowmin; row < rowmax; row++ ) {
204 // fit = fopen("fit.dat", "w");
205 // fit1 = fopen("fit1.dat", "w");
209 // cout << " fitting row = " << row << endl;
211 while ( start < colmax - 1 ) {
215 x[0] = start * col_step;
216 y[0] = in_data[start][row];
218 x[1] = end * col_step;
219 y[1] = in_data[end][row];
223 // cout << "Least square of first 2 points" << endl;
224 least_squares(x, y, n, &m, &b);
228 while ( (end < colmax) && good_fit ) {
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);
238 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
239 start, end, ave_error, max_error, m, b);
241 f = fopen("gnuplot.dat", "w");
242 for ( j = 0; j <= end; j++) {
243 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
246 for ( j = start; j <= end; j++) {
247 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
252 printf("Please hit return: "); gets(junk);
255 if ( max_error > error_sq ) {
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
267 // we popped out of the above loop while still within
268 // the error threshold, so we must be at the end of
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);
279 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
280 start, end, ave_error, max_error, m, b);
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);
287 if ( start > colmin ) {
288 // skip this for the first line segment
290 add_fit_node( start, row, (lasty + cury) / 2 );
291 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
294 lasty = m * x[end-start] + b;
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 ),
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);
317 // outputmesh_output_nodes(fg_root, p);
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
326 double xlocal, ylocal, dx, dy, zA, zB, elev;
327 int x1, x2, x3, y1, y2, y3;
331 /* determine if we are in the lower triangle or the upper triangle
339 then calculate our end points
342 xlocal = (lon - originx) / col_step;
343 ylocal = (lat - originy) / row_step;
345 xindex = (int)(xlocal);
346 yindex = (int)(ylocal);
348 // printf("xindex = %d yindex = %d\n", xindex, yindex);
350 if ( xindex + 1 == cols ) {
354 if ( yindex + 1 == rows ) {
358 if ( (xindex < 0) || (xindex + 1 >= cols) ||
359 (yindex < 0) || (yindex + 1 >= rows) ) {
360 cout << "WARNING: Attempt to interpolate value outside of array!!!"
365 dx = xlocal - xindex;
366 dy = ylocal - yindex;
370 // printf(" Lower triangle\n");
374 z1 = in_data[x1][y1];
378 z2 = in_data[x2][y2];
382 z3 = in_data[x3][y3];
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);
389 zA = dx * (z2 - z1) + z1;
390 zB = dx * (z3 - z1) + z1;
392 // printf(" zA = %.2f zB = %.2f\n", zA, zB);
394 if ( dx > FG_EPSILON ) {
395 elev = dy * (zB - zA) / dx + zA;
401 // printf(" Upper triangle\n");
405 z1 = in_data[x1][y1];
409 z2 = in_data[x2][y2];
413 z3 = in_data[x3][y3];
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);
420 zA = dy * (z2 - z1) + z1;
421 zB = dy * (z3 - z1) + z1;
423 // printf(" zA = %.2f zB = %.2f\n", zA, zB );
424 // printf(" xB - xA = %.2f\n", col_step * dy / row_step);
426 if ( dy > FG_EPSILON ) {
427 elev = dx * (zB - zA) / dy + zA;
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 )
444 double exnodes[MAX_EX_NODES][3];
445 struct stat stat_buf;
447 char file[256], exfile[256];
454 int colmin, colmax, rowmin, rowmax;
455 int i, j, count, excount, result;
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";
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;
471 // stat() directory and create if needed
473 result = stat(dir.c_str(), &stat_buf);
474 if ( result != 0 && errno == ENOENT ) {
475 cout << " Creating directory\n";
477 command = "mkdir -p " + dir + "\n";
478 system( command.c_str() );
480 // assume directory exists
483 // get index and generate output file name
484 index = p.gen_index();
485 sprintf(file, "%s/%ld.node", dir.c_str(), index);
487 // get (optional) extra node file name (in case there is matching
489 strcpy(exfile, file);
490 strcat(exfile, ".ex");
492 // load extra nodes if they exist
494 if ( (fd = fopen(exfile, "r")) != NULL ) {
496 fscanf(fd, "%d %d %d %d", &excount, &junki, &junki, &junki);
498 if ( excount > MAX_EX_NODES - 1 ) {
499 printf("Error, too many 'extra' nodes, increase array size\n");
502 printf(" Expecting %d 'extra' nodes\n", excount);
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]);
514 printf("Creating node file: %s\n", file);
515 fd = fopen(file, "w");
517 // first count regular nodes to generate header
519 for ( j = rowmin; j <= rowmax; j++ ) {
520 for ( i = colmin; i <= colmax; i++ ) {
521 if ( out_data[i][j] > -9000.0 ) {
525 // printf(" count = %d\n", count);
527 fprintf(fd, "%d 2 1 0\n", count + excount);
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]);
535 // write out actual node data
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",
542 originx + (double)i * col_step,
543 originy + (double)j * row_step,
547 // printf(" count = %d\n", count);
555 FGArray::~FGArray( void ) {
556 // printf("class FGArray DEstructor called.\n");
558 // delete [] out_data;
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.
567 // Revision 1.3 1999/03/17 23:48:17 curt
568 // Removed forced -g compile flag.
569 // Fixed a couple compiler warnings.
571 // Revision 1.2 1999/03/13 23:50:26 curt
572 // Tweaked output formatting a bit.
574 // Revision 1.1 1999/03/13 18:45:02 curt
575 // Initial revision. (derived from libDEM.a code.)