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