]> git.mxchange.org Git - flightgear.git/blob - DEM/dem.cxx
Use c++ streams (fg_gzifstream). Also converted many character arrays to
[flightgear.git] / DEM / dem.cxx
1 // -*- Mode: C++ -*-
2 //
3 // dem.c -- DEM management class
4 //
5 // Written by Curtis Olson, started March 1998.
6 //
7 // Copyright (C) 1998  Curtis L. Olson  - curt@me.umn.edu
8 //
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.
13 //
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.
18 //
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.
22 //
23 // $Id$
24 // (Log is kept at end of this file)
25
26
27 #ifdef HAVE_CONFIG_H
28 #  include <config.h>
29 #endif
30
31 #include <ctype.h>    // isspace()
32 #include <stdlib.h>   // atoi()
33 #include <math.h>     // rint()
34 #include <stdio.h>
35 #include <string.h>
36 #include <sys/stat.h> // stat()
37 #include <unistd.h>   // stat()
38
39 #include <string>
40
41 // #include <zlib/zlib.h>
42 #include <Misc/fgstream.hxx>
43 #include <Misc/strutils.hxx>
44
45 #include "dem.hxx"
46 #include "leastsqs.hxx"
47
48 #include <Include/fg_constants.h>
49
50
51 #define MAX_EX_NODES 10000
52
53 #ifdef WIN32
54 #  define MKDIR(a) mkdir(a,S_IRWXU)     // I am just guessing at this flag (NHV)
55 #endif // WIN32
56
57
58 fgDEM::fgDEM( void ) {
59     // printf("class fgDEM CONstructor called.\n");
60     dem_data = new float[DEM_SIZE_1][DEM_SIZE_1];
61     output_data = new float[DEM_SIZE_1][DEM_SIZE_1];
62 }
63
64
65 #ifdef WIN32
66
67 // return the file path name ( foo/bar/file.ext = foo/bar )
68 static void extract_path (char *in, char *base) {
69     int len, i;
70     
71     len = strlen (in);
72     strcpy (base, in);
73
74     i = len - 1;
75     while ( (i >= 0) && (in[i] != '/') ) {
76         i--;
77     }
78
79     base[i] = '\0';
80 }
81
82
83 // Make a subdirectory
84 static int my_mkdir (char *dir) {
85     struct stat stat_buf;
86     int result;
87
88     printf ("mk_dir() ");
89
90     result = stat (dir, &stat_buf);
91
92     if (result != 0) {
93         MKDIR (dir);
94         result = stat (dir, &stat_buf);
95         if (result != 0) {
96             printf ("problem creating %s\n", dir);
97         } else {
98             printf ("%s created\n", dir);
99         }
100     } else {
101         printf ("%s already exists\n", dir);
102     }
103
104     return (result);
105 }
106
107 #endif // WIN32
108
109
110 // open a DEM file
111 int fgDEM::open ( const string& file ) {
112     // open input file (or read from stdin)
113     if ( file ==  "-" ) {
114         printf("Loading DEM data file: stdin\n");
115         // fd = stdin;
116         // fd = gzdopen(STDIN_FILENO, "r");
117         printf("Not yet ported ...\n");
118         return 0;
119     } else {
120         in = new fg_gzifstream( file );
121         if ( !in ) {
122             cout << "Cannot open " + file + "\n";
123             return 0;
124         }
125         cout << "Loading DEM data file: " + file + "\n";
126     }
127
128     return 1;
129 }
130
131
132 // close a DEM file
133 int fgDEM::close () {
134     // the fg_gzifstream doesn't seem to have a close()
135
136     delete(in);
137
138     return 1;
139 }
140
141
142 // return next token from input stream
143 string fgDEM::next_token() {
144     string token;
145
146     in->stream() >> token;
147
148     cout << "    returning " + token + "\n";
149
150     return token;
151 }
152
153
154 // return next integer from input stream
155 int fgDEM::next_int() {
156     int result;
157
158     in->stream() >> result;
159
160     return result;
161 }
162
163
164 // return next double from input stream
165 double fgDEM::next_double() {
166     double result;
167
168     in->stream() >> result;
169
170     return result;
171 }
172
173
174 // return next exponential num from input stream
175 int fgDEM::next_exp() {
176     string token;
177     double mantissa;
178     int exp, acc;
179     int i;
180
181     token = next_token();
182
183     sscanf(token.c_str(), "%lfD%d", &mantissa, &exp);
184
185     cout << "    Mantissa = " << mantissa << "  Exp = " << exp << "\n";
186
187     acc = 1;
188     if ( exp > 0 ) {
189         for ( i = 1; i <= exp; i++ ) {
190             acc *= 10;
191         }
192     } else if ( exp < 0 ) {
193         for ( i = -1; i >= exp; i-- ) {
194             acc /= 10;
195         }
196     }
197
198     return( (int)rint(mantissa * (double)acc) );
199 }
200
201
202 // read and parse DEM "A" record
203 int fgDEM::read_a_record() {
204     int i, inum;
205     double dnum;
206     string name, token;
207     char c;
208     char *ptr;
209
210     // get the name field (144 characters)
211     for ( i = 0; i < 144; i++ ) {
212         in->get(c);
213         name += c;
214     }
215   
216     // clean off the trailing whitespace
217     name = trimleft(name);
218     cout << "    Quad name field: " + name + "\n";
219
220     // DEM level code, 3 reflects processing by DMA
221     inum = next_int();
222     cout << "    DEM level code = " << inum << "\n";
223
224     if ( inum > 3 ) {
225         return 0;
226     }
227
228     // Pattern code, 1 indicates a regular elevation pattern
229     inum = next_int();
230     cout << "    Pattern code = " << inum << "\n";
231
232     // Planimetric reference system code, 0 indicates geographic
233     // coordinate system.
234     inum = next_int();
235     cout << "    Planimetric reference code = " << inum << "\n";
236
237     // Zone code
238     inum = next_int();
239     cout << "    Zone code = " << inum << "\n";
240
241     // Map projection parameters (ignored)
242     for ( i = 0; i < 15; i++ ) {
243         dnum = next_double();
244         // printf("%d: %f\n",i,dnum);
245     }
246
247     // Units code, 3 represents arc-seconds as the unit of measure for
248     // ground planimetric coordinates throughout the file.
249     inum = next_int();
250     if ( inum != 3 ) {
251         cout << "    Unknown (X,Y) units code = " << inum << "!\n";
252         exit(-1);
253     }
254
255     // Units code; 2 represents meters as the unit of measure for
256     // elevation coordinates throughout the file.
257     inum = next_int();
258     if ( inum != 2 ) {
259         cout << "    Unknown (Z) units code = " << inum << "!\n";
260         exit(-1);
261     }
262
263     // Number (n) of sides in the polygon which defines the coverage of
264     // the DEM file (usually equal to 4).
265     inum = next_int();
266     if ( inum != 4 ) {
267         cout << "    Unknown polygon dimension = " << inum << "!\n";
268         exit(-1);
269     }
270
271     // Ground coordinates of bounding box in arc-seconds
272     dem_x1 = originx = next_exp();
273     dem_y1 = originy = next_exp();
274     cout << "    Origin = (" << originx << "," << originy << ")\n";
275
276     dem_x2 = next_exp();
277     dem_y2 = next_exp();
278
279     dem_x3 = next_exp();
280     dem_y3 = next_exp();
281
282     dem_x4 = next_exp();
283     dem_y4 = next_exp();
284
285     // Minimum/maximum elevations in meters
286     dem_z1 = next_exp();
287     dem_z2 = next_exp();
288     cout << "    Elevation range " << dem_z1 << dem_z2 << "\n";
289
290     // Counterclockwise angle from the primary axis of ground
291     // planimetric referenced to the primary axis of the DEM local
292     // reference system.
293     token = next_token();
294
295     // Accuracy code; 0 indicates that a record of accuracy does not
296     // exist and that no record type C will follow.
297
298     // DEM spacial resolution.  Usually (3,3,1) (3,6,1) or (3,9,1)
299     // depending on latitude
300
301     // I will eventually have to do something with this for data at
302     // higher latitudes */
303     token = next_token();
304     cout << "    accuracy & spacial resolution string = " + token + "\n";
305     i = token.length();
306     cout << "    length = " << i << "\n";
307
308     ptr = token.c_str() + i - 12;
309     cout << "    last field = " << ptr << " = " << atof(ptr) << "\n";
310     ptr[0] = '\0';
311
312     ptr = ptr - 12;
313     col_step = atof(ptr);
314     cout << "    last field = " << ptr << " = " << col_step << "\n";
315     ptr[0] = '\0';
316
317     ptr = ptr - 12;
318     row_step = atof(ptr);
319     cout << "    last field = " << ptr << " = " << row_step << "\n";
320     ptr[0] = '\0';
321
322     // accuracy code = atod(token)
323     ptr = ptr - 12;
324     inum = atoi(ptr);
325     cout << "    Accuracy code = " << inum << "\n";
326
327     cout << "    column step = " << col_step << 
328         "  row step = " << row_step << "\n";
329
330     // dimension of arrays to follow (1)
331     token = next_token();
332
333     // number of profiles
334     dem_num_profiles = cols = next_int();
335     cout << "    Expecting " << dem_num_profiles << " profiles\n";
336
337     return 1;
338 }
339
340
341 // read and parse DEM "B" record
342 void fgDEM::read_b_record( ) {
343     string token;
344     int i;
345
346     // row / column id of this profile
347     prof_row = next_int();
348     prof_col = next_int();
349     // printf("col id = %d  row id = %d\n", prof_col, prof_row);
350
351     // Number of columns and rows (elevations) in this profile
352     prof_num_rows = rows = next_int();
353     prof_num_cols = next_int();
354     // printf("    profile num rows = %d\n", prof_num_rows);
355
356     // Ground planimetric coordinates (arc-seconds) of the first
357     // elevation in the profile
358     prof_x1 = next_exp();
359     prof_y1 = next_exp();
360     // printf("    Starting at %.2f %.2f\n", prof_x1, prof_y1);
361
362     // Elevation of local datum for the profile.  Always zero for
363     // 1-degree DEM, the reference is mean sea level.
364     token = next_token();
365
366     // Minimum and maximum elevations for the profile.
367     token = next_token();
368     token = next_token();
369
370     // One (usually) dimensional array (prof_num_cols,1) of elevations
371     for ( i = 0; i < prof_num_rows; i++ ) {
372         prof_data = next_int();
373         dem_data[cur_col][i] = (float)prof_data;
374     }
375 }
376
377
378 // parse dem file
379 int fgDEM::parse( ) {
380     int i;
381
382     cur_col = 0;
383
384     if ( !read_a_record() ) {
385         return(0);
386     }
387
388     for ( i = 0; i < dem_num_profiles; i++ ) {
389         // printf("Ready to read next b record\n");
390         read_b_record();
391         cur_col++;
392
393         if ( cur_col % 100 == 0 ) {
394             cout << "    loaded " << cur_col << "profiles of data\n";
395         }
396     }
397
398     cout << "    Done parsing\n";
399
400     return 1;
401 }
402
403
404 // return the current altitude based on mesh data.  We should rewrite
405 // this to interpolate exact values, but for now this is good enough
406 double fgDEM::interpolate_altitude( double lon, double lat ) {
407     // we expect incoming (lon,lat) to be in arcsec for now
408
409     double xlocal, ylocal, dx, dy, zA, zB, elev;
410     int x1, x2, x3, y1, y2, y3;
411     float z1, z2, z3;
412     int xindex, yindex;
413
414     /* determine if we are in the lower triangle or the upper triangle 
415        ______
416        |   /|
417        |  / |
418        | /  |
419        |/   |
420        ------
421
422        then calculate our end points
423      */
424
425     xlocal = (lon - originx) / col_step;
426     ylocal = (lat - originy) / row_step;
427
428     xindex = (int)(xlocal);
429     yindex = (int)(ylocal);
430
431     // printf("xindex = %d  yindex = %d\n", xindex, yindex);
432
433     if ( xindex + 1 == cols ) {
434         xindex--;
435     }
436
437     if ( yindex + 1 == rows ) {
438         yindex--;
439     }
440
441     if ( (xindex < 0) || (xindex + 1 >= cols) ||
442          (yindex < 0) || (yindex + 1 >= rows) ) {
443         return(-9999);
444     }
445
446     dx = xlocal - xindex;
447     dy = ylocal - yindex;
448
449     if ( dx > dy ) {
450         // lower triangle
451         // printf("  Lower triangle\n");
452
453         x1 = xindex; 
454         y1 = yindex; 
455         z1 = dem_data[x1][y1];
456
457         x2 = xindex + 1; 
458         y2 = yindex; 
459         z2 = dem_data[x2][y2];
460                                   
461         x3 = xindex + 1; 
462         y3 = yindex + 1; 
463         z3 = dem_data[x3][y3];
464
465         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
466         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
467         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
468         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
469
470         zA = dx * (z2 - z1) + z1;
471         zB = dx * (z3 - z1) + z1;
472         
473         // printf("  zA = %.2f  zB = %.2f\n", zA, zB);
474
475         if ( dx > FG_EPSILON ) {
476             elev = dy * (zB - zA) / dx + zA;
477         } else {
478             elev = zA;
479         }
480     } else {
481         // upper triangle
482         // printf("  Upper triangle\n");
483
484         x1 = xindex; 
485         y1 = yindex; 
486         z1 = dem_data[x1][y1];
487
488         x2 = xindex; 
489         y2 = yindex + 1; 
490         z2 = dem_data[x2][y2];
491                                   
492         x3 = xindex + 1; 
493         y3 = yindex + 1; 
494         z3 = dem_data[x3][y3];
495
496         // printf("  dx = %.2f  dy = %.2f\n", dx, dy);
497         // printf("  (x1,y1,z1) = (%d,%d,%d)\n", x1, y1, z1);
498         // printf("  (x2,y2,z2) = (%d,%d,%d)\n", x2, y2, z2);
499         // printf("  (x3,y3,z3) = (%d,%d,%d)\n", x3, y3, z3);
500  
501         zA = dy * (z2 - z1) + z1;
502         zB = dy * (z3 - z1) + z1;
503         
504         // printf("  zA = %.2f  zB = %.2f\n", zA, zB );
505         // printf("  xB - xA = %.2f\n", col_step * dy / row_step);
506
507         if ( dy > FG_EPSILON ) {
508             elev = dx * (zB - zA) / dy    + zA;
509         } else {
510             elev = zA;
511         }
512     }
513
514     return(elev);
515 }
516
517
518 // Use least squares to fit a simpler data set to dem data
519 void fgDEM::fit( double error, fgBUCKET *p ) {
520     double x[DEM_SIZE_1], y[DEM_SIZE_1];
521     double m, b, ave_error, max_error;
522     double cury, lasty;
523     int n, row, start, end, good_fit;
524     int colmin, colmax, rowmin, rowmax;
525     // FILE *dem, *fit, *fit1;
526
527     printf("Initializing output mesh structure\n");
528     outputmesh_init();
529
530     // determine dimensions
531     colmin = p->x * ( (cols - 1) / 8);
532     colmax = colmin + ( (cols - 1) / 8);
533     rowmin = p->y * ( (rows - 1) / 8);
534     rowmax = rowmin + ( (rows - 1) / 8);
535     printf("Fitting region = %d,%d to %d,%d\n", colmin, rowmin, colmax, rowmax);
536     
537     // include the corners explicitly
538     outputmesh_set_pt(colmin, rowmin, dem_data[colmin][rowmin]);
539     outputmesh_set_pt(colmin, rowmax, dem_data[colmin][rowmax]);
540     outputmesh_set_pt(colmax, rowmax, dem_data[colmax][rowmax]);
541     outputmesh_set_pt(colmax, rowmin, dem_data[colmax][rowmin]);
542
543     printf("Beginning best fit procedure\n");
544
545     for ( row = rowmin; row <= rowmax; row++ ) {
546         // fit  = fopen("fit.dat",  "w");
547         // fit1 = fopen("fit1.dat", "w");
548
549         start = colmin;
550
551         // printf("    fitting row = %d\n", row);
552
553         while ( start < colmax ) {
554             end = start + 1;
555             good_fit = 1;
556
557             x[(end - start) - 1] = 0.0 + ( start * col_step );
558             y[(end - start) - 1] = dem_data[start][row];
559
560             while ( (end <= colmax) && good_fit ) {
561                 n = (end - start) + 1;
562                 // printf("Least square of first %d points\n", n);
563                 x[end - start] = 0.0 + ( end * col_step );
564                 y[end - start] = dem_data[end][row];
565                 least_squares(x, y, n, &m, &b);
566                 ave_error = least_squares_error(x, y, n, m, b);
567                 max_error = least_squares_max_error(x, y, n, m, b);
568
569                 /*
570                 printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n", 
571                 start, end, ave_error, max_error, m, b);
572                 
573                 f = fopen("gnuplot.dat", "w");
574                 for ( j = 0; j <= end; j++) {
575                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ), 
576                             dem_data[row][j]);
577                 }
578                 for ( j = start; j <= end; j++) {
579                     fprintf(f, "%.2f %.2f\n", 0.0 + ( j * col_step ), 
580                             dem_data[row][j]);
581                 }
582                 fclose(f);
583
584                 printf("Please hit return: "); gets(junk);
585                 */
586
587                 if ( max_error > error ) {
588                     good_fit = 0;
589                 }
590                 
591                 end++;
592             }
593
594             if ( !good_fit ) {
595                 // error exceeded the threshold, back up
596                 end -= 2;  // back "end" up to the last good enough fit
597                 n--;       // back "n" up appropriately too
598             } else {
599                 // we popped out of the above loop while still within
600                 // the error threshold, so we must be at the end of
601                 // the data set
602                 end--;
603             }
604             
605             least_squares(x, y, n, &m, &b);
606             ave_error = least_squares_error(x, y, n, m, b);
607             max_error = least_squares_max_error(x, y, n, m, b);
608
609             /*
610             printf("\n");
611             printf("%d - %d  ave error = %.2f  max error = %.2f  y = %.2f*x + %.2f\n", 
612                    start, end, ave_error, max_error, m, b);
613             printf("\n");
614
615             fprintf(fit1, "%.2f %.2f\n", x[0], m * x[0] + b);
616             fprintf(fit1, "%.2f %.2f\n", x[end-start], m * x[end-start] + b);
617             */
618
619             if ( start > colmin ) {
620                 // skip this for the first line segment
621                 cury = m * x[0] + b;
622                 outputmesh_set_pt(start, row, (lasty + cury) / 2);
623                 // fprintf(fit, "%.2f %.2f\n", x[0], (lasty + cury) / 2);
624             }
625
626             lasty = m * x[end-start] + b;
627             start = end;
628         }
629
630         /*
631         fclose(fit);
632         fclose(fit1);
633
634         dem = fopen("gnuplot.dat", "w");
635         for ( j = 0; j < DEM_SIZE_1; j++) {
636             fprintf(dem, "%.2f %.2f\n", 0.0 + ( j * col_step ), 
637                     dem_data[j][row]);
638         } 
639         fclose(dem);
640         */
641
642         // NOTICE, this is for testing only.  This instance of
643         // output_nodes should be removed.  It should be called only
644         // once at the end once all the nodes have been generated.
645         // newmesh_output_nodes(&nm, "mesh.node");
646         // printf("Please hit return: "); gets(junk);
647     }
648
649     // outputmesh_output_nodes(fg_root, p);
650 }
651
652
653 // Initialize output mesh structure
654 void fgDEM::outputmesh_init( void ) {
655     int i, j;
656     
657     for ( j = 0; j < DEM_SIZE_1; j++ ) {
658         for ( i = 0; i < DEM_SIZE_1; i++ ) {
659             output_data[i][j] = -9999.0;
660         }
661     }
662 }
663
664
665 // Get the value of a mesh node
666 double fgDEM::outputmesh_get_pt( int i, int j ) {
667     return ( output_data[i][j] );
668 }
669
670
671 // Set the value of a mesh node
672 void fgDEM::outputmesh_set_pt( int i, int j, double value ) {
673     // printf("Setting data[%d][%d] = %.2f\n", i, j, value);
674    output_data[i][j] = value;
675 }
676
677
678 // Write out a node file that can be used by the "triangle" program.
679 // Check for an optional "index.node.ex" file in case there is a .poly
680 // file to go along with this node file.  Include these nodes first
681 // since they are referenced by position from the .poly file.
682 void fgDEM::outputmesh_output_nodes( const string& fg_root, fgBUCKET *p ) {
683     double exnodes[MAX_EX_NODES][3];
684     double junk1, junk2, junk3;
685     struct stat stat_buf;
686     string dir;
687     char base_path[256], file[256], exfile[256];
688 #ifdef WIN32
689     char tmp_path[256];
690 #endif
691     string command;
692     FILE *fd;
693     long int index;
694     int colmin, colmax, rowmin, rowmax;
695     int i, j, count, excount, result;
696
697     // determine dimensions
698     colmin = p->x * ( (cols - 1) / 8);
699     colmax = colmin + ( (cols - 1) / 8);
700     rowmin = p->y * ( (rows - 1) / 8);
701     rowmax = rowmin + ( (rows - 1) / 8);
702     cout << "  dumping region = " << colmin << "," << rowmin << " to " <<
703         colmax << "," << rowmax << "\n";
704
705     // generate the base directory
706     fgBucketGenBasePath(p, base_path);
707     cout << "fg_root = " + fg_root + "  Base Path = " + base_path + "\n";
708     dir = fg_root + "/Scenery/" + base_path;
709     cout << "Dir = " + dir + "\n";
710     
711     // stat() directory and create if needed
712     result = stat(dir.c_str(), &stat_buf);
713     if ( result != 0 ) {
714         cout << "Stat error need to create directory\n";
715
716 #ifndef WIN32
717
718         command = "mkdir -p " + dir + "\n";
719         system( command.c_str() );
720
721 #else // WIN32
722
723         // Cygwin crashes when trying to output to node file
724         // explicitly making directory structure seems OK on Win95
725
726         extract_path (base_path, tmp_path);
727
728         dir = fg_root + "/Scenery";
729         if (my_mkdir ( dir.c_str() )) { exit (-1); }
730
731         dir += fg_root + "/Scenery/" + tmp_path;
732         if (my_mkdir ( dir.c_str() )) { exit (-1); }
733
734         dir += fg_root + "/Scenery/" + base_path;
735         if (my_mkdir ( dir.c_str() )) { exit (-1); }
736
737 #endif // WIN32
738
739     } else {
740         // assume directory exists
741     }
742
743     // get index and generate output file name
744     index = fgBucketGenIndex(p);
745     sprintf(file, "%s/%ld.node", dir.c_str(), index);
746
747     // get (optional) extra node file name (in case there is matching
748     // .poly file.
749     strcpy(exfile, file);
750     strcat(exfile, ".ex");
751
752     // load extra nodes if they exist
753     excount = 0;
754     if ( (fd = fopen(exfile, "r")) != NULL ) {
755         fscanf(fd, "%d %d %d %d", &excount, &junk1, &junk2, &junk3);
756
757         if ( excount > MAX_EX_NODES - 1 ) {
758             printf("Error, too many 'extra' nodes, increase array size\n");
759             exit(-1);
760         } else {
761             printf("    Expecting %d 'extra' nodes\n", excount);
762         }
763
764         for ( i = 1; i <= excount; i++ ) {
765             fscanf(fd, "%d %lf %lf %lf\n", &junk1, 
766                    &exnodes[i][0], &exnodes[i][1], &exnodes[i][2]);
767             printf("(extra) %d %.2f %.2f %.2f\n", 
768                     i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
769         }
770         fclose(fd);
771     }
772
773     printf("Creating node file:  %s\n", file);
774     fd = fopen(file, "w");
775
776     // first count regular nodes to generate header
777     count = 0;
778     for ( j = rowmin; j <= rowmax; j++ ) {
779         for ( i = colmin; i <= colmax; i++ ) {
780             if ( output_data[i][j] > -9000.0 ) {
781                 count++;
782             }
783         }
784         // printf("    count = %d\n", count);
785     }
786     fprintf(fd, "%d 2 1 0\n", count + excount);
787
788     // now write out extra node data
789     for ( i = 1; i <= excount; i++ ) {
790         fprintf(fd, "%d %.2f %.2f %.2f\n", 
791                 i, exnodes[i][0], exnodes[i][1], exnodes[i][2]);
792     }
793
794     // write out actual node data
795     count = excount + 1;
796     for ( j = rowmin; j <= rowmax; j++ ) {
797         for ( i = colmin; i <= colmax; i++ ) {
798             if ( output_data[i][j] > -9000.0 ) {
799                 fprintf(fd, "%d %.2f %.2f %.2f\n", 
800                         count++, 
801                         originx + (double)i * col_step, 
802                         originy + (double)j * row_step,
803                         output_data[i][j]);
804             }
805         }
806         // printf("    count = %d\n", count);
807     }
808
809     fclose(fd);
810 }
811
812
813 fgDEM::~fgDEM( void ) {
814     // printf("class fgDEM DEstructor called.\n");
815     delete(dem_data);
816     delete(output_data);
817 }
818
819
820 // $Log$
821 // Revision 1.14  1998/09/19 17:59:45  curt
822 // Use c++ streams (fg_gzifstream).  Also converted many character arrays to
823 // the string class.
824 //
825 // Revision 1.13  1998/09/09 16:24:04  curt
826 // Fixed a bug in the handling of exclude files which was causing
827 // a crash by calling fclose() on an invalid file handle.
828 //
829 // Revision 1.12  1998/08/24 20:03:31  curt
830 // Eliminated a possible memory overrun error.
831 // Use the proper free() rather than the incorrect delete().
832 //
833 // Revision 1.11  1998/07/20 12:46:11  curt
834 // When outputing to a .node file, first check for an optional
835 // "index.node.ex" file in case there is a .poly file to go along with this
836 // node file.  Include these nodes first since they are referenced by position
837 // from the .poly file.  This is my first pass at adding an area "cutout"
838 // feature to the terrain generation pipeline.
839 //
840 // Revision 1.10  1998/07/13 20:58:02  curt
841 // .
842 //
843 // Revision 1.9  1998/07/13 15:29:49  curt
844 // Added #ifdef HAVE_CONFIG_H
845 //
846 // Revision 1.8  1998/07/04 00:47:18  curt
847 // typedef'd struct fgBUCKET.
848 //
849 // Revision 1.7  1998/06/05 18:14:39  curt
850 // Abort out early when reading the "A" record if it doesn't look like
851 // a proper DEM file.
852 //
853 // Revision 1.6  1998/05/02 01:49:21  curt
854 // Fixed a bug where the wrong variable was being initialized.
855 //
856 // Revision 1.5  1998/04/25 15:00:32  curt
857 // Changed "r" to "rb" in gzopen() options.  This fixes bad behavior in win32.
858 //
859 // Revision 1.4  1998/04/22 13:14:46  curt
860 // Fixed a bug in zlib usage.
861 //
862 // Revision 1.3  1998/04/18 03:53:05  curt
863 // Added zlib support.
864 //
865 // Revision 1.2  1998/04/14 02:43:27  curt
866 // Used "new" to auto-allocate large DEM parsing arrays in class constructor.
867 //
868 // Revision 1.1  1998/04/08 22:57:22  curt
869 // Adopted Gnu automake/autoconf system.
870 //
871 // Revision 1.3  1998/04/06 21:09:41  curt
872 // Additional win32 support.
873 // Fixed a bad bug in dem file parsing that was causing the output to be
874 // flipped about x = y.
875 //
876 // Revision 1.2  1998/03/23 20:35:41  curt
877 // Updated to use FG_EPSILON
878 //
879 // Revision 1.1  1998/03/19 02:54:47  curt
880 // Reorganized into a class lib called fgDEM.
881 //
882 // Revision 1.1  1998/03/19 01:46:28  curt
883 // Initial revision.
884 //