/*
 * Read a 1-degree USGS Digital Elevation Model file (available on
 * ftp://edcftp.cr.usgs.gov/pub/data/DEM/250/), and convert to
 * various formats.
 *
 * This code assumes the DEM files are in what appears to be standard
 * format (1201x1201 grids.)
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <ctype.h>
#include <math.h>
#include <assert.h>

#define N_ROWS 1201
#define N_COLS 1201

typedef unsigned short Elevation;
#define ELEVATION_MAX 65535

typedef struct tagFileFormat FileFormat;
struct tagFileFormat {
    char* name;
    void (*init)(FILE* out, int width, int height);
    void (*write_row)(FILE* out, Elevation* row, int width);
};

static void Die(char* format, ...);
static int FindProfile(FILE* in, int row);
static int ScanInt(FILE* in, int* result);
static void ShortInit(FILE* out, int width, int height);
static void ShortWriteRow(FILE* out, Elevation* row, int width);
static void PGMInit(FILE* out, int width, int height);
static void PGMWriteRow(FILE* out, Elevation* row, int width);
static void PPMHFInit(FILE* out, int width, int height);
static void PPMHFWriteRow(FILE* out, Elevation* row, int width);

static FileFormat file_formats[] = {
    { "short", ShortInit, ShortWriteRow },
    { "pgm", PGMInit, PGMWriteRow },
    { "ppmhf", PPMHFInit, PPMHFWriteRow }
};
#define N_FORMATS (sizeof(file_formats)/sizeof(file_formats[0]))

int main(int argc, char** argv)
{
    int j, i, e, x, y, width, height, max_e, found_max_e;
    double scale;
    Elevation* data;
    Elevation* dp;
    Elevation* row;
    const char* backspace = "\b\b\b\b";
    FILE* in;
    FILE* out;
    FileFormat* ff;
    
    if (argc < 7 || argc > 9)
	Die("usage: %s format x y width height max-elevation ?dem-file? ?ppm-file?\n",
	    argv[0]);
    in = stdin;
    out = stdout;
    if (argc >= 8) {
	in = fopen(argv[7], "r");
	if (in == 0) Die("couldn't open DEM file \"%s\"\n", argv[7]);
    }
    if (argc == 9) {
	out = fopen(argv[8], "w");
	if (out == 0) Die("couldn't open PPM file \"%s\"\n", argv[8]);
    }
    for (i = 0; i < N_FORMATS; i++) {
	if (!strcmp(file_formats[i].name, argv[1]))
	    break;
    }
    if (i == N_FORMATS)
	Die("unknown file format \"%s\"\n", argv[1]);
    ff = &file_formats[i];
    x = atoi(argv[2]);
    y = atoi(argv[3]);
    width = atoi(argv[4]);
    height = atoi(argv[5]);
    if (!strcmp(argv[6], "auto")) {
	max_e = -1;
    } else {
	max_e = atoi(argv[6]);
    }
    data = (Elevation*) malloc(N_COLS * N_ROWS * sizeof(Elevation));
    row = (Elevation*) malloc(width * sizeof(Elevation));
    ff->init(out, width, height);
    dp = data;
    found_max_e = -1;
    for (j = 0; j < N_COLS; j++) {
	if (!FindProfile(in, j))
	    Die("couldn't find profile %d\n", j);
	for (i = 0; i < N_ROWS; i++) {
	    if (ScanInt(in, &e) != 1)
		Die("couldn't scan elevation, row = %d, column = %d\n", j, i);
	    if (e > found_max_e) found_max_e = e;
	    if (e > ELEVATION_MAX) e = ELEVATION_MAX;
	    else if (e < 0) e = 0;
	    *dp++ = e;
	}
	if (!(j & 0xf))
	    fprintf(stderr, "%4d%s", j, backspace);
    }
    if (max_e == -1)
	max_e = found_max_e;
    scale = ELEVATION_MAX / (double) max_e;
    fprintf(stderr, "\nMax elevation was %d, scale is %g\n", max_e, scale);
    fprintf(stderr, "Reversing rows & columns...");
    for (i = N_ROWS-y-1; i >= N_ROWS-(y+height); --i) {
	dp = data + i + x * N_COLS;
	for (j = 0; j < width; j++) {
	    e = (int) floor(scale * *dp);
	    if (e > ELEVATION_MAX) e = ELEVATION_MAX;
	    row[j] = e;
	    dp += N_COLS;
	}
	ff->write_row(out, row, width);
    }
    fclose(in);
    fclose(out);
    fprintf(stderr, "\n");
    exit(0);
}

static char* NextToken(FILE* in)
{
    static char buf[256];
    char* bp;
    int ch;
    /* skip whitespace */
    while ((ch = getc(in)) != EOF && isspace(ch))
	;
    if (ch == EOF) return 0;
    bp = buf;
    *bp++ = ch;
    /* copy until whitespace */
    while ((ch = getc(in)) != EOF && !isspace(ch))
	*bp++ = ch;
    if (ch != EOF) ungetc(ch, in);
    *bp = '\0';
    /* printf("NextToken() -> %s\n", buf); */
    return buf;
}

static int FindProfile(FILE* in, int row)
{
    char row_s[20], cols_s[20], one_s[2];
    char* token;
    char* looking_for;
    enum { ST_01, ST_row, ST_N_COLS, ST_11 } state;
    sprintf(row_s, "%d", row + 1);
    sprintf(cols_s, "%d", N_COLS);
    strcpy(one_s, "1");
    /* look for pattern "1 row N_COLS 1" */
    state = ST_01;
    looking_for = one_s;
    while ((token = NextToken(in)) != 0) {
	if (!strcmp(looking_for, token)) {
	    switch (state) {
	    case ST_01:
		looking_for = row_s;
		state = ST_row;
		break;
	    case ST_row:
		looking_for = cols_s;
		state = ST_N_COLS;
		break;
	    case ST_N_COLS:
		looking_for = one_s;
		state = ST_11;
		break;
	    case ST_11:
		goto found1;
	    }
	} else {
	    state = ST_01;
	    looking_for = one_s;
	}
    }
 found1:
    /* now skip coordinates, base elevation of profile, min & max elevations
     * of profile
     */
    if (!NextToken(in) || !NextToken(in) || !NextToken(in) ||
	!NextToken(in) || !NextToken(in))
	return 0;
    /* we should be at the start of the profile */
    return 1;
}

static void Die(char* format, ...)
{
    va_list args;
    va_start(args, format);
    vfprintf(stderr, format, args);
    va_end(args);
    exit(1);
}

static int ScanInt(FILE* in, int* result)
{
    int ch, n, sign;
    n = 0;
    while ((ch = getc(in)) != EOF && isspace(ch))
	;
    if (ch == EOF) return 0;
    if (ch == '-') {
	sign = -1;
    } else {
	ungetc(ch, in);
	sign = 1;
    }
    while ((ch = getc(in)) != EOF && isdigit(ch))
	n = ch - '0' + n * 10;
    if (ch != EOF) ungetc(ch, in);
    *result = sign * n;
    return 1;
}

/* conversion to 16-bit samples */

static void ShortInit(FILE* out, int width, int height)
{
}

static void ShortWriteRow(FILE* out, Elevation* row, int width)
{
    fwrite(row, sizeof(short), width, out);
}

/* conversion to 8-bit PGM */

static void PGMInit(FILE* out, int width, int height)
{
    fprintf(out, "P5\n%d %d\n255\n", width, height);
}

static void PGMWriteRow(FILE* out, Elevation* row, int width)
{
    unsigned char val;
    int i;
    for (i = 0; i < width; i++) {
	if (row[i] == 0) {
	    val = 0;
	} else {
	    val = ((row[i] - 1) >> 8) + 1;
	}
	putc(val, out);
    }
}

/* conversion to PPM stuffing elevation into red (high) and green (low)
 * bytes.  Using ppmtotga -rgb you can generate a TGA file suitable
 * for use as a POV height_field.
 */

static void PPMHFInit(FILE* out, int width, int height)
{
    fprintf(out, "P6\n%d %d\n255\n", width, height);
}

static void PPMHFWriteRow(FILE* out, Elevation* row, int width)
{
    int i;
    unsigned char val;
    unsigned short e;
    for (i = 0; i < width; i++) {
	e = row[i];
	val = (e >> 8) & 0xff;
	putc(val, out);
	val = e & 0xff;
	putc(val, out);
	val = 0;
	putc(val, out);
    }
}