3 // dem.c -- DEM management class
5 // Written by Curtis Olson, started March 1998.
7 // Copyright (C) 1998 Curtis L. Olson - curt@me.umn.edu
9 // This program is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU General Public License as
11 // published by the Free Software Foundation; either version 2 of the
12 // License, or (at your option) any later version.
14 // This program is distributed in the hope that it will be useful, but
15 // WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // General Public License for more details.
19 // You should have received a copy of the GNU General Public License
20 // along with this program; if not, write to the Free Software
21 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 // (Log is kept at end of this file)
31 #include <ctype.h> // isspace()
32 #include <stdlib.h> // atoi()
33 #include <math.h> // rint()
36 #include <sys/stat.h> // stat()
37 #include <unistd.h> // stat()
39 #include <zlib/zlib.h>
42 #include "leastsqs.hxx"
44 #include <Include/fg_constants.h>
47 #define MAX_EX_NODES 10000
50 # define MKDIR(a) mkdir(a,S_IRWXU) // I am just guessing at this flag (NHV)
54 fgDEM::fgDEM( void ) {
55 // printf("class fgDEM CONstructor called.\n");
56 dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
57 output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
63 // return the file path name ( foo/bar/file.ext = foo/bar )
64 void extract_path (char *in, char *base) {
71 while ( (i >= 0) && (in[i] != '/') ) {
79 // Make a subdirectory
80 int my_mkdir (char *dir) {
86 result = stat (dir, &stat_buf);
90 result = stat (dir, &stat_buf);
92 printf ("problem creating %s\n", dir);
94 printf ("%s created\n", dir);
97 printf ("%s already exists\n", dir);
107 int fgDEM::open ( char *file ) {
108 // open input file (or read from stdin)
109 if ( strcmp(file, "-") == 0 ) {
110 printf("Loading DEM data file: stdin\n");
112 fd = gzdopen(STDIN_FILENO, "r");
114 if ( (fd = gzopen(file, "rb")) == NULL ) {
115 printf("Cannot gzopen %s\n", file);
118 printf("Loading DEM data file: %s\n", file);
126 int fgDEM::close ( void ) {
133 // return next token from input stream
134 static void next_token(gzFile fd, char *token) {
141 while ( (c != -1) && (c == ' ') ) {
144 while ( (c != -1) && (c != ' ') && (c != '\n') ){
152 strcpy(token, "__END_OF_FILE__");
153 printf(" Warning: Reached end of file!\n");
156 // printf(" returning %s\n", token);
160 // return next integer from input stream
161 static int next_int(gzFile fd) {
164 next_token(fd, token);
165 return ( atoi(token) );
169 // return next double from input stream
170 static double next_double(gzFile fd) {
173 next_token(fd, token);
174 return ( atof(token) );
178 // return next exponential num from input stream
179 static int next_exp(gzFile fd) {
185 next_token(fd, token);
187 sscanf(token, "%lfD%d", &mantissa, &exp);
189 // printf(" Mantissa = %.4f Exp = %d\n", mantissa, exp);
193 for ( i = 1; i <= exp; i++ ) {
196 } else if ( exp < 0 ) {
197 for ( i = -1; i >= exp; i-- ) {
202 return( (int)rint(mantissa * (double)acc) );
206 // read and parse DEM "A" record
207 int fgDEM::read_a_record( void ) {
214 // get the name field (144 characters)
215 for ( i = 0; i < 144; i++ ) {
216 name[i] = gzgetc(fd);
220 // clean off the whitespace at the end
221 for ( i = strlen(name)-2; i > 0; i-- ) {
222 if ( !isspace(name[i]) ) {
228 printf(" Quad name field: %s\n", name);
230 // DEM level code, 3 reflects processing by DMA
232 printf(" DEM level code = %d\n", inum);
238 // Pattern code, 1 indicates a regular elevation pattern
240 printf(" Pattern code = %d\n", inum);
242 // Planimetric reference system code, 0 indicates geographic
243 // coordinate system.
245 printf(" Planimetric reference code = %d\n", inum);
249 printf(" Zone code = %d\n", inum);
251 // Map projection parameters (ignored)
252 for ( i = 0; i < 15; i++ ) {
253 dnum = next_double(fd);
254 // printf("%d: %f\n",i,dnum);
257 // Units code, 3 represents arc-seconds as the unit of measure for
258 // ground planimetric coordinates throughout the file.
261 printf(" Unknown (X,Y) units code = %d!\n", inum);
265 // Units code; 2 represents meters as the unit of measure for
266 // elevation coordinates throughout the file.
269 printf(" Unknown (Z) units code = %d!\n", inum);
273 // Number (n) of sides in the polygon which defines the coverage of
274 // the DEM file (usually equal to 4).
277 printf(" Unknown polygon dimension = %d!\n", inum);
281 // Ground coordinates of bounding box in arc-seconds
282 dem_x1 = originx = next_exp(fd);
283 dem_y1 = originy = next_exp(fd);
284 printf(" Origin = (%.2f,%.2f)\n", originx, originy);
286 dem_x2 = next_exp(fd);
287 dem_y2 = next_exp(fd);
289 dem_x3 = next_exp(fd);
290 dem_y3 = next_exp(fd);
292 dem_x4 = next_exp(fd);
293 dem_y4 = next_exp(fd);
295 // Minimum/maximum elevations in meters
296 dem_z1 = next_exp(fd);
297 dem_z2 = next_exp(fd);
298 printf(" Elevation range %.4f %.4f\n", dem_z1, dem_z2);
300 // Counterclockwise angle from the primary axis of ground
301 // planimetric referenced to the primary axis of the DEM local
303 next_token(fd, token);
305 // Accuracy code; 0 indicates that a record of accuracy does not
306 // exist and that no record type C will follow.
308 // DEM spacial resolution. Usually (3,3,1) (3,6,1) or (3,9,1)
309 // depending on latitude
311 // I will eventually have to do something with this for data at
312 // higher latitudes */
313 next_token(fd, token);
314 printf(" accuracy & spacial resolution string = %s\n", token);
316 printf(" length = %d\n", i);
318 ptr = token + i - 12;
319 printf(" last field = %s = %.2f\n", ptr, atof(ptr));
323 col_step = atof(ptr);
324 printf(" last field = %s = %.2f\n", ptr, col_step);
328 row_step = atof(ptr);
329 printf(" last field = %s = %.2f\n", ptr, row_step);
332 // accuracy code = atod(token)
334 printf(" Accuracy code = %d\n", inum);
336 printf(" column step = %.2f row step = %.2f\n",
338 // dimension of arrays to follow (1)
339 next_token(fd, token);
341 // number of profiles
342 dem_num_profiles = cols = next_int(fd);
343 printf(" Expecting %d profiles\n", dem_num_profiles);
349 // read and parse DEM "B" record
350 void fgDEM::read_b_record( void ) {
354 // row / column id of this profile
355 prof_row = next_int(fd);
356 prof_col = next_int(fd);
357 // printf("col id = %d row id = %d\n", prof_col, prof_row);
359 // Number of columns and rows (elevations) in this profile
360 prof_num_rows = rows = next_int(fd);
361 prof_num_cols = next_int(fd);
362 // printf(" profile num rows = %d\n", prof_num_rows);
364 // Ground planimetric coordinates (arc-seconds) of the first
365 // elevation in the profile
366 prof_x1 = next_exp(fd);
367 prof_y1 = next_exp(fd);
368 // printf(" Starting at %.2f %.2f\n", prof_x1, prof_y1);
370 // Elevation of local datum for the profile. Always zero for
371 // 1-degree DEM, the reference is mean sea level.
372 next_token(fd, token);
374 // Minimum and maximum elevations for the profile.
375 next_token(fd, token);
376 next_token(fd, token);
378 // One (usually) dimensional array (prof_num_cols,1) of elevations
379 for ( i = 0; i < prof_num_rows; i++ ) {
380 prof_data = next_int(fd);
381 dem_data[cur_col][i] = (float)prof_data;
387 int fgDEM::parse( void ) {
392 if ( !read_a_record() ) {
396 for ( i = 0; i < dem_num_profiles; i++ ) {
397 // printf("Ready to read next b record\n");
401 if ( cur_col % 100 == 0 ) {
402 printf(" loaded %d profiles of data\n", cur_col);
406 printf(" Done parsing\n");
412 // return the current altitude based on mesh data. We should rewrite
413 // this to interpolate exact values, but for now this is good enough
414 double fgDEM::interpolate_altitude( double lon, double lat ) {
415 // we expect incoming (lon,lat) to be in arcsec for now
417 double xlocal, ylocal, dx, dy, zA, zB, elev;
418 int x1, x2, x3, y1, y2, y3;
422 /* determine if we are in the lower triangle or the upper triangle
430 then calculate our end points
433 xlocal = (lon - originx) / col_step;
434 ylocal = (lat - originy) / row_step;
436 xindex = (int)(xlocal);
437 yindex = (int)(ylocal);
439 // printf("xindex = %d yindex = %d\n", xindex, yindex);
441 if ( xindex + 1 == cols ) {
445 if ( yindex + 1 == rows ) {
449 if ( (xindex < 0) || (xindex + 1 >= cols) ||
450 (yindex < 0) || (yindex + 1 >= rows) ) {
454 dx = xlocal - xindex;
455 dy = ylocal - yindex;
459 // printf(" Lower triangle\n");
463 z1 = dem_data[x1][y1];
467 z2 = dem_data[x2][y2];
471 z3 = dem_data[x3][y3];
473 // printf(" dx = %.2f dy = %.2f\n", dx, dy);
474 // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
475 // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
476 // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
478 zA = dx * (z2 - z1) + z1;
479 zB = dx * (z3 - z1) + z1;
481 // printf(" zA = %.2f zB = %.2f\n", zA, zB);
483 if ( dx > FG_EPSILON ) {
484 elev = dy * (zB - zA) / dx + zA;
490 // printf(" Upper triangle\n");
494 z1 = dem_data[x1][y1];
498 z2 = dem_data[x2][y2];
502 z3 = dem_data[x3][y3];
504 // printf(" dx = %.2f dy = %.2f\n", dx, dy);
505 // printf(" (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
506 // printf(" (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
507 // printf(" (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
509 zA = dy * (z2 - z1) + z1;
510 zB = dy * (z3 - z1) + z1;
512 // printf(" zA = %.2f zB = %.2f\n", zA, zB );
513 // printf(" xB - xA = %.2f\n", col_step * dy / row_step);
515 if ( dy > FG_EPSILON ) {
516 elev = dx * (zB - zA) / dy + zA;
526 // Use least squares to fit a simpler data set to dem data
527 void fgDEM::fit( double error, fgBUCKET *p ) {
528 double x[DEM_SIZE_1], y[DEM_SIZE_1];
529 double m, b, ave_error, max_error;
531 int n, row, start, end, good_fit;
532 int colmin, colmax, rowmin, rowmax;
533 // FILE *dem, *fit, *fit1;
535 printf("Initializing output mesh structure\n");
538 // determine dimensions
539 colmin = p->x * ( (cols - 1) / 8);
540 colmax = colmin + ( (cols - 1) / 8);
541 rowmin = p->y * ( (rows - 1) / 8);
542 rowmax = rowmin + ( (rows - 1) / 8);
543 printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
545 // include the corners explicitly
546 outputmesh_set_pt(colmin, rowmin, dem_data[colmin][rowmin]);
547 outputmesh_set_pt(colmin, rowmax, dem_data[colmin][rowmax]);
548 outputmesh_set_pt(colmax, rowmax, dem_data[colmax][rowmax]);
549 outputmesh_set_pt(colmax, rowmin, dem_data[colmax][rowmin]);
551 printf("Beginning best fit procedure\n");
553 for ( row = rowmin; row <= rowmax; row++ ) {
554 // fit = fopen("fit.dat", "w");
555 // fit1 = fopen("fit1.dat", "w");
559 // printf(" fitting row = %d\n", row);
561 while ( start < colmax ) {
565 x[(end - start) - 1] = 0.0 + ( start * col_step );
566 y[(end - start) - 1] = dem_data[start][row];
568 while ( (end <= colmax) && good_fit ) {
569 n = (end - start) + 1;
570 // printf("Least square of first %d points\n", n);
571 x[end - start] = 0.0 + ( end * col_step );
572 y[end - start] = dem_data[end][row];
573 least_squares(x, y, n, &m, &b);
574 ave_error = least_squares_error(x, y, n, m, b);
575 max_error = least_squares_max_error(x, y, n, m, b);
578 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
579 start, end, ave_error, max_error, m, b);
581 f = fopen("gnuplot.dat", "w");
582 for ( j = 0; j <= end; j++) {
583 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
586 for ( j = start; j <= end; j++) {
587 fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ),
592 printf("Please hit return: "); gets(junk);
595 if ( max_error > error ) {
603 // error exceeded the threshold, back up
604 end -= 2; // back "end" up to the last good enough fit
605 n--; // back "n" up appropriately too
607 // we popped out of the above loop while still within
608 // the error threshold, so we must be at the end of
613 least_squares(x, y, n, &m, &b);
614 ave_error = least_squares_error(x, y, n, m, b);
615 max_error = least_squares_max_error(x, y, n, m, b);
619 printf("%d - %d ave error = %.2f max error = %.2f y = %.2f*x + %.2f\n",
620 start, end, ave_error, max_error, m, b);
623 fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
624 fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
627 if ( start > colmin ) {
628 // skip this for the first line segment
630 outputmesh_set_pt(start, row, (lasty + cury) / 2);
631 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
634 lasty = m * x[end-start] + b;
642 dem = fopen("gnuplot.dat", "w");
643 for ( j = 0; j < DEM_SIZE_1; j++) {
644 fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ),
650 // NOTICE, this is for testing only. This instance of
651 // output_nodes should be removed. It should be called only
652 // once at the end once all the nodes have been generated.
653 // newmesh_output_nodes(&nm, "mesh.node");
654 // printf("Please hit return: "); gets(junk);
657 // outputmesh_output_nodes(fg_root, p);
661 // Initialize output mesh structure
662 void fgDEM::outputmesh_init( void ) {
665 for ( j = 0; j < DEM_SIZE_1; j++ ) {
666 for ( i = 0; i < DEM_SIZE_1; i++ ) {
667 output_data[i][j] = -9999.0;
673 // Get the value of a mesh node
674 double fgDEM::outputmesh_get_pt( int i, int j ) {
675 return ( output_data[i][j] );
679 // Set the value of a mesh node
680 void fgDEM::outputmesh_set_pt( int i, int j, double value ) {
681 // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
682 output_data[i][j] = value;
686 // Write out a node file that can be used by the "triangle" program.
687 // Check for an optional "index.node.ex" file in case there is a .poly
688 // file to go along with this node file. Include these nodes first
689 // since they are referenced by position from the .poly file.
690 void fgDEM::outputmesh_output_nodes( char *fg_root, fgBUCKET *p ) {
691 double exnodes[MAX_EX_NODES][3];
692 double junk1, junk2, junk3;
693 struct stat stat_buf;
694 char base_path[256], dir[256], file[256], exfile[256];
701 int colmin, colmax, rowmin, rowmax;
702 int i, j, count, excount, result;
704 // determine dimensions
705 colmin = p->x * ( (cols - 1) / 8);
706 colmax = colmin + ( (cols - 1) / 8);
707 rowmin = p->y * ( (rows - 1) / 8);
708 rowmax = rowmin + ( (rows - 1) / 8);
709 printf(" dumping region = %d,%d to %d,%d\n",
710 colmin, rowmin, colmax, rowmax);
712 // generate the base directory
713 fgBucketGenBasePath(p, base_path);
714 printf("fg_root = %s Base Path = %s\n", fg_root, base_path);
715 sprintf(dir, "%s/Scenery/%s", fg_root, base_path);
716 printf("Dir = %s\n", dir);
718 // stat() directory and create if needed
719 result = stat(dir, &stat_buf);
721 printf("Stat error need to create directory\n");
725 sprintf(command, "mkdir -p %s\n", dir);
730 // Cygwin crashes when trying to output to node file
731 // explicitly making directory structure seems OK on Win95
733 extract_path (base_path, tmp_path);
735 sprintf (dir, "%s/Scenery", fg_root);
736 if (my_mkdir (dir)) { exit (-1); }
738 sprintf (dir, "%s/Scenery/%s", fg_root, tmp_path);
739 if (my_mkdir (dir)) { exit (-1); }
741 sprintf (dir, "%s/Scenery/%s", fg_root, base_path);
742 if (my_mkdir (dir)) { exit (-1); }
747 // assume directory exists
750 // get index and generate output file name
751 index = fgBucketGenIndex(p);
752 sprintf(file, "%s/%ld.node", dir, index);
754 // get (optional) extra node file name (in case there is matching
756 strcpy(exfile, file);
757 strcat(exfile, ".ex");
759 // load extra nodes if they exist
761 if ( (fd = fopen(exfile, "r")) != NULL ) {
762 fscanf(fd, "%d %d %d %d", &excount, &junk1, &junk2, &junk3);
764 if ( excount > MAX_EX_NODES - 1 ) {
765 printf("Error, too many 'extra' nodes, increase array size\n");
768 printf(" Expecting %d 'extra' nodes\n", excount);
771 for ( i = 1; i <= excount; i++ ) {
772 fscanf(fd, "%d %lf %lf %lf\n", &junk1,
773 &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]);
774 printf("(extra) %d %.2f %.2f %.2f\n",
775 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
780 printf("Creating node file: %s\n", file);
781 fd = fopen(file, "w");
783 // first count regular nodes to generate header
785 for ( j = rowmin; j <= rowmax; j++ ) {
786 for ( i = colmin; i <= colmax; i++ ) {
787 if ( output_data[i][j] > -9000.0 ) {
791 // printf(" count = %d\n", count);
793 fprintf(fd, "%d 2 1 0\n", count + excount);
795 // now write out extra node data
796 for ( i = 1; i <= excount; i++ ) {
797 fprintf(fd, "%d %.2f %.2f %.2f\n",
798 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
801 // write out actual node data
803 for ( j = rowmin; j <= rowmax; j++ ) {
804 for ( i = colmin; i <= colmax; i++ ) {
805 if ( output_data[i][j] > -9000.0 ) {
806 fprintf(fd, "%d %.2f %.2f %.2f\n",
808 originx + (double)i * col_step,
809 originy + (double)j * row_step,
813 // printf(" count = %d\n", count);
820 fgDEM::~fgDEM( void ) {
821 // printf("class fgDEM DEstructor called.\n");
828 // Revision 1.12 1998/08/24 20:03:31 curt
829 // Eliminated a possible memory overrun error.
830 // Use the proper free() rather than the incorrect delete().
832 // Revision 1.11 1998/07/20 12:46:11 curt
833 // When outputing to a .node file, first check for an optional
834 // "index.node.ex" file in case there is a .poly file to go along with this
835 // node file. Include these nodes first since they are referenced by position
836 // from the .poly file. This is my first pass at adding an area "cutout"
837 // feature to the terrain generation pipeline.
839 // Revision 1.10 1998/07/13 20:58:02 curt
842 // Revision 1.9 1998/07/13 15:29:49 curt
843 // Added #ifdef HAVE_CONFIG_H
845 // Revision 1.8 1998/07/04 00:47:18 curt
846 // typedef'd struct fgBUCKET.
848 // Revision 1.7 1998/06/05 18:14:39 curt
849 // Abort out early when reading the "A" record if it doesn't look like
850 // a proper DEM file.
852 // Revision 1.6 1998/05/02 01:49:21 curt
853 // Fixed a bug where the wrong variable was being initialized.
855 // Revision 1.5 1998/04/25 15:00:32 curt
856 // Changed "r" to "rb" in gzopen() options. This fixes bad behavior in win32.
858 // Revision 1.4 1998/04/22 13:14:46 curt
859 // Fixed a bug in zlib usage.
861 // Revision 1.3 1998/04/18 03:53:05 curt
862 // Added zlib support.
864 // Revision 1.2 1998/04/14 02:43:27 curt
865 // Used "new" to auto-allocate large DEM parsing arrays in class constructor.
867 // Revision 1.1 1998/04/08 22:57:22 curt
868 // Adopted Gnu automake/autoconf system.
870 // Revision 1.3 1998/04/06 21:09:41 curt
871 // Additional win32 support.
872 // Fixed a bad bug in dem file parsing that was causing the output to be
873 // flipped about x = y.
875 // Revision 1.2 1998/03/23 20:35:41 curt
876 // Updated to use FG_EPSILON
878 // Revision 1.1 1998/03/19 02:54:47 curt
879 // Reorganized into a class lib called fgDEM.
881 // Revision 1.1 1998/03/19 01:46:28 curt