/* dem2pgm.c copyright 1998 Mark Stock
 *
 * Original header:
 *
 * dem2tga.c (c) 1997 Jon Larimer
 *
 *  This is a program to convert USGS Digital Elevation Model (DEM) data
 *  into a .TGA file, viewable by many graphics applications.  This program
 *  is public domain, do whatever you want with it, except claim it as your
 *  own or sell it. If you find a way to make money from my work, I want a
 *  piece of the action.
 *
 *  This should compile with any ANSI standard compiler, including Unix
 *  and DOS versions of gcc and many other C compilers.  If you have any
 *  questions or comments email jonl@alltel.net.
 *
 * End original header
 *
 * Modified 1998-05-07 by Mark Stock, mstock@umich.edu
 *  -writes Poskanzer Portable Graymap files now
 *  -rewrote nextint, nextfloat routines, added nextdouble
 *  -took out elevation scaling feature, as PGMs go up to 16 bits
 *  -fixed reading of incomplete columns in 7.5' DEMs
 *  -sends a suggested command to stdout to rectify and scale the PGM
 *
 * Compile with  cc -o dem2pgm dem2pgm.c -lm
 */


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <sys/errno.h>

int writepgmheader(FILE *fptr, int r, int c, int max) {
    /* PGMs have a very simple header */
    fprintf(fptr,"P2\n%d %d\n%d\n",r,c,max);
    fflush(fptr);
    return 0;
}

int findint(FILE * fptr) {
    /* this function searches forward to find the next int,  
     * as opposed to nextint(), which takes the next 6 chars
     * and makes it an int. */
    char ch;
    int i = 0;
    char in[7];

    while(isspace(ch = fgetc(fptr))) { }
    (void)ungetc(ch,fptr);
    while(!isspace(ch = fgetc(fptr))) {
        in[i++] = ch;
    }
    (void)ungetc(' ',fptr);

    in[i] = '\0';
    return(atoi(in));
}

int nextint(FILE * fptr) {
    char ch;
    int i = 0;
    int count;
    char in[7] = "";

    for (count=0; count<6; count++) {
        ch = fgetc(fptr);
        if (isspace(ch)) {
            /* check for LF in 7.5' DEMs */
            if (ch == 10) { return(-1); }
        } else { 
            in[i++] = ch;
        }
    }
    /* for 1-degree DEMs, there's no LFs */
    if (strlen(in) == 0) return(-2);

    in[i] = '\0';
    return(atoi(in));
} 
          
double nextfloat(FILE * fptr) {
    char ch;
    int i = 0;
    int count;
    char in[13];

    for (count=0; count<12; count++) {
        ch = fgetc(fptr);
        if (!isspace(ch)) in[i++] = ch;
    }

    in[i] = '\0';
    atof(in);
}

double nextdouble(FILE * fptr) {
    char ch;
    int i = 0;

    int count;
    char in[25];

    for (count=0; count<24; count++) {
        ch = fgetc(fptr);
        if (!isspace(ch)) {
            if (ch == 'D') ch = 'E';
            in[i++] = ch;
        }
    }

    in[i] = '\0';
    atof(in);
}

int main(int argc, char **argv) {
    FILE * demfile;
    FILE * outfile;
    int i, j, k;
    int maxelev, minelev;
    int maxrows, rows, columns;
    int rowoffset;
    int r, c;
    int cur_r, cur_c;
    double row_ymin, ymin, yres, yval;
    double xres, xval;
    char name[144];
  
    fprintf(stdout, "\noriginal: dem2tga v1.0 (c) 1997 Jon Larimer\n");
    fprintf(stdout, "dem2pgm changes by Mark Stock, mstock@umich.edu\n");
    fprintf(stdout, "This program is public domain.\n\n");

    if(argc < 2) {
        fprintf(stdout, "usage: dem2pgm input_file output_file\n");
        fprintf(stdout, "  input_file    is the DEM you want to convert\n");
        fprintf(stdout, "  output_file   is the name you want to save the image under.\n");
        fprintf(stdout, "You should supply a .pgm to the end of the output file name.\n");
        exit(1);
    }

    if((demfile = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "%s: fopen: %s", argv[0], strerror(0));
        exit(1);
    }
  
    /****** DEM TYPE A RECORDS ******/
    /* get the name field (144 characters) */
    for (i=0; i<144; i++) { name[i] = fgetc(demfile); }
    name[i+1] = '\0';

    /* clean off the whitespace at the end */
    for (i=strlen(name)-2; i>0; i--) {
        if (!isspace(name[i])) i=0;
        else name[i] = '\0';
    }

    fprintf(stdout, "Quad name field: %s\n", name);

    /* don't need the next 19 items for anything */
    for (i=0; i<4; i++) (void)nextint(demfile);
    for (i=0; i<15; i++) (void)nextdouble(demfile);
    fprintf(stdout, "Units code (ground planametric coordinates): %i\n", nextint(demfile));
    fprintf(stdout, "Units code (elevation coordinates): %i\n", nextint(demfile));
    (void)nextint(demfile);

    printf("Ground coordinates of 4 corners of DEM: (in arc-seconds) \n");
    for (i=0; i<4; i++) {
        xval = nextdouble(demfile);
        yval = nextdouble(demfile);
        fprintf(stdout, "  %.4f  %.4f\n", xval, yval);
        if (i==0) ymin = yval;
        if (yval < ymin) ymin = yval;
    }

    minelev = (int)nextdouble(demfile);
    maxelev = (int)nextdouble(demfile);
    fprintf(stdout, "\nMinimum elevation: %i\n", minelev);
    fprintf(stdout, "Maximum elevation: %i\n", maxelev);

    /* only need one of the next items */
    (void)nextdouble(demfile);
    (void)nextint(demfile);
    xres = (double)nextfloat(demfile);
    yres = (double)nextfloat(demfile);
    (void)nextfloat(demfile);

    cur_r = nextint(demfile);
    columns = nextint(demfile);
    if ((columns-1)%300 == 0) maxrows = 1201; /* all 1-deg DEMs are 1201x1201 */
        /* but some alaskan DEMs are 601, or 301 columns wide */
    else maxrows = 470;	/* just to be safe, 7.5' DEMs are 463 rows high */

    /* DEM TYPE B RECORDS */
    fprintf(stdout, "Writing image...");
    fflush(stdout);
    outfile = fopen(argv[2], "wb+");

    /* we need to see how many rows of data before we write the header */
    /* FIND the next int, its not just the next 6 chars */
    cur_r = findint(demfile);
    cur_c = nextint(demfile);
    /* Note: this number of rows is not the total number of rows! */
    rows = nextint(demfile);
    /* writetgaheader(outfile, rows, columns); */
    writepgmheader(outfile, maxrows, columns, maxelev);

    /* now ready to write the image data */
    while(cur_c <= columns) {

        (void)nextint(demfile);
        (void)nextdouble(demfile);
        row_ymin = nextdouble(demfile);
        for (i=0; i<3; i++) (void)nextdouble(demfile);

        /* determine row offset if the data begins some rows up from 
         * the lower edge, but only if the file is not a 1-deg DEM */
        if ((columns-1)%300 != 0) {
            rowoffset = (int)(row_ymin - ymin)/yres;
            if (rowoffset > 0) {
                for (i=0; i<rowoffset; i++) fprintf(outfile,"0\n");
            }
        } else {
            rowoffset = 0;
        }

        /* read the data present */
        for (i=0; i<rows; i++) { 
            j = nextint(demfile); 
            if (j < 0) {
                if (j == -1) {
                    j = nextint(demfile);
                } else if (j == -2) {
                    j = findint(demfile);
                }
            }
            fprintf(outfile,"%d\n",j);
        }

        /* fill in all of the rest of the rows with zero */
        /* fprintf(stderr,"cur_c=%d: i=%d+%d; i<%d; i++\n",cur_c,rows,rowoffset,maxrows); */
        for (i=rows+rowoffset; i<maxrows; i++) fprintf(outfile,"0\n");

        if (cur_c < columns) {
            cur_r = findint(demfile);
            cur_c = nextint(demfile);
            rows = nextint(demfile);
        } else cur_c++;
    }

    fprintf(stdout, "\ndone.\n");
    fprintf(stdout, "Reccommend running\n");
    if ((columns-1)%300 == 0) {
        fprintf(stdout, "pnmscale -yscale %lf %s | pnmflip -ccw\n",cos(1.5708*(ymin+1800.0)/324000.0)*xres/yres,argv[2]);
    } else {
        fprintf(stdout, "pnmflip -ccw %s\n",argv[2]);
    }
    fclose(outfile);
    fclose(demfile);
    return 0;
}