/*
 * =========================================================================
 * drawmap - A program to draw maps using data from USGS geographic data files.
 * Copyright (c) 1997  Fred M. Erickson
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 * =========================================================================
 *
 *
 * Program to process 250K Digital Elevation Model (DEM),
 * 100K (optional-format) Digital Line Graph (DLG),
 * and Geographic Names Information System (GNIS)
 * files and produce colored maps in SUN Rasterfile format.
 *
 * At the time this program was written, DEM and DLG files were available
 * for free download from http://edcwww.cr.usgs.gov/ and GNIS files
 * were available from http://mapping.usgs.gov/
 */


#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <math.h>
#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <errno.h>
#include <time.h>
#include "raster.h"
#include "font_5x8.h"
#include "font_6x10.h"



#define VERSION "Version 1.3"



struct point  {
	double x;
	double y;
	struct point *point;
};


ssize_t buf_read(int, void *, size_t);
ssize_t buf_write(int, const void *, size_t);
ssize_t get_a_line(int, void *, size_t);
ssize_t buf_read_z(int, void *, size_t);
ssize_t get_a_line_z(int, void *, size_t);
double lat_conv(unsigned char *);
double lon_conv(unsigned char *);
void add_text(char *, long, long, long, unsigned char *, long, long, long, long);
void fill_area(double, double, long);
long get_factor(double);

void process_optional_dlg_file(int, int);
void draw_lines(struct point *, long);

/* #define COPYRIGHT_NAME	"Fred M. Erickson" */	/* Now defined in the Makefile */

#define DEM_SIZE	1201

/* The borders should be at least 60, if possible. */
#define TOP_BORDER	60
#define BOTTOM_BORDER	80
#define LEFT_BORDER	60
#define RIGHT_BORDER	60

#define	NUM_DEM	20	/* Number of DEM files allowed on input */
#define	NUM_DLG	400	/* Number of DLG files allowed on input */

#define HIGHEST_ELEVATION	32000	/* Elevation higher than any elevation expected in the DEM data */

#define SMOOTH_MAX	8	/* maximum radius of smoothing kernel */

#define BUF_SIZE	16384

#define OMIT_NEATLINES	1	/* If this is non-zero, then neatlines won't be drawn on the image. */


/*
 * Boundaries of image in latitude and longitude.
 * East Longitude is Positive, and North Latitude is positive.
 * West Longitude is Negative, and South Latitude is negative.
 */
double latitude_low;
double latitude_high;
double longitude_low;
double longitude_high;

/* double hlat = 0.0, hlong = 0.0, llat = 1.0e20, llong = 1.0e20; */

long zone;

long x, y;	/* width and height of output image (in pixels) */
long x_prime;

unsigned char *image_out;	/* Output image. */

/* long histogram[256]; */	/* For debugging. */
/* long angle_hist[100000]; */	/* For debugging. */
/* long total; */	/* For debugging. */

/*
 * Storage for attribute types.
 */
#define MAX_A_ATTRIB	100
#define MAX_L_ATTRIB	100
long num_A_attrib;
long num_L_attrib;
struct  {
	long major;
	long minor;
} attributes_A[MAX_A_ATTRIB];
struct  {
	long major;
	long minor;
} attributes_L[MAX_L_ATTRIB];

/*
 * Color map index values for the beginning of the various color blocks.
 *
 * The color map is built in blocks of 16.
 * In each block, the most intense version of a given color is first,
 * followed by decreasingly intense versions, with the last being black.
 * Thus any color plus 15 is black.  Any color + (some integer less than 15)
 * is a less intense version of that color.
 *
 * The 16th bank of 16 is reserved for bright, saturated colors, for use
 * in drawing roads, streams, and such.
 */
#define	GREEN		0
#define	YELLOW		16
#define	ORANGE		32	/* Really more like brown */
#define	RED		48
#define	MAGENTA		64
#define	CYAN		80
#define	L_GREEN		96
#define	L_YELLOW	112
#define	L_ORANGE	128
#define	L_RED		144
#define	L_MAGENTA	160
#define	L_BLUE		176
#define	L_CYAN		192
#define	WHITE		208
#define	UNUSED		224
#define COLOR_CHART_WIDTH	UNUSED	/* The color chart includes everything up through WHITE. */
#define	MISC		240
#define	BLUE		240
#define	B_RED		241
#define	B_GREEN		242
#define	B_BLUE		243
#define	B_CYAN		244
#define	B_MAGENTA	245
#define	B_YELLOW	246
#define	BLACK		GREEN + 15

#define BREAKPOINT        (M_PI/2.0)
#define BREAKINDEX        (3)

long round(double);
long max(long, long);
double max3(double, double, double);
double min3(double, double, double);

void
usage(char *program_name)
{
	fprintf(stderr, "\nDrawmap, %s.\n\n", VERSION);
	fprintf(stderr, "Usage:  %s [-o output_file.sun] [-l latitude1,longitude1,latitude2,longitude2] [-L]\n", program_name);
	fprintf(stderr, "          [-d dem_file1 [-d dem_file2 [...]]] [-g gnis_file] [-a attribute_file]\n");
	fprintf(stderr, "          [-x x_size] [-y y_size] [dlg_file1 [dlg_file2 [...]]]\n");
	fprintf(stderr, "\nNote that the DLG files are processed in order, and each one overlays the\n");
	fprintf(stderr, "last.  If you want (for example) roads on top of streams, put the\n");
	fprintf(stderr, "transportation data after the hydrography data.  Note also that\n");
	fprintf(stderr, "latitude/longitude values are in decimal degrees, and that east and north\n");
	fprintf(stderr, "are positive while west and south are negative.\n");
}



void
license(void)
{
	fprintf(stderr, "This program is free software; you can redistribute it and/or modify\n");
	fprintf(stderr, "it under the terms of the GNU General Public License as published by\n");
	fprintf(stderr, "the Free Software Foundation; either version 2, or (at your option)\n");
	fprintf(stderr, "any later version.\n\n");

	fprintf(stderr, "This program is distributed in the hope that it will be useful,\n");
	fprintf(stderr, "but WITHOUT ANY WARRANTY; without even the implied warranty of\n");
	fprintf(stderr, "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n");
	fprintf(stderr, "GNU General Public License for more details.\n\n");

	fprintf(stderr, "You should have received a copy of the GNU General Public License\n");
	fprintf(stderr, "along with this program; if not, write to the Free Software\n");
	fprintf(stderr, "Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.\n");
}



main(int argc, char *argv[])
{
	long i, j, k, l, m, n;
	int interp_size;
	long tick_width;
	double f, g;
	long file_index;
	long xx, yy;
	double red, green, blue;
	unsigned char a, b, c, d;
	long *lptr;
	long lsize;
	long smooth[SMOOTH_MAX + SMOOTH_MAX + 1][SMOOTH_MAX + SMOOTH_MAX + 1];
	long smooth_size;
	double tangent, tangent1, tangent2, tangent3;
	double fraction;
	double latitude;
	double longitude;
	long factor;
	long angle;
	long sum;
	long sum_count;
	struct rasterfile hdr;
	unsigned char  map[3][256];
	int dem_fdesc;
	int gnis_fdesc;
	int attribute_fdesc;
	int dlg_fdesc;
	int output_fdesc;
	ssize_t ret_val;
	long length;
	long start_x, start_y;
	unsigned char buf[BUF_SIZE];
	char *ptr;
	unsigned char *font;
	long font_width, font_height;
	time_t time_val;
	char latitude_code;
	char e_w_code;
	long location_code;
	unsigned char dem_name[100];
	long dem_flag;
	double latitude_dem, longitude_dem;
	double latitude1, longitude1, latitude2, longitude2;
	long tmp_width, tmp_height, tmp_x, tmp_y;
	char *dem_files[20];
	long num_dem, num_dlg;
	char *gnis_file;
	char *attribute_file;
	char *output_file;
	long option;
	long x_low, x_high, y_low, y_high;
	double reduction_y, reduction_xy;
	short *image_tmp;
	short *image_in;
	int gz_flag, lat_flag;
	long dem_size_x, dem_size_y;
	unsigned char ll_code[8];

	if (argc == 1)  {
		usage(argv[0]);
		exit(0);
	}


	/* Process arguments */
	x = -1;
	y = -1;
	latitude_low = 91.0;
	longitude_low = 181.0;
	latitude_high = -91.0;
	longitude_high = -181.0;
	gnis_file = (char *)0;
	attribute_file = (char *)0;
	output_file = (char *)0;
	num_dem = 0;
	dem_flag = 0;	/* When set to 1, this flag says that at least some DEM data was read in. */
	lat_flag = 0;	/* When set to 1, this flag says that the user explicitly specified the map boundaries. */
	opterr = 0;	/* Shut off automatic unrecognized-argument messages. */

	while ((option = getopt(argc, argv, "o:d:g:a:x:y:l:L")) != EOF)  {
		switch(option)  {
		case 'o':
			if (output_file != (char *)0)  {
				fprintf(stderr, "More than one output file specified with -o\n");
				usage(argv[0]);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No output file specified\n");
				usage(argv[0]);
				exit(0);
			}
			output_file = optarg;
			break;
		case 'd':
			if (num_dem >= NUM_DEM)  {
				fprintf(stderr, "Out of storage for DEM file names (max %d)\n", NUM_DEM);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No DEM file specified with -d\n");
				usage(argv[0]);
				exit(0);
			}
			dem_files[num_dem++] = optarg;
			break;
		case 'g':
			if (gnis_file != (char *)0)  {
				fprintf(stderr, "More than one GNIS file specified\n");
				usage(argv[0]);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No GNIS file specified with -g\n");
				usage(argv[0]);
				exit(0);
			}
			gnis_file = optarg;
			break;
		case 'a':
			if (attribute_file != (char *)0)  {
				fprintf(stderr, "More than one attribute file specified\n");
				usage(argv[0]);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No attribute file specified with -a\n");
				usage(argv[0]);
				exit(0);
			}
			attribute_file = optarg;
			break;
		case 'x':
			if (x >= 0)  {
				fprintf(stderr, "More than one -x value specified\n");
				usage(argv[0]);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No value specified with -x\n");
				usage(argv[0]);
				exit(0);
			}
			x = atoi(optarg);
			break;
		case 'y':
			if (y >= 0)  {
				fprintf(stderr, "More than one -y value specified\n");
				usage(argv[0]);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No value specified with -y\n");
				usage(argv[0]);
				exit(0);
			}
			y = atoi(optarg);
			break;
		case 'l':
			if ((latitude_low != 91.0) || (longitude_low != 181.0) ||
			    (latitude_high != -91.0) || (longitude_high != -181.0))  {
				fprintf(stderr, "More than one set of -l values specified\n");
				usage(argv[0]);
				exit(0);
			}
			if (optarg == (char *)0)  {
				fprintf(stderr, "No values specified with -l\n");
				usage(argv[0]);
				exit(0);
			}
			ptr = optarg;
			if (*ptr != '\0')  {
				latitude_low = strtod(ptr, &ptr);
			}
			ptr++;
			if (*ptr != '\0')  {
				longitude_low = strtod(ptr, &ptr);
			}
			ptr++;
			if (*ptr != '\0')  {
				latitude_high = strtod(ptr, &ptr);
			}
			ptr++;
			if (*ptr != '\0')  {
				longitude_high = strtod(ptr, &ptr);
			}
			if ((latitude_low == 91.0) || (longitude_low == 181.0) ||
			    (latitude_high == -91.0) || (longitude_high == -181.0))  {
				fprintf(stderr, "Incomplete set of -l values specified\n");
				usage(argv[0]);
				exit(0);
			}
			if ((latitude_low < -90.0) || (latitude_low > 90.0) || (latitude_high < -90.0) || (latitude_high > 90.0))  {
				fprintf(stderr, "Latitude must fall between -90 and 90 degrees, inclusive\n");
				usage(argv[0]);
				exit(0);
			}
			if ((longitude_low < -180.0) || (longitude_low > 180.0) || (longitude_high < -180.0) || (longitude_high > 180.0))  {
				fprintf(stderr, "Longitude must fall between -180 and 180 degrees, inclusive\n");
				usage(argv[0]);
				exit(0);
			}
			if (latitude_low > latitude_high)  {
				f = latitude_low;
				latitude_low = latitude_high;
				latitude_high = f;
			}
			if (longitude_low > longitude_high)  {
				f = longitude_low;
				longitude_low = longitude_high;
				longitude_high = f;
			}
			lat_flag = 1;
			break;
		case 'L':
			license();
			exit(0);
			break;
		default:
			usage(argv[0]);
			exit(0);
			break;
		}
	}
	num_dlg = argc - optind;

	/* Clean up the options. */
	if (output_file == (char *)0)  {
		output_file = "drawmap.sun";
	}
	if (x < 0)  {
		/* The user didn't specify an x value.  Provide one that is half of full resolution. */
		if (lat_flag != 0)  {
			x = round(0.5 * (longitude_high - longitude_low) * (double)(DEM_SIZE - 1));
		}
		else  {
			x = (DEM_SIZE - 1) >> 1;
		}
	}
	if (y < 0)  {
		/* The user didn't specify a y value.  Provide one that is half of full resolution. */
		if (lat_flag != 0)  {
			y = round(0.5 * (latitude_high - latitude_low) * (double)(DEM_SIZE - 1));
		}
		else  {
			y = (DEM_SIZE - 1) >> 1;
		}
	}
	if ((num_dem != 1) && (lat_flag == 0))  {
		fprintf(stderr, "The -l option is required unless there is exactly one -d option given.\n");
		usage(argv[0]);
		exit(0);
	}


	/* If we can't create the output file, then there is no point in continuing. */
	if ((output_fdesc = open(output_file, O_WRONLY | O_CREAT | O_TRUNC, 0644)) < 0)  {
		fprintf(stderr, "Can't create %s for writing, errno = %d\n", output_file, errno);
		exit(0);
	}


	/* If an attribute file was specified, then parse it now. */
	num_A_attrib = 0;
	num_L_attrib = 0;
	if (attribute_file != (char *)0)  {
		if (strcmp(attribute_file + strlen(attribute_file) - 3, ".gz") == 0)  {
			gz_flag = 1;
			if ((attribute_fdesc = buf_open_z(attribute_file, O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", attribute_file, errno);
				exit(0);
			}
		}
		else  {
			gz_flag = 0;
			if ((attribute_fdesc = buf_open(attribute_file, O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", attribute_file, errno);
				exit(0);
			}
		}

		fprintf(stderr, "Processing Attribute file:  %s\n", attribute_file);

		while ( 1 )  {
			if (gz_flag == 0)  {
				if ((ret_val = get_a_line(attribute_fdesc, buf, BUF_SIZE)) <= 0)  {
					break;
				}
			}
			else  {
				if ((ret_val = get_a_line_z(attribute_fdesc, buf, BUF_SIZE)) <= 0)  {
					break;
				}
			}

			buf[ret_val - 1] = '\0';	/* Put a null in place of the newline */

			switch(buf[0])  {
			case '\0':
			case '\n':
			case ' ':
			case '\t':
				/* Blank line, or line that begins with white space.  Ignore. */
				break;
			case '#':
				/* Comment line.  Ignore. */
				break;
			case 'N':
				/* We don't currently use Node attributes, so do nothing with them. */
				fprintf(stderr, "Ignoring Node attribute:  %s\n", buf);
				break;
			case 'A':
				/* Area attribute. */
				if (num_A_attrib >= MAX_A_ATTRIB)  {
					fprintf(stderr, "Out of space for Area attributes, ignoring:  %s\n", buf);
					break;
				}
				attributes_A[num_A_attrib].major = strtol(&buf[1], &ptr, 10);
				attributes_A[num_A_attrib].minor = strtol(ptr, &ptr, 10);
				num_A_attrib++;
				break;
			case 'L':
				/* Line attribute. */
				if (num_L_attrib >= MAX_L_ATTRIB)  {
					fprintf(stderr, "Out of space for Line attributes, ignoring:  %s\n", buf);
					break;
				}
				attributes_L[num_L_attrib].major = strtol(&buf[1], &ptr, 10);
				attributes_L[num_L_attrib].minor = strtol(ptr, &ptr, 10);
				num_L_attrib++;
				break;
			default:
				fprintf(stderr, "Ignoring unknown attribute type:  %s\n", buf);
				break;
			}
		}

		if (gz_flag == 0)  {
			buf_close(attribute_fdesc);
		}
		else  {
			buf_close_z(attribute_fdesc);
		}
	}


	/*
	 * Get memory for the DEM data.
	 *
	 * As all of the DEM files are read in, their data
	 * eventually get combined into this storage area.
	 * On the way, the data may get cropped, smoothed, or
	 * subsampled.
	 */
	image_in = (short *)malloc(2 * (y + 1) * (x + 1));
	if (image_in == (short *)0)  {
		fprintf(stderr, "malloc of image_in failed\n");
		exit(0);
	}


	/* Initialize SUN rasterfile header. */
	hdr.magic = MAGIC;
	hdr.width = x + LEFT_BORDER + RIGHT_BORDER;
	hdr.height = y + TOP_BORDER + BOTTOM_BORDER;
	hdr.depth = 8;
	hdr.length = (x + LEFT_BORDER + RIGHT_BORDER) * (y + TOP_BORDER + BOTTOM_BORDER);
	hdr.type = STANDARD;
	hdr.maptype = EQUAL_RGB;
	hdr.maplength = 768;

	/*
	 * Set up the color map.  A maximum of 15 colors (to denote elevation), each
	 * of which is modulated by 16 brightness levels (to denote topographical gradient).
	 *
	 * In general, a spiral on the CIE chromaticity diagram is considered
	 * a good way to encode sequential data.  It begins with a relatively
	 * saturated version of a color and spirals inward to reach brighter, whiter
	 * colors.  Since relief maps generally start with a greenish color at sea
	 * level, the spiral should go something like:  green, yellow, orange, red,
	 * magenta, blue, cyan, lighter green, lighter yellow, ligher orange, lighter
	 * red, lighter magenta, lighter blue, lighter cyan, white.
	 *
	 * After trying the above color sequence, it appears prudent to eliminate blue
	 * from the rotation, because it makes it look like there are big lakes
	 * everywhere.
	 * 
	 * CMYK printers and computer monitors roughly (but not completely) coincide
	 * in the portion of the CIE diagram that they occupy.  Their region occupies
	 * the center of the diagram, but doesn't extend to the edges.  Thus, if we
	 * want to be able to display and print on a wide variety of hardware,
	 * fully-saturated colors are not a good idea for the outer portions
	 * of the spiral.  For elevations below sea level, it seems reasonable to run
	 * the spiral backwards through cyan, blue, magenta, and black, as the
	 * elevations proceed from 0 to the trench bottoms.  For now, we aren't going
	 * to have any really deep depths, so we will just use cyan for any elevations
	 * below zero.
	 *
	 * Since most Americans prefer to deal in feet, we will use 1000-foot
	 * intervals, which are pretty close to 300-meter intervals.  The DEM
	 * data is in meters.
	 *
	 * We choose RGB values for each color in multiples of 15,
	 * because they get divided by 15 to fill out the color map.
	 *
	 *
	 *       Feet		Color		[R,G,B]			Meters
	 *   --------------	----------	-------------		----------------------
	 *         Below 0	Cyan		[ 60,255,255]		(center of Earth to 0)
	 *       0 to 1000	Green		[ 60,255, 60]		(   0 to  305)
	 *    1000 to 2000	Yellow		[255,255, 60]		( 305 to  610)
	 *    2000 to 3000	Orange		[255,165, 60]		( 610 to  914)
	 *    3000 to 4000	Red		[255, 60, 60]		( 914 to 1219)
	 *    4000 to 5000	Magenta		[255, 60,255]		(1219 to 1524)
	 *    5000 to 6000	Cyan		[ 60,255,255]		(1524 to 1829)
	 *    6000 to 7000	Light Green	[165,255,165]		(1829 to 2134)
	 *    7000 to 8000	Light Yellow	[255,255,165]		(2134 to 2438)
	 *    8000 to 9000 	Light Orange	[255,210,165]		(2438 to 2743)
	 *    9000 to 10000	Light Red	[255,165,165]		(2743 to 3048)
	 *   10000 to 11000	Light Magenta	[255,165,255]		(3048 to 3353)
	 *   11000 to 12000	Light Blue	[165,165,255]		(3353 to 3658)
	 *   12000 to 13000	Light Cyan	[165,255,255]		(3658 to 3962)
	 *     Above 13000	White		[255,255,255]		(4267 to infinity)
	 */
	map[0][GREEN] = 60;
	map[1][GREEN] = 255;
	map[2][GREEN] = 60;
	map[0][YELLOW] = 255;
	map[1][YELLOW] = 255;
	map[2][YELLOW] = 60;
	map[0][ORANGE] = 255;
	map[1][ORANGE] = 165;
	map[2][ORANGE] = 60;
	map[0][RED] = 255;
	map[1][RED] = 60;
	map[2][RED] = 60;
	map[0][MAGENTA] = 255;
	map[1][MAGENTA] = 60;
	map[2][MAGENTA] = 255;
	map[0][CYAN] = 60;
	map[1][CYAN] = 255;
	map[2][CYAN] = 255;
	map[0][L_GREEN] = 165;
	map[1][L_GREEN] = 255;
	map[2][L_GREEN] = 165;
	map[0][L_YELLOW] = 255;
	map[1][L_YELLOW] = 255;
	map[2][L_YELLOW] = 165;
	map[0][L_ORANGE] = 255;
	map[1][L_ORANGE] = 210;
	map[2][L_ORANGE] = 165;
	map[0][L_RED] = 255;
	map[1][L_RED] = 165;
	map[2][L_RED] = 165;
	map[0][L_MAGENTA] = 255;
	map[1][L_MAGENTA] = 165;
	map[2][L_MAGENTA] = 255;
	map[0][L_BLUE] = 165;
	map[1][L_BLUE] = 165;
	map[2][L_BLUE] = 255;
	map[0][L_CYAN] = 165;
	map[1][L_CYAN] = 255;
	map[2][L_CYAN] = 255;
	map[0][WHITE] = 255;
	map[1][WHITE] = 255;
	map[2][WHITE] = 255;
	/* Put black into the miscellaneous part of the table, for later update. */
	map[0][MISC] = 0;
	map[1][MISC] = 0;
	map[2][MISC] = 0;
	/* Put black into the rest of the table, which is unused for now. */
	map[0][UNUSED] = 0;
	map[1][UNUSED] = 0;
	map[2][UNUSED] = 0;

	/*
	 * We have the most intense color values inserted into the table.
	 * Now insert decreasingly less intense versions of each color.
	 * Each color decreases to black.
	 */
	for (i = 0; i < 16; i++)  {
		red = (double)map[0][i * 16] / 15.0;
		blue = (double)map[1][i * 16] / 15.0;
		green = (double)map[2][i * 16] / 15.0;

		for (j = 1; j < 15; j++)  {
			map[0][(i * 16) + j] = map[0][i * 16] - (unsigned char)round(((double)j * red));
			map[1][(i * 16) + j] = map[1][i * 16] - (unsigned char)round(((double)j * blue));
			map[2][(i * 16) + j] = map[2][i * 16] - (unsigned char)round(((double)j * green));
		}
		map[0][(i * 16) + 15] = 0;
		map[1][(i * 16) + 15] = 0;
		map[2][(i * 16) + 15] = 0;
	}

	/* Generate miscellaneous colors. */
	map[0][BLUE] = 60;
	map[1][BLUE] = 60;
	map[2][BLUE] = 255;
	map[0][B_RED] = 255;
	map[1][B_RED] = 0;
	map[2][B_RED] = 0;
	map[0][B_GREEN] = 0;
	map[1][B_GREEN] = 255;
	map[2][B_GREEN] = 0;
	map[0][B_BLUE] = 0;
	map[1][B_BLUE] = 0;
	map[2][B_BLUE] = 255;
	map[0][B_CYAN] = 0;
	map[1][B_CYAN] = 255;
	map[2][B_CYAN] = 255;
	map[0][B_MAGENTA] = 255;
	map[1][B_MAGENTA] = 0;
	map[2][B_MAGENTA] = 255;
	map[0][B_YELLOW] = 255;
	map[1][B_YELLOW] = 255;
	map[2][B_YELLOW] = 0;

	/*
	 * Write SUN rasterfile header and color map.
	 * My X86 Linux machine requires some swabbing.
	 * Your mileage may vary.
	 */
	lsize = sizeof(struct rasterfile) / 4;
	lptr = (long *)&hdr;
	for (i = 0; i < lsize; i++)  {
		a = ((*lptr) >> 24) & 0xff;
		b = ((*lptr) >> 16) & 0xff;
		c = ((*lptr) >> 8) & 0xff;
		d = (*lptr) & 0xff;

		*lptr = d << 8;
		*lptr = (*lptr | c) << 8;
		*lptr = (*lptr | b) << 8;
		*lptr = *lptr | a;

		lptr++;
	}
	write(output_fdesc, &hdr, sizeof(struct rasterfile));
	write(output_fdesc, map, sizeof(map));


	/*
	 * Before we begin processing map data, here is a short lecture on the
	 * Universal Transverse Mercator (UTM) coordinate system, which is commonly
	 * used in topographical maps, and by the military (it has been adopted by
	 * NATO and is used by the US military for ground operations).  UTM coordinates
	 * take the place of latitude and longitude, which can be cumbersome to deal
	 * with in the context of a small-area map.
	 *
	 * (UTM coordinates are used in the optional-format DLG files,
	 * and there is some reference to them in the DEM files.  GNIS files
	 * use latitude and longitude, in DDDMMSS format.)
	 *
	 * The UTM system begins by dividing the earth into 60 zones, each of
	 * which represents a slice (like a colored panel in a beach ball) that
	 * spans 6 degrees of longitude.  Zone 1 runs from 180 degrees West
	 * Longitude to 174 degrees West Longitude.  Zone 2 runs from 175W to
	 * 168W.  Zone 60 runs from 174E to 180E.
	 *
	 * UTM is only used from 84N to 80S.  At the poles, the Universal Polar
	 * Stereographic (UPS) projection is used.
	 *
	 * In each zone, points are represented by rectangular (x,y) coordinates
	 * that give distances, in meters, from the zone reference point.  This
	 * reference point is at the intersection of the Equator and the Central
	 * Meridian (the longitude line that runs down the center of the zone).
	 * The (x,y) coordinates give the distance in meters to the east and north
	 * of the reference point.
	 *
	 * In order to avoid having negative values for the UTM coordinates,
	 * some adjustments are made.  In the northern hemisphere, the y
	 * coordinate is simply measured from zero at the equator, but the
	 * Central Meridian is assigned a value of 500,000 meters (called a
	 * false easting), meaning that the distance (to the east) of a
	 * given point in the zone is the UTM x coordinate minus 500,000.
	 * In the southern hemisphere, the Central Meridian is again assigned
	 * a false easting of 500,000 meters, but the equator is no longer
	 * assigned a value of 0, but rather is assigned a value of 10,000,000
	 * meters north (called a false northing).
	 *
	 * When the zone is further divided into 4-degree latitude bands, the
	 * resulting pseudo-rectangular areas are called UTM grids.
	 *
	 * Note that a Mercator projection can be visualized by imagining
	 * a cylinder, sitting on one of its ends, with a globe inside.
	 * If a light is shined from, say, the center of the globe, the longitude
	 * lines will be projected onto the cylinder as vertical lines, and the
	 * latitude lines will be projected as circles around the cylinder.
	 * The longitude lines will be evenly spaced, but the latitude lines
	 * will be farther apart as the latitude increases.  One advantage
	 * of this projection is that it is conformal, meaning that angles and
	 * shapes are preserved during the transformation, for any given small
	 * region of the map.
	 *
	 * The Transverse Mercator projection is the same deal, except that the
	 * cylinder is tipped on its side and the desired Central Meridian is
	 * aligned so that it is vertical and tangent to the cylinder wall.
	 * Because of this orientation, shapes and areas near the Central
	 * Meridian are preserved, while shapes and areas distant from it
	 * are less accurate, especially when the top and/or bottom of the map
	 * is close to one of the poles so that the zone slice must be considerably
	 * stretched to form a rectangular grid.  Within a given UTM zone, however,
	 * the distortion is small.
	 *
	 * UTM is a Transverse Mercator projection, standardized for international
	 * use.
	 */


	/*
	 * Before reading in the DEM data, initialize the entire image to
	 * HIGHEST_ELEVATION, which will eventually be translated to the color WHITE.
	 * This is because we don't know, in advance, which parts of the image will
	 * be covered with data from the various DEM files.
	 */
	for (i = 0; i <= y; i++)  {
		for (j = 0; j <= x; j++)  {
			*(image_in + i * (x + 1) + j) = HIGHEST_ELEVATION;
		}
	}


	/*
	 * Read in DEM data, and rotate so North is on top.
	 * DEM files have a lot of header information, which we throw away,
	 * except for the code that tells us the latitude and longitude of the
	 * block, the x,y sample size of the block, and the name of the DEM block.
	 *
	 * Before getting started, malloc() some memory for temporary storage.
	 * Each DEM file is read into this temporary area, and then transferred to the
	 * appropriate place in the image_in storage area, after any needed
	 * smoothing, cropping, or subsampling.
	 */
	image_tmp = (short *)malloc(2 * DEM_SIZE * DEM_SIZE);
	if (image_tmp == (short *)0)  {
		fprintf(stderr, "malloc of image_tmp failed\n");
		exit(0);
	}
	dem_name[0] = '\0';
	file_index = 0;
	while (file_index < num_dem)  {
		if (strcmp(&dem_files[file_index][strlen(dem_files[file_index]) - 3], ".gz") == 0)  {
			gz_flag = 1;
			if ((dem_fdesc = buf_open_z(dem_files[file_index], O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", dem_files[file_index], errno);
				exit(0);
			}
		}
		else  {
			gz_flag = 0;
			if ((dem_fdesc = open(dem_files[file_index], O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", dem_files[file_index], errno);
				exit(0);
			}
		}

		fprintf(stderr, "Processing DEM file:  %s\n", dem_files[file_index]);

		dem_flag = 1;

		file_index++;

		if (gz_flag == 0)  {
			if ((ret_val = read(dem_fdesc, buf, 1024)) != 1024)  {	/* First strip the big header from the front */
				fprintf(stderr, "read from DEM file returns %d, expected 1024\n", ret_val);
				exit(0);
			}
		}
		else  {
			if ((ret_val = buf_read_z(dem_fdesc, buf, 1024)) != 1024)  {	/* First strip the big header from the front */
				fprintf(stderr, "zread from DEM file returns %d, expected 1024\n", ret_val);
				exit(0);
			}
		}

		/*
		 * Apparently, when the USGS digitized the DEM data, they didn't decide
		 * on a consistent format for the file headers.  Some files have the
		 * latitude/longitude code (described in a later comment) at column 49
		 * in the header, some files have it right at the beginning of the header,
		 * and others may have it somewhere else (although I have no examples of
		 * such at this time).
		 *
		 * It appears that the code always appears in the first 144-byte record,
		 * so we will just search the whole thing for something that looks like the
		 * code.  This is a pain, but it should be reliable.
		 */
		for (i = 0; i < 137; i++)  {
			if ((buf[i] != 'N') && (buf[i] != 'S'))  {
				continue;
			}

			if ((buf[i + 1] < 'A') || (buf[i + 1] > 'Z') ||
			    (buf[i + 2] < '0') || (buf[i + 2] > '9') ||
			    (buf[i + 3] < '0') || (buf[i + 3] > '9') ||
			    (buf[i + 4] != '-')  ||
			    (buf[i + 5] < '0') || (buf[i + 5] > '9') ||
			    (buf[i + 6] < '0') || (buf[i + 6] > '9') ||
			    ((buf[i + 7] != 'E') && (buf[i + 7] != 'W')))  {
				continue;
			}

			strncpy(ll_code, &buf[i], 8);
			strncpy(&buf[i], "        ", 8);
		}

	
		/*
		 * Save the name of the DEM block for later use.
		 * Strip out excessive spaces.
		 */
		if (dem_name[0] == '\0')  {
			i = 0;

			for (j = 0; j < 144; j++)  {
				if (buf[j] != ' ')  {
					/* If the character is not a space, just copy it. */
					dem_name[i++] = buf[j];
				}
				else  {
					/* Allow a maximum of two spaces in a row */
					if (i < 3)  continue;

					if ((dem_name[i - 2] != ' ') || (dem_name[i - 1] != ' '))  {
						dem_name[i++] = ' ';
					}
				}
			}

			dem_name[i] = '\0';
		}
		else  {
			strcpy(dem_name, "Map area spans multiple 1 degree blocks");
		}
	
		/*
		 * The longitude/latitude code in a DEM file is cryptic, and apparently
		 * is the name of corresponding map sheet.  It basically encodes the description
		 * of a UTM grid.  It takes a form that is illustrated by the following example:
		 *
		 *        NL12-08W
		 *
		 * where I think the N simply means "Northern Hemisphere".
		 * The L12 is a code that gives a 4 degree by 6 degree block.
		 * Starting at the equator, with 'A', the letter represents 4 degree
		 * chunks of latitude.  Thus 'L' represents the block from 44N to 48N.
		 * The 12 is the UTM zone number.  The calculation "-186 + (6 * zone)"
		 * gives the lower longitude of a 6 degree zone.  Thus, zone 12 represents
		 * longitudes from -114 to -108 (108W to 114W).
		 *
		 * The 4 degree by 6 degree block is divided into 12 rectangular areas,
		 * each of which contains an east and west chunk.  (The W tells us that
		 * this is the western 1 degree by 1 degree block.)  The areas are numbered
		 * as follows:
		 *
		 *       1       2        3
		 *       4       5        6
		 *       7       8        9
		 *      10      11       12
		 *
		 * In this case, area 08 is at 45N-46N and 110W-112W,
		 * so NL12-08W is at 45N-46N and 111W-112W.
		 * Area 1 defines the highest-latitude, highest-longitude
		 * block which, in this case, spans 47N-48N and 112W-114W.
		 * Area 12 defines the lowest-latitude, lowest-longitude
		 * block which, in this case, spans 44N-46N and 108W-110W.
		 */
		latitude_code = ll_code[1];
		zone = strtol(&ll_code[2], (char **)0, 10);
		location_code = strtol(&ll_code[5], (char **)0, 10);
		e_w_code = ll_code[7];
		latitude_dem = (double)((latitude_code - 'A') * 4);
		longitude_dem = -186.0 + (double)(zone * 6);
		i = (location_code - 1) / 3;
		j = (location_code + 2) % 3;
		latitude_dem = latitude_dem + 3.0 - (double)i;
		longitude_dem = longitude_dem + (double)j * 2.0 + (e_w_code == 'W' ? 0.0 : 1.0);

		/* If user did not provide the -l option, then initialize image boundary specifications. */
		if (lat_flag == 0)  {
			if (latitude_low > latitude_dem)  {
				latitude_low = latitude_dem;
			}
			if (longitude_low > longitude_dem)  {
				longitude_low = longitude_dem;
			}
			if (latitude_high < (latitude_dem + 1.0))  {
				latitude_high = latitude_dem + 1.0;
			}
			if (longitude_high < (longitude_dem + 1.0))  {
				longitude_high = longitude_dem + 1.0;
			}
		}

		/*
		 * If the DEM data don't overlap the image, then ignore them.
		 */
		if ((latitude_dem >= latitude_high) || ((latitude_dem + 1.0) <= latitude_low) ||
		    (longitude_dem >= longitude_high) || ((longitude_dem + 1.0) <= longitude_low))  {
			continue;
		}

		/*
		 * Get the number of columns in the data set.
		 * For all states in the USA, except Alaska, this value should be 1201.
		 * In Alaska, it can be 401 or 601.
		 *
		 * We use the number of columns to calculate an interpolation step size,
		 * which will be used to interpolate to fill out the dataset to 1201 by 1201 samples.
		 */
		dem_size_x = strtol(&buf[858], (char **)0, 10);
		if ((dem_size_x != 401) && (dem_size_x != 601) && (dem_size_x != 1201))  {
			fprintf(stderr, "Unexpected number of south-north columns in DEM data: %d\n", dem_size_x);
			exit(0);
		}
		interp_size = (DEM_SIZE - 1) / (dem_size_x - 1);
	
		/*
		 * Read in the entire DEM file into image_tmp.  We will select
		 * the chunks we need and process them into image_in later.
		 *
		 * Each record we read is a south-to-north slice of the DEM block.  Successive records move from
		 * west to east.  Thus, we read each record into a one-dimensional array, and then copy it
		 * into the desired two-dimensional storage area, simulaneously rotating the data so that north
		 * is at row zero and west is at column zero.
		 */
		dem_size_y = -1;
		for (i = 0; i < DEM_SIZE; i = i + interp_size)  {
			if (gz_flag == 0)  {
				if ((ret_val = read(dem_fdesc, buf, 8192)) != 8192)  {
					fprintf(stderr, "read from DEM file returns %d, expected 8192\n", ret_val);
					exit(0);
				}
			}
			else  {
				if ((ret_val = buf_read_z(dem_fdesc, buf, 8192)) != 8192)  {
					fprintf(stderr, "zread from DEM file returns %d, expected 8192\n", ret_val);
					exit(0);
				}
			}

			if (dem_size_y < 0)  {
				dem_size_y = strtol(&buf[12], (char **)0, 10);
				if (dem_size_y != DEM_SIZE)  {
					fprintf(stderr, "Number of rows in DEM file is %d, and should be 1201.\n", dem_size_y);
					exit(0);
				}

			}

			ptr = &buf[144];	/* Ignore header information on each block */

			for (j = DEM_SIZE - 1; j >=0; j--)  {
				*(image_tmp + j * DEM_SIZE + i) = strtol(ptr, &ptr, 10);

				/*
				 * If there are less than 1201 south-north columns, then interpolate to form
				 * a full 1201x1201 dataset.  That way the program only needs to handle one
				 * dataset size, and things are a lot easier.
				 */
				if ((interp_size > 1) && (i > 0))  {
					if (interp_size == 2)  {
						*(image_tmp + j * DEM_SIZE + i - 1) = round(0.5 * (double)(*(image_tmp + j * DEM_SIZE + i) +
										       *(image_tmp + j * DEM_SIZE + i - 2)));
					}
					else  {
						f = (double)(*(image_tmp + j * DEM_SIZE + i) -
							     *(image_tmp + j * DEM_SIZE + i - 3)) / 3.0;
						g = (double)*(image_tmp + j * DEM_SIZE + i - 3);
						*(image_tmp + j * DEM_SIZE + i - 2) = round(g + f);
						*(image_tmp + j * DEM_SIZE + i - 1) = round(g + f + f);
					}
				}
			}
		}
	
	
		/*
		 * Figure out the area of the image that will be covered by this set of DEM file data.
		 * Fill in that area with data from image_tmp.
		 *
		 * Because the relative sizes can take any ratio (in either the x or y direction)
		 * we simply choose the point from image_tmp that lies closest to the relative
		 * location in the covered area.  The exception to this is when the image is
		 * being subsampled, in which case we smooth the data to get representative data points.
		 *
		 * Note:  All 250K DEM files represent 1-degree by 1-degree blocks (as far as I know).
		 */
		latitude1 = max3(-91.0, latitude_dem, latitude_low);
		longitude1 = max3(-181.0, longitude_dem, longitude_low);
		latitude2 = min3(91.0, latitude_dem + 1.0, latitude_high);
		longitude2 = min3(181.0, longitude_dem + 1.0, longitude_high);
		tmp_width = round((double)(DEM_SIZE - 1) * (longitude2 - longitude1));
		tmp_height = round((double)(DEM_SIZE - 1) * (latitude2 - latitude1));
		tmp_x = round((double)(DEM_SIZE - 1) * (longitude1 - longitude_dem));
		tmp_y = DEM_SIZE - 1 - round((double)(DEM_SIZE - 1) * (latitude2 - latitude_dem));

		x_low = round((double)x * (longitude1 - longitude_low) / (longitude_high - longitude_low));
		x_high = round((double)(x + 1) * (longitude2 - longitude_low) / (longitude_high - longitude_low));
		y_low = y - round((double)y * (latitude2 - latitude_low) / (latitude_high - latitude_low));
		y_high = y + 1 - round((double)y * (latitude1 - latitude_low) / (latitude_high - latitude_low));

		if ((x_low < 0) || (x_high > (x + 1)) || (y_low < 0) || (y_high > (y + 1)))  {
			fprintf(stderr, "One of x_low=%d, x_high=%d, y_low=%d, y_high=%d out of range\n", x_low, x_high, y_low, y_high);
			exit(0);
		}


		/*
		 * Prepare a smoothing kernel in case we have more data than pixels to display it.
		 * The kernel is a square, a maximum of 2*SMOOTH_MAX+1 on a side.
		 *
		 * Here is one possible kernel, that I have tried:
		 *    If a kernel element is a distance of x from the
		 *    center, then its weight is 10*1.5^(-x/2)
		 *    Implemented by:
		 *       smooth[k + smooth_size][l + smooth_size] = round(10.0 * pow(1.0, - sqrt(k * k + l * l) / 2.0));
		 *
		 * For now, we just take the straight average over the kernel, since it seems to work reasonbly
		 * well.
		 */
		smooth_size = min3(SMOOTH_MAX,
				   (long)(600.0 * (longitude_high - longitude_low) / (double)x),
				   (long)(600.0 * (latitude_high - latitude_low) / (double)y));
		for (k = -smooth_size; k <= smooth_size; k++)  {
			for (l = -smooth_size; l <= smooth_size; l++)  {
				smooth[k + smooth_size][l + smooth_size] = 1;
			}
		}


		/*
		 * This is the loop that transfers the data for a single DEM into the image_in array.
		 * The image_in array will eventually hold the data from all DEM files given by the user.
		 */
		for (i = y_low; i < y_high; i++)  {
			for (j = x_low; j < x_high; j++)  {
				k = tmp_y + round((double)(tmp_height * (i - y_low)) / (double)(y_high - 1 - y_low));
				l = tmp_x + round((double)(tmp_width * (j - x_low)) / (double)(x_high - 1 - x_low));

				if ((l < 0) || (l > (DEM_SIZE - 1)) || (k < 0) || (k > (DEM_SIZE - 1)))  {
					fprintf(stderr, "One of l=%d, k=%d out of range, (i=%d, j=%d, tmp_y=%d, tmp_x=%d, tmp_height=%d, tmp_width=%d)\n",
						l, k, i, j, tmp_y, tmp_x, tmp_height, tmp_width);
					exit(0);
				}

				/*
				 * Note:  In the denominator of the left-hand-side of each comparison,
				 * we should have simply "x" or "y", but we add 1 to x or y
				 * to avoid having the comparison succeed due to quantization error.
				 */
				if ((((latitude_high - latitude_low) / (double)(x + 1)) > (1.0 / (double)(DEM_SIZE - 1))) ||
				    (((latitude_high - latitude_low) / (double)(y + 1)) > (1.0 / (double)(DEM_SIZE - 1))))  {
					/*
					 * We have an image that has less than (DEM_SIZE - 1)
					 * pixels per degree of longitude.
					 *
					 * Do some smoothing.
					 */
					sum = 0;
					sum_count = 0;
					for (m = -smooth_size; m <= smooth_size; m++)  {
						for (n = -smooth_size; n <= smooth_size; n++)  {
							if (((k + m) < 0) || ((k + m) >= DEM_SIZE) || ((l + n) < 0) || ((l + n) >= DEM_SIZE))  {
								continue;
							}

							sum += *(image_tmp + (k + m) * DEM_SIZE + l + n) * smooth[m + smooth_size][n + smooth_size];
							sum_count += smooth[m + smooth_size][n + smooth_size];
						}
					}
					*(image_in + i * (x + 1) + j) = round((double)sum / (double)sum_count);
				}
				else  {
					/*
					 * We have an image that is either one-to-one, or has
					 * more than (DEM_SIZE - 1) pixels per degree of longitude.
					 *
					 * Don't do any smoothing.  Simply pick the nearest
					 * point from image_tmp.
					 *
					 * If the x and y image size, given by the user, is
					 * not related by an integer factor to the number of elevation samples
					 * in the available data, then the image will contain some
					 * stripy anomalies because the rounding (above) to arrive
					 * at the k and l values will periodically give two k or
					 * l values in a row that have the same value.  Since
					 * the image color at a given point depends on changes in
					 * elevation around that point, having repeated elevation
					 * values can result in anomalous flat areas (with a neutral
					 * color) in an area of generally steep terrain (with generally
					 * bright or dark colors).
					 */
					*(image_in + i * (x + 1) + j) = *(image_tmp + k * DEM_SIZE + l);
				}
			}
		}
		if (gz_flag == 0)  {
			close(dem_fdesc);
		}
		else  {
			buf_close_z(dem_fdesc);
		}
	}
	free(image_tmp);
	

	/*
	 * Get memory for the actual output image.
	 *
	 * The index values for image_out, in the "y" and "x" directions, can each
	 * be -1 (when the latitude goes to latitude_high or the longitude goes to
	 * longitude_low, respectively).
	 * We allow this negative index because it makes it conceptually easier to
	 * translate latitudes and longitudes into x and y index values.
	 * We depend on the borders around the image to absorb these negative values
	 * so that we don't scribble memory outside the memory assigned to image_out.
	 *
	 *	(Say that we have an image that covers a 1x1 degree block, and
	 *	 1200x1200 pixels.  It makes sense that we assign either latitude_low
	 *	 or latitude_high to the zeroeth pixel in the x direction.  In this
	 *	 program, latitude_low/longitude_low is placed in the lower left
	 *	 corner of the image (at pixel (-1,1199)), and latitude increases going
	 *	 toward the top while longitude increases going toward the right.
	 *       (Remember that, in the western hemisphere, longitude is negative.)
	 *	 Thus, in theory, latitude_high/longitude_high is at the top right
	 *	 corner of the image (pixel (1199,-1)).  If you think about it,
	 *	 it makes sense that latitude_high and longitude_low are actually
	 *	 outside the image space since, as we reach a new integral value
	 *	 for these two values, we should begin cycling over at 0 again
	 *	 for a pixel index.  In other words, latitude_high and longitude_low
	 *	 would be part of the three images corresponding to the three 1x1 degree
	 *	 blocks that border the top left corner of the image we are working on.
	 *	 In order to take this into account, we set up the calculations so that
	 *	 longitude_low corresponds to an x pixel index of -1, and latitude_high
	 *	 corresponds to a y pixel index of -1.  Thus, as we generate the image,
	 *	 there is an extra strip of pixels along the top edge and the left edge.
	 *	 We come back later and change these to WHITE so that they merge back
	 *	 into the WHITE border around the edge of the actual map area.
	 *	 If you look closely at a map for a 1 degree by 1 degree block, you
	 *	 will note that the tick marks for the low longitude and high latitude
	 *	 are actually one pixel beyond the edge of the map area, at the locations
	 *	 that would be specified by horizontal and vertical indexes of -1.
	 *	 This may seem convoluted and silly, but it made sense
	 *	 to me at the time, and it does make it conceptually easier to
	 *	 verify that everything is working properly, and to join blocks together
	 *	 when a map spans more than one 1-degree DEM block.)
	 *
	 * Note that these negative index values aren't used for the DEM data, since
	 * the gradient calculations throw the extra data away.
	 */
	x_prime = x + LEFT_BORDER + RIGHT_BORDER;
	image_out = (unsigned char *)malloc((y + TOP_BORDER + BOTTOM_BORDER) * x_prime);
	if (image_out == (unsigned char *)0)  {
		fprintf(stderr, "malloc of image_out failed\n");
		exit(0);
	}


	/*
	 * Do the big calculation for processing the DEM data.
	 *
	 * This is where we transform elevation data into pixel colors
	 * in the output image.
	 *
	 * Note that the zeroeth row and zeroeth column of the elevation data
	 * (in image_in) are discarded during this process.  They consist of the
	 * data that would be plotted at the -1 horizontal and vertical index values
	 * in image_out.
	 *
	 * In order to do the gradient calculations, we need to have some
	 * reduction/expansion factors to figure out the approximate ground distance
	 * between image points.  We could normalize all of this out, but then it would be
	 * difficult to understand the geometrical implications when trying to modify the
	 * algorithms.
	 */
	reduction_y = (double)(DEM_SIZE - 1) / (double)y;
	reduction_xy = (double)DEM_SIZE * sqrt(2.0 / ((double)x * (double)x + (double)y * (double)y));

	for (i = 1; i <= y; i++)  {
		for (j = 1; j <= x; j++)  {
			/*
			 * Note:  spacing of sample points is nominally about 92.76 meters
			 * in the continental U.S.  Note the ridiculously precise value for a number that
			 * is inherently approximate.  Actually, this number is reasonably accurate
			 * for north-south spacing, except for the errors introduced by deviations of the
			 * earth from a perfect sphere.  But, for east-west spacing, the number is
			 * only correct at the equator and approaches zero as we approach the poles.
			 * Although the spacing obviously varies as we move away from the equator (and
			 * a 1x1 degree block starts to look more trapezoidal) it
			 * doesn't matter all that much to the gradient calculations
			 * since the number is only used to get a gross idea of the shape
			 * of the terrain so that we can assign a color to it.  At the extreme
			 * ends of the latitude range, we probably get a little distortion of the
			 * gradient shading, but it shouldn't be a really big deal.  If it becomes
			 * a big deal, it wouldn't be all that difficult to change the number a bit,
			 * based on the latitude.  At present, it doesn't seem worthwhile.
			 */

			/*
			 * We vary the shade of the DEM data to correspond to the slope of the
			 * terrain at each point.  The gradient calculations determine the slope
			 * in two directions and choose the largest of the two.
			 *
			 * The basic idea is to assume that the sun is shining from the northwest
			 * corner of the image.  Then, terrain with a negative gradient (toward the northwest
			 * or west) will be brightly colored, and terrain with a positive gradient will be dark, and
			 * level terrain will be somewhere in between.
			 */
			tangent1 = (((double)*(image_in + (i - 1) * (x + 1) + j - 1)) - ((double)*(image_in + i * (x + 1) + j))) / (92.76425 * reduction_xy);
			tangent2 = (((double)*(image_in + (i - 1) * (x + 1) + j)) - ((double)*(image_in + i * (x + 1) + j))) / (92.76425 * reduction_y);
			tangent3 = -10000000000.0;
			tangent = max3(tangent1, tangent2, tangent3);

			factor = get_factor(tangent);
/*			histogram[factor]++; */	/* Information for debugging. */


			/*
			 * Set the color based on the elevation and the factor
			 * retrieved from the gradient calculations.
			 * This is called a "factor" for historical reasons.
			 * At one time, I experimented with finding a multiplicative
			 * factor instead of the current additive modifier.
			 * It wasn't worth going through the code and changing the name.
			 * Besides, I might want to try a factor again someday.
			 *
			 * Elevations are in meters.  The intervals are chosen
			 * so that each color represents 1000 feet.
			 */
			if (*(image_in + i * (x + 1) + j) < 0)  {
				/* Elevations can theoretically be less than 0, but it's unusual, so report it. */
				/* fprintf(stderr, "An elevation was less than 0:  %d\n", *(image_in + i * (x + 1) + j)); */

				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = CYAN + factor;
			}
			else if (*(image_in + i * (x + 1) + j) == 0)  {
				/*
				 * Special case for sea level.  If things are totally flat,
				 * assume it's water.  Otherwise treat it like it's Death Valley.
				 *
				 * If this proves to be a problem, then get rid of the water assumption.
				 * The reason it is there is that the DLG files for coastal regions
				 * don't appear to treat oceans as bodies of water.  This was resulting
				 * in the ocean areas being set to GREEN.  Thus, I kludged in this special
				 * check; and it appears to work fine, in general.
				 *
				 * I later commented it out because, for example, sacramento-w.gz gets colored
				 * oddly, because there are areas below sea level within areas that meet the
				 * criterion for ocean.
				 */
/*
				if (tangent == 0.0)  {
					*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = B_BLUE;
				}
				else  {
					*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = CYAN + factor;
				}
*/
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = GREEN + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 305)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = GREEN + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 610)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = YELLOW + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 914)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = ORANGE + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 1219)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = RED + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 1524)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = MAGENTA + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 1829)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = CYAN + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 2134)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_GREEN + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 2438)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_YELLOW + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 2743)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_ORANGE + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 3048)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_RED + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 3353)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_MAGENTA + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 3658)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_BLUE + factor;
			}
			else if (*(image_in + i * (x + 1) + j) <= 3962)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = L_CYAN + factor;
			}
			else if (*(image_in + i * (x + 1) + j) < HIGHEST_ELEVATION)  {
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = WHITE + factor;
			}
			else  {
				/*
				 * Special case for creating WHITE areas by setting the
				 * DEM elevation data to HIGHEST_ELEVATION.
				 */
				*(image_out + (i - 1 + TOP_BORDER) * x_prime + j - 1 + LEFT_BORDER) = WHITE;
			}
		}
	}
	free(image_in);


	/*
	 * Process any DLG files.
	 * These files contain line and area information, for drawing
	 * things like streams, roads, boundaries, lakes, and such.
	 */
	file_index = 0;
	while (file_index < num_dlg)  {
		if (strcmp(argv[optind + file_index] + strlen(argv[optind + file_index]) - 3, ".gz") == 0)  {
			gz_flag = 1;
			if ((dlg_fdesc = buf_open_z(argv[optind + file_index], O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", argv[optind + file_index], errno);
				exit(0);
			}
		}
		else  {
			gz_flag = 0;
			if ((dlg_fdesc = buf_open(argv[optind + file_index], O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", argv[optind + file_index], errno);
				exit(0);
			}
		}

		fprintf(stderr, "Processing DLG file:  %s\n", argv[optind + file_index]);

		process_optional_dlg_file(dlg_fdesc, gz_flag);

		file_index++;

		if (gz_flag == 0)  {
			buf_close(dlg_fdesc);
		}
		else  {
			buf_close_z(dlg_fdesc);
		}
	}


	/* Select a font size, based on the image size. */
	if ((x >= 1000) && (y >= 1000))  {
		font_width = 6;
		font_height = 10;
		font = &font_6x10[0][0];
	}
	else  {
		font_width = 5;
		font_height = 8;
		font = &font_5x8[0][0];
	}


	/*
	 * Process any GNIS data.
	 * GNIS data consists of place names, with specific latitude/longitude
	 * coordinates, and other data.  We put a cursor at each given location
	 * and add the place name beside it.
	 */
	if (gnis_file != (char *)0)  {
		if (strcmp(gnis_file + strlen(gnis_file) - 3, ".gz") == 0)  {
			gz_flag = 1;
			if ((gnis_fdesc = buf_open_z(gnis_file, O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", gnis_file, errno);
				exit(0);
			}
		}
		else  {
			gz_flag = 0;
			if ((gnis_fdesc = buf_open(gnis_file, O_RDONLY)) < 0)  {
				fprintf(stderr, "Can't open %s for reading, errno = %d\n", gnis_file, errno);
				exit(0);
			}
		}

		fprintf(stderr, "Processing GNIS file:  %s\n", gnis_file);

		while ( 1 )  {
			/*
			 * Note that there is nothing special about the "170" in the following "if"
			 * statements.  I chose it because it is large enough to guarantee that
			 * we get enough GNIS data to include both the place name and the latitude/longitude.
			 * I didn't make it as large as a full-sized GNIS record (239 bytes with a newline, 240
			 * bytes with a newline and carriage return), because I thought it might be useful to
			 * be able to chop the unused data from the end of each line in a GNIS file and save
			 * some storage space.
			 */
			if (gz_flag == 0)  {
				if ((ret_val = get_a_line(gnis_fdesc, buf, BUF_SIZE)) <= 170)  {
					break;
				}
			}
			else  {
				if ((ret_val = get_a_line_z(gnis_fdesc, buf, BUF_SIZE)) <= 170)  {
					break;
				}
			}

			/*
			 * Note:  We assume latitude_low, longitude_low, latitude_high, and longitude_high
			 * were entered in decimal degrees.
			 * latitude and longitude from the GNIS files, however are in DDDMMSS format, and
			 * require special conversion functions.
			 */
			latitude = lat_conv(&buf[149]);
			longitude = lon_conv(&buf[156]);
	
			/* Ignore this entry if it is out of the map area. */
			if ((latitude < latitude_low) || (latitude > latitude_high))  {
				continue;
			}
			if ((longitude < longitude_low) || (longitude > longitude_high))  {
				continue;
			}
	
			/* draw a cursor at the specified point */
			xx = - 1 + round((longitude - longitude_low) * (double)x / (longitude_high - longitude_low));
			yy = y - 1 - round((latitude - latitude_low) * (double)y / (latitude_high - latitude_low));
	
			for (i = -3; i <= 3; i++)  {
				if (((xx + i) >= 0) && ((xx + i) <= (x - 1)))  {
					*(image_out + (yy + TOP_BORDER) * x_prime + xx + LEFT_BORDER + i) = dem_flag == 0 ? BLACK : WHITE;
				}
				if (((yy + i) >= 0) && ((yy + i) <= (y - 1)))  {
					*(image_out + (yy + TOP_BORDER + i) * x_prime + xx + LEFT_BORDER) = dem_flag == 0 ? BLACK : WHITE;
				}
			}
	
			/* Check for a feature name and find it's length */
			for (length = 79; length >= 0; length--)  {
				if (buf[length] != ' ')  {
					break;
				}
			}
			length++;
			buf[length] = '\0';
	
			/* If there was a feature name, then put it into the image near the cursor */
			if (length > 0)  {
				if ((xx + 5 + length * font_width) >= x)  {
					start_x = xx - 5 - length * font_width;
				}
				else  {
					start_x = xx + 5;
				}
				if ((yy + (font_height >> 1) - 1) >= y)  {
					start_y = y - font_height;
				}
				else if ((yy - (font_height >> 1)) < 0)  {
					start_y = 0;
				}
				else  {
					start_y = yy - (font_height >> 1);
				}
	
				add_text(buf, length, start_x + LEFT_BORDER,
					 start_y + TOP_BORDER, font, font_width, font_height, WHITE, -1);
			}
		}
		if (gz_flag == 0)  {
			buf_close(gnis_fdesc);
		}
		else  {
			buf_close_z(gnis_fdesc);
		}
	}


	/*
	 * Put a white border around the edges of the output image.
	 * Note that this will cover up the one-pixel slop over the left
	 * and top edges that is the result of the fact that we
	 * set the latitude and longitude to whole-number values, while
	 * the pixels don't quite cover that whole area.
	 * This was discussed at length in a previous comment.
	 *
	 * Note that the DEM file data don't slop over the edges because,
	 * when we process them, they are already in the form of an array of
	 * points, and we can cleanly discard the data we don't need.
	 * However, the DLG and GNIS data are in the form of
	 * latitude/longitude or UTM grid coordinates, and it is possible
	 * for array index values of -1 to crop up at the image edges.
	 * (In the case of GNIS data, we explicitly check for this and
	 * don't slop over.  For DLG data, we don't bother because it is
	 * cheaper in CPU time to just null out the border here.)
	 */
	for (i = 0; i < TOP_BORDER; i++)  {
		for (j = 0; j < (x + LEFT_BORDER + RIGHT_BORDER); j++)  {
			*(image_out + i * x_prime + j) = WHITE;
		}
	}
	for (i = y + TOP_BORDER; i < (y + TOP_BORDER + BOTTOM_BORDER); i++)  {
		for (j = 0; j < (x + LEFT_BORDER + RIGHT_BORDER); j++)  {
			*(image_out + i * x_prime + j) = WHITE;
		}
	}
	for (i = TOP_BORDER; i < (y + TOP_BORDER); i++)  {
		for (j = 0; j < LEFT_BORDER; j++)  {
			*(image_out + i * x_prime + j) = WHITE;
		}
	}
	for (i = TOP_BORDER; i < (y + TOP_BORDER); i++)  {
		for (j = x + LEFT_BORDER; j < (x + LEFT_BORDER + RIGHT_BORDER); j++)  {
			*(image_out + i * x_prime + j) = WHITE;
		}
	}


	/*
	 * Add a copyright notice to the image if, when the program was compiled, the
	 * makefile contained a non-null COPYRIGHT_NAME.
	 */
	if (COPYRIGHT_NAME[0] != '\0')  {
		time_val = time((time_t *)0);
		sprintf(buf, "Copyright (c) %4.4s  %s", ctime(&time_val) + 20, COPYRIGHT_NAME);
		length = strlen(buf);
		add_text(buf, length, x + LEFT_BORDER + RIGHT_BORDER - (length * font_width + 4),
			 y + TOP_BORDER + BOTTOM_BORDER - font_height - 4, font, font_width, font_height, BLACK, WHITE);
	}


	/*
	 * Put some latitude/longitude tick marks on the edges of the image.
	 *
	 * The purpose of the 0.049999999999 is to round the latitude/longitude up to
	 * the nearest tenth.  Since we put a tick mark every tenth of a degree,
	 * we need to find the first round tenth above latitude_low/longitude_low.
	 */
	i = (long)round((latitude_low + 0.049999999999) * 10.0);
	for (; i <= (latitude_high * 10.0); i++)  {
		k = TOP_BORDER - 1 + y - round((double)y * ((double)i * 0.1 - latitude_low) / (latitude_high - latitude_low));
		if (((i % 10) == 0) || ((i % 10) == 5) || ((i % 10) == -5))  {
			tick_width = 6;

			sprintf(buf, "%.2f %c", fabs((double)i / 10.0), i < 0 ? 'S' : 'N');
			length = strlen(buf);
			add_text(buf, length, x + LEFT_BORDER + 7, k - (font_height >> 1), font, font_width, font_height, BLACK, WHITE);
			add_text(buf, length, LEFT_BORDER - 8 - font_width * length, k - (font_height >> 1), font, font_width, font_height, BLACK, WHITE);
		}
		else  {
			tick_width = 4;
		}

		for (j = LEFT_BORDER - 1; j > (LEFT_BORDER - 1 - tick_width); j--)  {	/* Left side */
			*(image_out + k * x_prime + j) = BLACK;
		}
		for (j = x + LEFT_BORDER; j < (x + LEFT_BORDER + tick_width); j++)  {	/* Right side */
			*(image_out + k * x_prime + j) = BLACK;
		}
	}
	i = (long)round((longitude_low + 0.049999999999) * 10.0);
	for (; i <= (longitude_high * 10.0); i++)  {
		k = LEFT_BORDER - 1 + round((double)x * ((double)i * 0.1 - longitude_low) / (longitude_high - longitude_low));

		if (((i % 10) == 0) || ((i % 10) == 5) || ((i % 10) == -5))  {
			tick_width = 6;

			sprintf(buf, "%.2f %c", fabs((double)i / 10.0), i < 0 ? 'W' : 'E');
			length = strlen(buf);
			add_text(buf, length, k - ((length * font_width) >> 1), y + TOP_BORDER + 6, font, font_width, font_height, BLACK, WHITE);
			add_text(buf, length, k - ((length * font_width) >> 1), TOP_BORDER - 7 - font_height, font, font_width, font_height, BLACK, WHITE);
		}
		else  {
			tick_width = 4;
		}

		for (j = TOP_BORDER - 1; j > (TOP_BORDER - 1 - tick_width); j--)  {	/* Top */
			*(image_out + j * x_prime + k) = BLACK;
		}
		for (j = y + TOP_BORDER; j < (y + TOP_BORDER + tick_width); j++)  {	/* Bottom */
			*(image_out + j * x_prime + k) = BLACK;
		}
	}


	/* Add some information at the top of the image, as an image label (if there is room). */
	if (dem_name[0] != '\0')  {
		sprintf(buf, "%s --- ", dem_name);
	}
	else  {
		buf[0] = '\0';
	}
	sprintf(buf + strlen(buf), "%.2f %c, %.2f %c to %.2f %c, %.2f %c",
		fabs(latitude_low), latitude_low < 0 ? 'S' : 'N',
		fabs(longitude_low), longitude_low < 0 ? 'W' : 'E',
		fabs(latitude_high), latitude_high < 0 ? 'S' : 'N',
		fabs(longitude_high), longitude_high < 0 ? 'W' : 'E');
	length = strlen(buf);
	if ((length * font_width) <= x)  {
		add_text(buf, length, (x >> 1) + LEFT_BORDER - 1 - ((length * font_width) >> 1),
			 (TOP_BORDER >> 1) - 1 - (font_height >> 1), font, font_width,
			 font_height, BLACK, WHITE);
	}
	

	/* Add an elevation color chart at the bottom of the image, if there is room. */
	if ((num_dem > 0) && (x >= COLOR_CHART_WIDTH) && (BOTTOM_BORDER >= (30 + 3 * font_height)))  {
		for (i = 0; i < COLOR_CHART_WIDTH; i++)  {
			for (j = 0; j < 16; j++)  {
				/*
				 * To represent a given range of elevation, we draw a square of
				 * a given color.  We pick one of the 16 possible colors for each elevation.
				 * This is not perfect, but it at least gives the user some
				 * clue as to how to decode the image.  I tried filling in
				 * all 16 colors within each elevation square, but it didn't
				 * look all that good.
				 */
				*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) - ((16 + 4 + font_height * 2) >> 1) + j) * x_prime +
					LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + i) = (i & ~0xf) + 3;
			}
			if ((i & 0xf) == 0)  {
				/* Add a tick mark */
				*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 6 - font_height) * x_prime +
					LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0)) = BLACK;
				*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 7 - font_height) * x_prime +
					LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0)) = BLACK;
				*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 8 - font_height) * x_prime +
					LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0)) = BLACK;
				sprintf(buf, "%d", (i >> 4));
				length = strlen(buf);
				add_text(buf, length, LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0) - ((font_width * length) >> 1),
					TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 9 - font_height,
					font, font_width, font_height, BLACK, WHITE);
			}
		}

		/* Add a tick mark at the right end of the scale */
		*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 6 - font_height) * x_prime +
			LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0)) = BLACK;
		*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 7 - font_height) * x_prime +
			LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0)) = BLACK;
		*(image_out + (TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 8 - font_height) * x_prime +
			LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0)) = BLACK;

		/* Attempt to put in an "infinity" sign by jamming two 'o' characters together. */
		sprintf(buf, "o");
		length = strlen(buf);
		add_text(buf, length, LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0) - 1,
			TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 9 - font_height,
			font, font_width, font_height, BLACK, -2);
		add_text(buf, length, LEFT_BORDER + (x >> 1) - (COLOR_CHART_WIDTH >> 1) + (i & 0xf0) - ((font_width * length) >> 1) - 2,
			TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 9 - font_height,
			font, font_width, font_height, BLACK, -2);

		/* Add a line to describe the units. */
		sprintf(buf, "Thousands of feet.");
		length = strlen(buf);
		add_text(buf, length, (x >> 1) + LEFT_BORDER - 1 - ((length * font_width) >> 1),
			 TOP_BORDER + y + (BOTTOM_BORDER >> 1) + 9, font, font_width,
			 font_height, BLACK, WHITE);
	}


	/* Output the image data. */
	for (i = 0; i < (y + TOP_BORDER + BOTTOM_BORDER); i++)  {
		write(output_fdesc, image_out + i * x_prime, x + LEFT_BORDER + RIGHT_BORDER);
	}

	free(image_out);
	close(output_fdesc);


	/* For debugging. */
/*	for (i = 0; i < 256; i++)  {
/*		if (histogram[i] != 0)  {
/*			fprintf(stderr, "histogram[%3d] = %d\n", i, histogram[i]);
/*		}
/*	}
*/
}




/*
 * Convert elevation gradient information into an index that
 * can be used to select a color from the color table.
 * This routine was largely developed by trial and error.
 * There is no deep theory associated with the numeric values
 * contained herein.
 */
long
get_factor(double tangent)
{
	double angle, fraction;
	long i;

	/*
	 * A symmetrical table that works fairly well:
	 *
	 *	0.405,
	 *	0.445,
	 *	0.470,
	 *	0.485,
	 *	0.495,
	 *	0.497,
	 *	0.499,
	 *	0.500,
	 *	0.501,
	 *	0.503,
	 *	0.505,
	 *	0.515,
	 *	0.530,
	 *	0.555,
	 *	0.595,
	 *
	 * The table is duplicated in this comment so that we can
	 * play with the actual table without losing track of a set of
	 * values that work reasonably well.
	 */
	double table[15] =  {
		0.405,
		0.445,
		0.470,
		0.485,
		0.495,
		0.497,
		0.499,
		0.500,
		0.501,
		0.503,
		0.505,
		0.515,
		0.530,
		0.555,
		0.595,
	};

	/* One possible way to create the table automatically. */
/*	for (i = 0; i < 15; i++)  {
/*		table[i] = table[0] + (table[14] - table[0]) * pow((table[i] - table[0]) / (table[14] - table[0]), 0.9);
/*	}
*/

	angle = atan(tangent) + (M_PI/2.0);
	fraction = angle / (M_PI);
/*	angle_hist[round(fraction * 100000.0)]++; */	/* For debugging. */
/*	total++; */	/* For debugging. */

	if (fraction > 1.0)  {
		fprintf(stderr, "bad fraction in get_factor(%f):  %f\n", tangent, fraction);
	}

	if (fraction < table[0])  {
		return(0);
	}
	else if (fraction < table[1])  {
		return(1);
	}
	else if (fraction < table[2])  {
		return(2);
	}
	else if (fraction < table[3])  {
		return(3);
	}
	else if (fraction < table[4])  {
		return(4);
	}
	else if (fraction < table[5])  {
		return(5);
	}
	else if (fraction < table[6])  {
		return(6);
	}
	else if (fraction < table[7])  {
		return(7);
	}
	else if (fraction < table[8])  {
		return(8);
	}
	else if (fraction < table[9])  {
		return(9);
	}
	else if (fraction < table[10])  {
		return(10);
	}
	else if (fraction < table[11])  {
		return(11);
	}
	else if (fraction < table[12])  {
		return(12);
	}
	else if (fraction < table[13])  {
		return(13);
	}
	else if (fraction < table[14])  {
		return(14);
	}
	else  {
		return(15);
	}
}




/*
 * Write a text string into the image.
 */
void
add_text(char *text_string, long text_string_length, long top_left_x,
	 long top_left_y, unsigned char *font, long font_width, long font_height, long foreground, long background)
{
	long i, j, k;
	long bit;

	/*
	 * Cycle through the font table for each given character in the text string.
	 * Characters are represented as bit maps, with a 1 indicating part of the
	 * character, and a 0 indicating part of the background.
	 */
	for (i = 0; i < text_string_length; i++)  {
		for (j = 0; j < font_width; j++)  {
			for (k = 0; k < font_height; k++)  {
				bit = (*(font + k * 128 + *(text_string + i)) >> (font_width - 1 - j)) & 1;
				if (bit != 0)  {
					/* foreground */
					*(image_out  + (top_left_y + k) * x_prime + top_left_x + i * font_width + j) = foreground;
				}
				else  {
					/* background */
					if (background < 0)  {
						/*
						 * If the background color map index is -1, then
						 * we don't insert a specific background value, but rather
						 * reduce the existing background in intensity.
						 *
						 * If the background color map index is any other negative
						 * number, then we use a clear background.
						 */
						if (background == -1)  {
							*(image_out + (top_left_y + k) * x_prime + top_left_x + i * font_width + j) +=
								(16 - (*(image_out + (top_left_y + k) * x_prime + top_left_x + i * font_width + j) & 0xf)) >> 1;
						}
					}
					else  {
						*(image_out + (top_left_y + k) * x_prime + top_left_x + i * font_width + j) = background;
					}
				}
			}
		}
	}
}




/*
 * The code that processes DLG files is very spaghetti-like, since
 * it got squeezed and twisted and stretched while I figured out how
 * DLG files are put together.
 *
 * Because of this, and because I don't like to write functions that
 * take 35 arguments (the add_text() function, above, greatly offends
 * my sensibilities), there are a lot of global variables used by the
 * DLG code.  Most of them are accumulated here.
 */
#define	BOUNDARIES		90
#define	HYDROGRAPHY		50
#define	HYPSOGRAPHY		20
#define	PIPE_TRANS_LINES	190
#define	RAILROADS		180
#define	ROADS_AND_TRAILS	170
#define	PUBLIC_LAND_SURVEYS	300

struct attribute  {
	long major;
	long minor;
	struct attribute *attribute;
};
/*
 * The sizes of the nodes, areas, and lines arrays are their theoretical maximum values.
 * It would probably be cooler to malloc() these as we go, but coolness was not an
 * objective of this program.  It would still be cool to read the maximum values from
 * the DLG file headers and check them against the values below to verify that
 * the standards haven't changed and left this program behind.
 */
struct nodes  {
	long id;
	double x;
	double y;
} nodes[25960];
struct areas  {
	long id;
	double x;
	double y;
	long number_attrib;
	struct attribute *attribute;
} areas[25960];
struct lines  {
	long id;
	long start_node;
	long end_node;
	long left_area;
	long right_area;
	long number_coords;
	struct point *point;
	long number_attrib;
	struct attribute *attribute;
} lines[25938];

double lat_reference, long_reference;
double grid_x_reference, grid_y_reference;
double lat_se, long_se, lat_sw, long_sw, lat_ne, long_ne, lat_nw, long_nw;
double grid_x_se, grid_y_se, grid_x_sw, grid_y_sw, grid_x_ne, grid_y_ne, grid_x_nw, grid_y_nw;
double c1, c2, c3, c4, lat_offset, long_offset;
long dlg_x_low, dlg_y_low, dlg_x_high, dlg_y_high;


/*
 * Process the data from an optional-format 100K DLG file.
 * If you haven't read the DLG file guide and looked at a
 * DLG file, this code will probably be incomprehensible.
 */
void
process_optional_dlg_file(int fdesc, int gz_flag)
{
	long i, j, ret_val;
	long count;
	long color;
	char *end_ptr;
	char buf[144];
	char buf2[144];
	struct point **current_point;
	struct point *tmp_point;
	struct attribute **current_attrib;
	struct attribute *tmp_attrib;
	long attrib;
	long line_list;
	long num_nodes = 0;
	long num_areas = 0;
	long num_lines = 0;
	long data_type = 0;
	double latitude1, longitude1, latitude2, longitude2;
	ssize_t (* read_function)(int, void *, size_t);


	if (gz_flag == 0)  {
		read_function = buf_read;
	}
	else  {
		read_function = buf_read_z;
	}

	/*
	 * There is a lot of information in the file header.  We extract
	 * those items we care about and ignore the rest.
	 * We aren't interested in the first 10 records (for now), so ignore them.
	 */
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "1 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "2 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "3 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "4 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "5 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "6 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "7 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "8 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "9 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "10 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "11 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	else  {
		for (i = 0; i < 80; i++)  {
			/* The DLG files use 'D' for exponentiation.  strtod() expects 'E' or 'e'. */
			if (buf[i] == 'D') buf[i] = 'E';
		}
		i = 6;
		lat_sw = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		long_sw = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_x_sw = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_y_sw = strtod(&buf[i], &end_ptr);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "12 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	else  {
		for (i = 0; i < 80; i++)  {
			/* The DLG files use 'D' for exponentiation.  strtod() expects 'E' or 'e'. */
			if (buf[i] == 'D') buf[i] = 'E';
		}
		i = 6;
		lat_nw = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		long_nw = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_x_nw = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_y_nw = strtod(&buf[i], &end_ptr);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "13 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	else  {
		for (i = 0; i < 80; i++)  {
			/* The DLG files use 'D' for exponentiation.  strtod() expects 'E' or 'e'. */
			if (buf[i] == 'D') buf[i] = 'E';
		}
		i = 6;
		lat_ne = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		long_ne = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_x_ne = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_y_ne = strtod(&buf[i], &end_ptr);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "14 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	else  {
		for (i = 0; i < 80; i++)  {
			/* The DLG files use 'D' for exponentiation.  strtod() expects 'E' or 'e'. */
			if (buf[i] == 'D') buf[i] = 'E';
		}
		i = 6;
		lat_se = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		long_se = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_x_se = strtod(&buf[i], &end_ptr);
		i = i + end_ptr - &buf[i];

		grid_y_se = strtod(&buf[i], &end_ptr);
	}
	if ((ret_val = read_function(fdesc, buf, 80)) != 80)  {
		fprintf(stderr, "15 record DLG read returns %d\n", ret_val);
		exit(0);
	}
	else  {
		switch(buf[0])  {
		case 'B':	/* BOUNDARIES */
			color = WHITE + 12;
			data_type = BOUNDARIES;
			break;
		case 'H':
			if (buf[2] == 'D')  {
				/* HYDROGRAPHY */
				color = B_BLUE;
				data_type = HYDROGRAPHY;
				break;
			}
			else  {
				/* HYPSOGRAPHY */
				color = L_ORANGE;
				data_type = HYPSOGRAPHY;
				break;
			}
		case 'P':
			if (buf[1] == 'I')  {
				/* PIPE & TRANS LINES */
				color = BLACK;
				data_type = PIPE_TRANS_LINES;
				break;
			}
			else  {
				/* PUBLIC LAND SURVEYS */
				color = BLACK;
				data_type = PUBLIC_LAND_SURVEYS;
				break;
			}
		case 'R':
			if (buf[1] == 'A')  {
				/* RAILROADS */
				color = BLACK;
				data_type = RAILROADS;
				break;
			}
			else  {
				/* ROADS AND TRAILS */
				color = B_RED;
				data_type = ROADS_AND_TRAILS;
				break;
			}
		default:
			fprintf(stderr, "Unknown record type %20.20s\n", buf);
			exit(0);
			break;
		}
	}

	/*
	 * Do a quick check here to find out if the data is off the map boundaries.
	 * If so, then we can return now and save a lot of work.
	 */
	if ((lat_sw > latitude_high) ||
	    (long_sw > longitude_high) ||
	    (lat_ne < latitude_low) ||
	    (long_ne < longitude_low))  {
		return;
	}


	/*
	 * Within the Optional-format DLG file, locations are specified with pairs of
	 * Universal Transverse Mercator (x,y) coordinates.
	 *
	 * The header information at the top of the DLG file gives 4 reference
	 * points for the corners of the polygon represented by the DLG data.  Here is a
	 * typical set of them:
	 *
	 *	SW       45.750000 -112.000000         422218.03  5066539.80                    
	 *	NW       46.000000 -112.000000         422565.07  5094315.16                    
	 *	NE       46.000000 -111.750000         441923.83  5094103.38                    
	 *	SE       45.750000 -111.750000         441663.14  5066327.07                    
	 *
	 * Note that the latitude-longitude points form a square area in latitude/longitude
	 * space (if latitudes and longitudes on a pseudo-sphere can ever be thought of as
	 * forming a square).  The UTM (x,y) grid coordinates, however, form a quadrilateral
	 * in which no two sides have the same length.  Thus, if we are to convert the grid
	 * points in the DLG file into latitudes and longitudes, we need to develop a general
	 * transformation between these grid points and the desired latitudes and longitudes.
	 *
	 * We begin by using the latitude and longitude (and associated
	 * grid points) of the UTM origin (for this zone) as a reference,
	 * and subtracting the reference values from the
	 * coordinates for the four corners to arrive at:
	 *
	 *	lat_sw, long_sw,	grid_y_sw, grid_x_sw
	 *	lat_se, long_ne,	grid_y_se, grid_x_se
	 *	lat_nw, long_nw,	grid_y_nw, grid_x_nw
	 *	lat_ne, long_ne,	grid_y_ne, grid_x_ne
	 *
	 * As mentioned in the UTM tutorial above, the reference for latitude
	 * is the equator (zero degrees latitude), and the reference for
	 * longitude is the Central Meridian (in this example case, it is at
	 * 111 degrees west latitude for zone 12).  See the discussion of zones
	 * in an earlier comment.  The vertical (latitude) grid reference is zero
	 * at the equator, and the horizontal (longitude) grid reference is 500,000
	 * at the Central Meridian.
	 *
	 * Now we write the matrix equation:
	 *
	 * ---           ---   ---    ---          ---         ---        ---  ---
	 * |               |   |        |          |             |        |      |
	 * | c1         c2 |   | grid_y |          | lat_offset  |        | lat  |
	 * |               |   |        |    +     |             |   =    |      |
	 * |               |   |        |          |             |        |      |
	 * | c3         c4 |   | grid_x |          | long_offset |        | long |
	 * |               |   |        |          |             |        |      |
	 * ---           ---   ---    ---          ---         ---        ---  ---
	 *
	 * From which we can calculate a desired latitude and longitude by plugging
	 * given grid_x and grid_y values into the equations (after subtracting the
	 * reference values):
	 *
	 *	lat  = c1 * grid_y + c2 * grid_x + lat_offset
	 *	long = c3 * grid_y + c4 * grid_x + long_offset
	 *
	 * First, however, we must find values for c1, c2, c3, c4, lat_offset, and long_offset.
	 * Plugging in the data from three of the corners, we get six equations for
	 * six unknowns.  Our handy-dandy symbolic math program spits out equations
	 * for the unknowns.  These equations appear in the following code.
	 *
	 * Note:  This linear mapping from grid values into latitudes and longitudes is
	 * a linear approximation to the what actually occurs on the surface of the (approximately)
	 * spherical Earth.  Corrections could be applied to make the map more accurate, but we
	 * leave that for another day.
	 *
	 * Note further:  Enterprising readers will want to verify that the equations
	 * convert the unused pair of grid values into the correct latitude/longitude
	 * values for the fourth corner (the corner that wasn't used to find the 6 constants).
	 *
	 * For reference:
	 * 4.0076594e7 meters is the equatorial circumference of the earth.
	 * 3.9942e7 meters is the polar circumference of the earth.
	 * 4.0024e7 meters is the average circumference of the earth.
	 */
	lat_reference = 0.0;
	long_reference = -186.0 + (double)(zone * 6) + 3.0;
	grid_y_reference = 0.0;
	grid_x_reference = 500000.0;

	lat_nw = lat_nw - lat_reference;
	long_nw = long_nw - long_reference;
	lat_ne = lat_ne - lat_reference;
	long_ne = long_ne - long_reference;
	lat_sw = lat_sw - lat_reference;
	long_sw = long_sw - long_reference;
	lat_se = lat_se - lat_reference;
	long_se = long_se - long_reference;

	grid_y_nw = grid_y_nw - grid_y_reference;
	grid_x_nw = grid_x_nw - grid_x_reference;
	grid_y_ne = grid_y_ne - grid_y_reference;
	grid_x_ne = grid_x_ne - grid_x_reference;
	grid_y_sw = grid_y_sw - grid_y_reference;
	grid_x_sw = grid_x_sw - grid_x_reference;
	grid_y_se = grid_y_se - grid_y_reference;	/* Never used */
	grid_x_se = grid_x_se - grid_x_reference;	/* Never used */

	c1 = (lat_nw*grid_x_ne - lat_sw*grid_x_ne - lat_ne*grid_x_nw +
	      lat_sw*grid_x_nw + lat_ne*grid_x_sw - lat_nw*grid_x_sw) /
	     (-(grid_x_nw*grid_y_ne) + grid_x_sw*grid_y_ne + grid_x_ne*grid_y_nw - grid_x_sw*grid_y_nw -
	      grid_x_ne*grid_y_sw + grid_x_nw*grid_y_sw);
	c2 = (lat_nw*grid_y_ne - lat_sw*grid_y_ne - lat_ne*grid_y_nw +
	      lat_sw*grid_y_nw + lat_ne*grid_y_sw - lat_nw*grid_y_sw) /
	     (grid_x_nw*grid_y_ne - grid_x_sw*grid_y_ne - grid_x_ne*grid_y_nw + grid_x_sw*grid_y_nw +
	      grid_x_ne*grid_y_sw - grid_x_nw*grid_y_sw);
	c3 = (long_nw*grid_x_ne - long_sw*grid_x_ne - long_ne*grid_x_nw +
	      long_sw*grid_x_nw + long_ne*grid_x_sw - long_nw*grid_x_sw) /
	     (-(grid_x_nw*grid_y_ne) + grid_x_sw*grid_y_ne + grid_x_ne*grid_y_nw -
	      grid_x_sw*grid_y_nw - grid_x_ne*grid_y_sw + grid_x_nw*grid_y_sw);
	c4 = (long_nw*grid_y_ne - long_sw*grid_y_ne - long_ne*grid_y_nw +
	      long_sw*grid_y_nw + long_ne*grid_y_sw - long_nw*grid_y_sw) /
	     (grid_x_nw*grid_y_ne - grid_x_sw*grid_y_ne - grid_x_ne*grid_y_nw +
	      grid_x_sw*grid_y_nw + grid_x_ne*grid_y_sw - grid_x_nw*grid_y_sw);
	lat_offset = (lat_sw*grid_x_nw*grid_y_ne - lat_nw*grid_x_sw*grid_y_ne - lat_sw*grid_x_ne*grid_y_nw +
		      lat_ne*grid_x_sw*grid_y_nw + lat_nw*grid_x_ne*grid_y_sw - lat_ne*grid_x_nw*grid_y_sw) /
		     (grid_x_nw*grid_y_ne - grid_x_sw*grid_y_ne - grid_x_ne*grid_y_nw +
		      grid_x_sw*grid_y_nw + grid_x_ne*grid_y_sw - grid_x_nw*grid_y_sw);
	long_offset = (long_sw*grid_x_nw*grid_y_ne - long_nw*grid_x_sw*grid_y_ne - long_sw*grid_x_ne*grid_y_nw +
		       long_ne*grid_x_sw*grid_y_nw + long_nw*grid_x_ne*grid_y_sw - long_ne*grid_x_nw*grid_y_sw) /
		      (grid_x_nw*grid_y_ne - grid_x_sw*grid_y_ne - grid_x_ne*grid_y_nw +
		       grid_x_sw*grid_y_nw + grid_x_ne*grid_y_sw - grid_x_nw*grid_y_sw);


	/*
	 * Following the header information, there is a sequence of data records for
	 * Nodes, Areas, and Lines.
	 * Parse these data records and put the data into the appropriate arrays.
	 * At present, we make absolutely no use of the Node information.
	 */
	while ((ret_val = read_function(fdesc, buf, 80)) == 80)  {
		switch(buf[0])  {
		case 'N':
			i = 1;
			nodes[num_nodes].id = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			nodes[num_nodes].x = strtod(&buf[i], &end_ptr);
			i = i + end_ptr - &buf[i];

			nodes[num_nodes].y = strtod(&buf[i], &end_ptr);

			i = 36;
			line_list = strtol(&buf[i], &end_ptr, 10);

			i = 48;
			attrib = strtol(&buf[i], &end_ptr, 10);

			if (line_list != 0)  {
				while(line_list > 0)  {
					if ((ret_val = read_function(fdesc, buf2, 80)) != 80)  {
						fprintf(stderr, "Line_list read 1 returns %d\n", ret_val);
						fprintf(stderr, "%80.80s\n", buf);
						exit(0);
					}

					line_list = line_list - 12;
				}
			}

			if (attrib != 0)  {
				while (attrib > 0)  {
					if ((ret_val = read_function(fdesc, buf2, 80)) != 80)  {
						fprintf(stderr, "Attribute read 1 returns %d\n", ret_val);
						fprintf(stderr, "%80.80s\n", buf);
						exit(0);
					}

					attrib = attrib - 6;
				}
			}

			num_nodes++;
			break;

		case 'A':
			i = 1;
			areas[num_areas].id = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			areas[num_areas].x = strtod(&buf[i], &end_ptr);
			i = i + end_ptr - &buf[i];

			areas[num_areas].y = strtod(&buf[i], &end_ptr);

			i = 36;
			line_list = strtol(&buf[i], &end_ptr, 10);

			i = 48;
			attrib = strtol(&buf[i], &end_ptr, 10);
			areas[num_areas].number_attrib = attrib;

			if (line_list != 0)  {
				while (line_list > 0)  {
					if ((ret_val = read_function(fdesc, buf2, 80)) != 80)  {
						fprintf(stderr, "Line_list read 2 returns %d\n", ret_val);
						fprintf(stderr, "%80.80s\n", buf);
						exit(0);
					}

					line_list = line_list - 12;
				}
			}

			if (attrib != 0)  {
				while (attrib > 0)  {
					if ((ret_val = read_function(fdesc, buf2, 80)) != 80)  {
						fprintf(stderr, "Attribute read 2 returns %d\n", ret_val);
						fprintf(stderr, "%80.80s\n", buf);
						exit(0);
					}

					current_attrib = &areas[num_areas].attribute;

					if (attrib > 6)  {
						i = 6;
						attrib = attrib - 6;
					}
					else  {
						i = attrib;
						attrib = 0;
					}

					end_ptr = buf2;

					while (i > 0)  {
						*current_attrib = (struct attribute *)malloc(sizeof(struct attribute));
						if (*current_attrib == (struct attribute *)0)  {
							fprintf(stderr, "malloc failed\n");
							exit(0);
						}

						(*current_attrib)->major = strtol(end_ptr, &end_ptr, 10);
						(*current_attrib)->minor = strtol(end_ptr, &end_ptr, 10);

						current_attrib = &((*current_attrib)->attribute);
						i--;
					}
					*current_attrib = (struct attribute *)0;
				}
			}

			num_areas++;
			break;

		case 'L':
			i = 1;
			lines[num_lines].id = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			lines[num_lines].start_node = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			lines[num_lines].end_node = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			lines[num_lines].left_area = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			lines[num_lines].right_area = strtol(&buf[i], &end_ptr, 10);

			i = 42;
			lines[num_lines].number_coords = strtol(&buf[i], &end_ptr, 10);
			i = i + end_ptr - &buf[i];

			attrib = strtol(&buf[i], &end_ptr, 10);
			lines[num_lines].number_attrib = attrib;

			current_point = &lines[num_lines].point;
			count = lines[num_lines].number_coords;
			while (count != 0)  {
				if ((ret_val = read_function(fdesc, buf2, 80)) != 80)  {
					fprintf(stderr, "Coordinate read returns %d\n", ret_val);
					fprintf(stderr, "%80.80s\n", buf);
					exit(0);
				}

				i = 0;
				while (i < 80)  {
					while ((i < 80) && (buf2[i] == ' '))  {
						i++;
					}
					if (i >= 80)  {
						continue;
					}

					*current_point = (struct point *)malloc(sizeof(struct point));
					if (*current_point == (struct point *)0)  {
						fprintf(stderr, "malloc failed\n");
						exit(0);
					}

					(*current_point)->x = (long)strtod(&buf2[i], &end_ptr);
					i = i + end_ptr - &buf2[i];
					(*current_point)->y = (long)strtod(&buf2[i], &end_ptr);
					i = i + end_ptr - &buf2[i];

					current_point = &((*current_point)->point);
					count--;
				}
			}
			*current_point = (struct point *)0;

			if (attrib != 0)  {
				while (attrib > 0)  {
					if ((ret_val = read_function(fdesc, buf2, 80)) != 80)  {
						fprintf(stderr, "Attribute read 3 returns %d\n", ret_val);
						fprintf(stderr, "%80.80s\n", buf);
						exit(0);
					}

					current_attrib = &lines[num_lines].attribute;

					if (attrib > 6)  {
						i = 6;
						attrib = attrib - 6;
					}
					else  {
						i = attrib;
						attrib = 0;
					}

					end_ptr = buf2;
					while (i > 0)  {
						*current_attrib = (struct attribute *)malloc(sizeof(struct attribute));
						if (*current_attrib == (struct attribute *)0)  {
							fprintf(stderr, "malloc failed\n");
							exit(0);
						}

						(*current_attrib)->major = strtol(end_ptr, &end_ptr, 10);
						(*current_attrib)->minor = strtol(end_ptr, &end_ptr, 10);

						current_attrib = &((*current_attrib)->attribute);
						i--;
					}
					*current_attrib = (struct attribute *)0;
				}
			}

			num_lines++;
			break;

		default:
			fprintf(stderr, "Unknown record type: %c  (hexadecimal: %x)\n", buf[0], buf[0]);
			fprintf(stderr, "%80.80s\n", buf);
			exit(0);
			break;
		}
	}


	/*
	 * All of the useful data is parsed.
	 * Now do something with it.
	 *
	 * First find the x and y image coordinates that border this DLG chunk.
	 *
	 * Then draw the lines that we have found, but don't go outside
	 * the x-y border.
	 *
	 * Then fill in all of the areas for which we have
	 * appropriate attribute codes stored, but don't go outside
	 * the x-y border.
	 */
	latitude1 = lat_se + lat_reference;
	longitude1 = long_se + long_reference;
	latitude2 = lat_nw + lat_reference;
	longitude2 = long_nw + long_reference;
	dlg_x_low = -1 + round((longitude2 - longitude_low) * (double)x / (longitude_high - longitude_low));
	dlg_y_low = y - 1 - round((latitude2 - latitude_low) * (double)y / (latitude_high - latitude_low));
	dlg_x_high = -1 + round((longitude1 - longitude_low) * (double)x / (longitude_high - longitude_low));
	dlg_y_high = y - 1 - round((latitude1 - latitude_low) * (double)y / (latitude_high - latitude_low));
	if (dlg_x_low < -1)  {
		dlg_x_low = -1;
	}
	if (dlg_y_low < -1)  {
		dlg_y_low = -1;
	}
	if (dlg_x_high >= x)  {
		dlg_x_high = x - 1;
	}
	if (dlg_y_high >= y)  {
		dlg_y_high = y - 1;
	}

	/*
	 * Cycle through all of the line data and draw all of the appropriate lines
	 * onto the image (overlaying any previous data).
	 */
	for (i = 0; i < num_lines; i++)  {
		/*
		 * In the DLG-3 format, the first area element listed
		 * represents the universe outside of the map area.
		 * Thus, lines that have area 1 as a boundary should be
		 * "neatlines" that bound the map area.
		 * Since these clutter up a map, we normally discard them.
		 * (If you want to keep them, then change the #define of OMIT_NEATLINES
		 * so that it is zero, rather than non-zero.)
		 *
		 * Here are relevant quotes from the DLG-3 guide:
		 *
		 *	expressed by network data is that of connectivity.  The network case
		 *	differs from the area case in that, irrespective of the number of closed
		 *	areas forming the graph, only two areas are encoded:  (1) the area out-
		 *	side the graph, termed the outside area; and (2) the area within the
		 *	graph, termed the background area.  All lines except the graph boundary,
		 *	or neatline, are considered to be contained within the background area.
		 *
		 *	map border.  There is one outside area for each DLG-3. It is always the
		 *	first area encountered (its ID is 1) and has the attribute code 000 0000.
		 */

		/*
		 * If the user provided a file full of attributes, then
		 * use them to control whether or not the lines are drawn.
		 * If not, then just go ahead and draw everything.
		 *
		 * Note:  If a major or minor attribute code in the attribute
		 *        file (supplied by the user) is less than
		 *        zero, it is treated as a wild card and matches
		 *        anything.
		 */
		if ((num_A_attrib > 0) || (num_L_attrib > 0))  {
			if ((OMIT_NEATLINES == 0) || ((lines[i].left_area != 1) && (lines[i].right_area != 1)))  {
				current_attrib = &lines[i].attribute;
				if (*current_attrib != (struct attribute *)0)  {
					while (*current_attrib != (struct attribute *)0)  {
						for (j = 0; j < num_L_attrib; j++)  {
							if (((attributes_L[j].major < 0) ||
							     (attributes_L[j].major == ((*current_attrib)->major))) &&
							    ((attributes_L[j].minor < 0) ||
							     (attributes_L[j].minor == ((*current_attrib)->minor))))  {
								draw_lines(lines[i].point, color);
								goto FIN1;
							}
						}
						current_attrib = &((*current_attrib)->attribute);
					}
				}
				else  {
					/*
					 * If the feature had no attribute codes, then check if
					 * it is covered by a wild card in the attributes file.
					 */
					for (j = 0; j < num_L_attrib; j++)  {
						if (((attributes_L[j].major < 0) ||
						     (attributes_L[j].major == data_type)) &&
						    (attributes_L[j].minor < 0))  {
							draw_lines(lines[i].point, color);
							goto FIN1;
						}
					}
				}
			}

			/*
			 * For those (hopefully rare) occasions in which something
			 * goes wrong, we provide the capability for a user to
			 * specifically request a single line from a DLG file so that
			 * the cause of the problem can be isolated.
			 * The user specifies a specific line by providing a major
			 * attribute number of 10000, and a minor attribute number
			 * equal to the desired line ID number.  Since no
			 * valid attribute (as far as I know) is ever as large as
			 * 10,000, such user-specified attribute pairs will not
			 * affect the search for legitimate attributes above (since
			 * they can't possibly match anything).  If we reach this point,
			 * then we failed to draw a line due to the legitimate-attribute
			 * checks above; so we give it one more try here, based on
			 * user-requested ID numbers.
			 *
			 * Note:  If you are using this feature, then it doesn't make
			 *        a lot of sense to process more than one DLG file,
			 *        since the ID number you give (as the minor attribute)
			 *        will be matched in every DLG file that has a
			 *        Line with that ID.  If you are trying to isolate
			 *        one (or a few) Line(s), then you probably want to
			 *        be certain which file is the source of the data.
			 */
			for (j = 0; j < num_L_attrib; j++)  {
				if ((attributes_L[j].major == 10000) &&
				     (attributes_L[j].minor == lines[i].id))  {
					draw_lines(lines[i].point, color);
					goto FIN1;
				}
			}
		}
		else  {
			if ((OMIT_NEATLINES == 0) || ((lines[i].left_area != 1) && (lines[i].right_area != 1)))  {
				draw_lines(lines[i].point, color);
			}
		}
FIN1:
	}

	/*
	 * Now we fill in each interesting area on the map with the
	 * same color that bounds the area.  (For example,
	 * lakes (attribute code 050 0421) might be filled in.)
	 * However, sometimes areas might be filled in improperly.
	 * The code assumes that the reference point for an area falls
	 * within the polygon of lines that define that area.
	 * According to the DLG guide, this isn't guaranteed
	 * to always be the case, but the assumption has nonetheless
	 * worked well in practice.
	 *
	 * Area attributes are processed a bit differently than the
	 * attributes for lines:  no areas are filled in automatically.
	 * If the user did not specify any Area attributes in the attribute
	 * file, then no areas are filled in.  This is because the area-fill
	 * algorithm can occasionally run amok, and therefore the appropriate
	 * default is to not give it a chance.  For extensive details on the
	 * area-fill algorithm, see the comments at the top of fill_area().
	 */
	if (num_A_attrib > 0)  {
		for (i = 0; i < num_areas; i++)  {
			if (areas[i].number_attrib <= 0)  {
				continue;
			}

			current_attrib = &areas[i].attribute;
			while (*current_attrib != (struct attribute *)0)  {
				for (j = 0; j < num_A_attrib; j++)  {
					if (((attributes_A[j].major < 0) ||
					     (attributes_A[j].major == ((*current_attrib)->major))) &&
					    ((attributes_A[j].minor < 0) ||
					     (attributes_A[j].minor == ((*current_attrib)->minor))))  {
						fill_area(areas[i].x, areas[i].y, color);
						goto FIN2;
					}
				}
				current_attrib = &((*current_attrib)->attribute);
			}

			/*
			 * As with the Line attributes, we provide an interface
			 * for the user to select specific areas, via their IDs.
			 */
			for (j = 0; j < num_A_attrib; j++)  {
				if ((attributes_A[j].major == 10000) &&
				     (attributes_A[j].minor == areas[i].id))  {
					fill_area(areas[i].x, areas[i].y, color);
					goto FIN2;
				}
			}
FIN2:
		}
	}


	/* Free up all of the malloc() memory */
	for (i = 0; i < num_lines; i++)  {
		if (lines[i].number_coords > 0)  {
			current_point = &lines[i].point;

			while (*current_point != (struct point *)0)  {
				tmp_point = (*current_point)->point;
				free(*current_point);
				*current_point = tmp_point;
			}
		}
		if (lines[i].number_attrib > 0)  {
			current_attrib = &lines[i].attribute;

			while (*current_attrib != (struct attribute *)0)  {
				tmp_attrib = (*current_attrib)->attribute;
				free(*current_attrib);
				*current_attrib = tmp_attrib;
			}
		}
	}
	for (i = 0; i < num_areas; i++)  {
		if (areas[i].number_attrib > 0)  {
			current_attrib = &areas[i].attribute;

			while (*current_attrib != (struct attribute *)0)  {
				tmp_attrib = (*current_attrib)->attribute;
				free(*current_attrib);
				*current_attrib = tmp_attrib;
			}
		}
	}
}



/*
 * Draw a series of line segments, as defined by a linked list of
 * points from an optional-format DLG file.
 *
 * This routine is recursive, not because it has to be, but because
 * it was slightly simpler that way.  Since it doesn't recurse very
 * far (on average), it isn't a performance or memory problem.
 *
 * It is a nasty routine to understand, because it has a generalized
 * interpolation algorithm to capture line segments that go beyond the
 * image boundaries.
 */
void
draw_lines(struct point *cur_point, long color)
{
	double latitude1, longitude1;
	double latitude2, longitude2;
	long xx1, yy1;
	long xx2, yy2;
	double fxx, fyy;
	double delta_x, delta_y;
	long steps;
	long i;
	double m_lat, m_long, b_lat, b_long;
	double p_lat1, p_long1, p_lat2, p_long2;
	double d_lat, d_long;
	long pointflags = 0;
	long bothflag = 0;

	/*
	 * We recurse to the end of the linked list, and then draw line
	 * segments as we pop back up the recursion stack.
	 */
	if (cur_point->point != (struct point *)0)  {
		draw_lines(cur_point->point, color);

		/*
		 * Draw a segment between this point and the next one down the linked list.
		 *
		 * Begin by figuring out the latitude and longitude of the endpoints.
		 */
		latitude1 = lat_reference + lat_offset +
			    c1 * (cur_point->y - grid_y_reference) +
			    c2 * (cur_point->x - grid_x_reference);
		longitude1 = long_reference + long_offset +
			     c3 * (cur_point->y - grid_y_reference) +
			     c4 * (cur_point->x - grid_x_reference);
		latitude2 = lat_reference + lat_offset +
			    c1 * ((cur_point->point)->y - grid_y_reference) +
			    c2 * ((cur_point->point)->x - grid_x_reference);
		longitude2 = long_reference + long_offset +
			     c3 * ((cur_point->point)->y - grid_y_reference) +
			     c4 * ((cur_point->point)->x - grid_x_reference);


		/*
		 * Find out whether only one endpoint, or both of them, fall
		 * outside the map area.
		 */
		if ((latitude1 < latitude_low) || (latitude1 > latitude_high) ||
		    (longitude1 < longitude_low) || (longitude1 > longitude_high))  {
			bothflag++;
		}
		if ((latitude2 < latitude_low) || (latitude2 > latitude_high) ||
		    (longitude2 < longitude_low) || (longitude2 > longitude_high))  {
			bothflag++;
		}


		/*
		 * If at least one endpoint of a line segment is outside of the area
		 * covered by the map image, then interpolate the segment.
		 *
		 * This isn't just to catch errors in a DLG file.  Since the user
		 * can specify arbitrary latitude/longitude boundaries for the
		 * map image, either or both endpoints of a segment can easily
		 * be outside of the map boundaries.
		 */
		if (bothflag > 0)  {
			/*
			 * Construct two equations for the line passing through the two
			 * endpoints.  These equations can be solved for four potential
			 * intercepts with the edge of the map area, only zero or two of
			 * which should be actual intercepts.  (In theory, there can
			 * be a single intercept at a corner, but this code should find
			 * it twice.)
			 *
			 * We construct the two lines using the classic Y = m * X + b formula,
			 * where, in one case, we let Y be the latitude and X be the longitude,
			 * and in the other case they switch roles.
			 */
			m_lat = (latitude2 - latitude1) / (longitude2 - longitude1);
			b_lat = latitude1 - m_lat * longitude1;
			m_long = 1.0 / m_lat;
			b_long = longitude1 - m_long * latitude1;

			/*
			 * We need the distance (in the Manhattan (city-block) metric) between
			 * the two endpoints.
			 * It will be used to determine whether one of the intercepts with
			 * the map edges falls between the two given endpoints.
			 */
			d_lat = fabs(latitude1 - latitude2);
			d_long = fabs(longitude1 - longitude2);

			/*
			 * Solve the two equations for the four possible intercepts, and check
			 * that they are truly intercepts.
			 * Set a flag to remember which points turned out to be intercepts.
			 */
			p_lat1 = m_lat * longitude_low + b_lat;
			if ((p_lat1 >= latitude_low) && (p_lat1 <= latitude_high))  {
				if ((fabs(longitude_low - longitude1) <= d_long) && (fabs(longitude_low - longitude2) <= d_long))  {
					pointflags |= 1;
				}
			}
			p_lat2 = m_lat * longitude_high + b_lat;
			if ((p_lat2 >= latitude_low) && (p_lat2 <= latitude_high))  {
				if ((fabs(longitude_high - longitude1) <= d_long) && (fabs(longitude_high - longitude2) <= d_long))  {
					pointflags |= 2;
				}
			}
			p_long1 = m_long * latitude_low + b_long;
			if ((p_long1 >= longitude_low) && (p_long1 <= longitude_high))  {
				if ((fabs(latitude_low - latitude1) <= d_lat) && (fabs(latitude_low - latitude2) <= d_lat))  {
					pointflags |= 4;
				}
			}
			p_long2 = m_long * latitude_high + b_long;
			if ((p_long2 >= longitude_low) && (p_long2 <= longitude_high))  {
				if ((fabs(latitude_high - latitude1) <= d_lat) && (fabs(latitude_high - latitude2) <= d_lat))  {
					pointflags |= 8;
				}
			}

			/*
			 * If both endpoints fall outside the map area, and there aren't exactly two
			 * intercepts, then there should be none.  (In theory, when a segment
			 * just touches a corner of the map area, then there is only one intercept,
			 * but the above code will find the same intercept twice.)
			 */
			if ((bothflag == 2) && (pointflags != 3) && (pointflags != 5) && (pointflags != 6) &&
			    (pointflags != 9) && (pointflags != 10) && (pointflags != 12))  {
				if (pointflags != 0)  {
			    		fprintf(stderr, " should have had exactly two intercepts:  0x%x  (%f %f) (%f %f)\n",
						pointflags, latitude1, longitude1, latitude2, longitude2);
				}
				return;
			}

			/* If the first endpoint is out of range, then replace it with an intercept. */
			if ((latitude1 < latitude_low) || (latitude1 > latitude_high) ||
			    (longitude1 < longitude_low) || (longitude1 > longitude_high))  {
				if (pointflags & 1)  {
					latitude1 = p_lat1;
					longitude1 = longitude_low;
					pointflags &= ~1;
					goto DONE1;
				}
				if (pointflags & 2)  {
					latitude1 = p_lat2;
					longitude1 = longitude_high;
					pointflags &= ~2;
					goto DONE1;
				}
				if (pointflags & 4)  {
					latitude1 = latitude_low;
					longitude1 = p_long1;
					pointflags &= ~4;
					goto DONE1;
				}
				if (pointflags & 8)  {
					latitude1 = latitude_high;
					longitude1 = p_long2;
					pointflags &= ~8;
					goto DONE1;
				}
			}
DONE1:

			/* If the second endpoint is out of range, then replace it with an intercept. */
			if ((latitude2 < latitude_low) || (latitude2 > latitude_high) ||
			    (longitude2 < longitude_low) || (longitude2 > longitude_high))  {
				if (pointflags & 1)  {
					latitude2 = p_lat1;
					longitude2 = longitude_low;
					goto DONE2;
				}
				if (pointflags & 2)  {
					latitude2 = p_lat2;
					longitude2 = longitude_high;
					goto DONE2;
				}
				if (pointflags & 4)  {
					latitude2 = latitude_low;
					longitude2 = p_long1;
					goto DONE2;
				}
				if (pointflags & 8)  {
					latitude2 = latitude_high;
					longitude2 = p_long2;
					goto DONE2;
				}
			}
DONE2:
		}



		/*
		 * Convert the latitude/longitude pairs into pixel locations within the image.
		 *
		 * Note:  because there may be small errors in longitude1, latitude1, longitude2,
		 * and latitude2, the values of xx1, yy1, xx2, or yy2 may occasionally be off by
		 * one pixel.
		 * This appears to be acceptable in the middle of the image, since one pixel
		 * doesn't amount to much linear distance in the image.  At the edges, one might
		 * worry that the discrepancy would cause us to go over the image edges.
		 * However, the interpolation code above should successfully eliminate this
		 * potential problem.
		 *
		 * As noted above, it is okay for the array index values to go to -1, since that
		 * is the appropriate value for longitude_low or latitude_high.
		 */
		xx1 = -1 + round((longitude1 - longitude_low) * (double)x / (longitude_high - longitude_low));
		yy1 = y - 1 - round((latitude1 - latitude_low) * (double)y / (latitude_high - latitude_low));
		xx2 = -1 + round((longitude2 - longitude_low) * (double)x / (longitude_high - longitude_low));
		yy2 = y - 1 - round((latitude2 - latitude_low) * (double)y / (latitude_high - latitude_low));
		if ((xx1 < -1) || (yy1 < -1) || (xx1 >= x) || (yy1 >= y))  {
			fprintf(stderr, "In drawlines(), a coordinate exceeds the image boundaries, %f %f   %f %f\n", xx1, yy1, xx2, yy2);
			exit(0);
		}


		/*
		 * Now all that remains is to draw the line segment.
		 * We begin by deciding whether x or y is the fastest-changing
		 * coordinate.
		 */
		delta_x = xx2 - xx1;
		delta_y = yy2 - yy1;

		if (fabs(delta_x) < fabs(delta_y))  {
			steps = (long)fabs(delta_y) - 1;

			if (delta_y > 0.0)  {
				delta_x = delta_x / delta_y;
				delta_y = 1.0;
			}
			else if (delta_y < 0.0)  {
				delta_x = -delta_x / delta_y;
				delta_y = -1.0;
			}
			else  {
				delta_x = 1.0;
			}
		}
		else  {
			steps = (long)fabs(delta_x) - 1;

			if (delta_x > 0.0)  {
				delta_y = delta_y / delta_x;
				delta_x = 1.0;
			}
			else if (delta_x < 0.0)  {
				delta_y = -delta_y / delta_x;
				delta_x = -1.0;
			}
			else  {
				delta_y = 1.0;
			}
		}

		/* Put dots at the two endpoints. */
		*(image_out + (yy1 + TOP_BORDER) * x_prime + xx1 + LEFT_BORDER) = color;
		*(image_out + (yy2 + TOP_BORDER) * x_prime + xx2 + LEFT_BORDER) = color;

		/* Fill in pixels between the two endpoints. */
		fxx = xx1;
		fyy = yy1;
		for (i = 0; i < steps; i++)  {
			fxx = fxx + delta_x;
			fyy = fyy + delta_y;
			*(image_out + (round(fyy) + TOP_BORDER) * x_prime + round(fxx) + LEFT_BORDER) = color;
		}
	}
}



/*
 * Fill in an area bounded by a polygon of the given color, beginning at the
 * given representative point.  (The polygon was previously created by the
 * line-drawing algorithm.)  The algorithm does this by filling in a given
 * point and then recursively calling itself to fill in the four nearest neighbors
 * (to the left, right, top, and bottom).
 *
 * Two functions handle this:  fill_area() sets things up, and then
 * fill_small_area() recursively does the work.  An enterprising reader might
 * want to convert the recursion into something less likely to consume all
 * computing resources on the planet.  However, these routines generally
 * work well unless somehow the representative point falls outside of a bounded
 * polygon (which appears not to be a common event, based on my limited testing).
 * If this happens, then, as the routine attempts to fill large swaths of the image,
 * the recursion chomps up all available stack memory and the program goes kaboom.
 * Less resources would be gobbled if, instead of using recursion, we simply built
 * a stack datatype, and pushed and popped coordinates onto/from it.  No program is
 * so perfect that it can't be improved.
 *
 * One other problem with the approach taken here is that, if a lake has a narrow
 * neck, the line segments at the sides of the neck may touch.  If that is the case,
 * then only one side of the lake will be filled in (the side containing the
 * representative point) because the neck forms a solid boundary.
 *
 * Yet another problem is that the representative point may be off the map boundaries
 * if, say, a lake is at the edge of the map and the whole lake doesn't show up on
 * the map.  In such a case, the lake won't get filled in because the representative
 * point is rejected by the sanity-checking code.
 *
 * This algorithm is very crude at this point.  We assume that the given
 * coordinates actually do fall within the bounded area that they represent,
 * something that the DLG guide says is normal for these points, but not guaranteed.
 * It would appear that a general solution not relying on this assumption would be
 * difficult to produce.  For a convex bounding polygon, one can determine if the
 * representative point is within the bounding polygon by following the line segments
 * around the boundaries of the area and checking that the point is always on the same side
 * of the line segment (relative to the direction of motion).  However, this wouldn't
 * do us a whole lot of good.  First, the polygons are not, in general, convex.
 * Second, unless we change the area fill algorithm in some fundamental way,
 * knowing a single point (that is guaranteed to be within the boundaries of the area)
 * still won't guarantee that the area gets filled properly (see the discussion of a
 * lake with a neck, above).  Third, knowing that a point is within the boundaries of
 * the area is not adequate to guarantee that it is within the boundaries drawn on
 * the image.  The lines drawn around the boundaries are "jagged", because we try
 * to draw slanted lines using pixels that
 * can only be placed on a square grid.  (This problem is often called "aliasing,"
 * which is a reference to Nyquist Sampling Theory; but that is a subject far
 * beyond the scope of this long, rambling comment block.)  It is possible for
 * the representative point to land on a pixel that falls outside the drawn
 * boundaries, because it just happens to fall at a place where a slanted line
 * segment "jags."  This problem can be exacerbated when the image is stretched
 * (for example, when a map area that is 2 degrees of longitude by 1 degree of
 * latitude is plotted on a 600 by 600 pixel grid, thus stretching the latitude
 * direction by a factor of 2).
 *
 * We also assume that the area is totally bounded on the right, left, top, and
 * bottom by points of the given color (or the map-area edges).  The line-drawing
 * algorithm, above, should ensure this, as long as the line segments given in the
 * DLG file don't leave gaps (which they normally don't appear to do).
 *
 * There may be some cool, sexy way to write an area-fill algorithm that would
 * be completely general and would run fast.  However, without giving it a massive
 * amount of thought, the only truly general algorithms I have come up with are very
 * slow, involving a detailed check of each candidate point to verify that it is indeed
 * withing the given area.  As an example, here is a very clunky algorithm that is
 * "guaranteed" to work without running amok:
 *
 * Determine which collection(s) of line segments is associated with the given area.
 *     (Multiple multi-segment linear features can bound an area, including the neatlines
 *      that bound the entire area encompassed by the DLG file.)
 * Follow the line segments around the bounding polygon and break the polygon into
 *     multiple polygons, each of which is convex.  This can be done by examining the
 *     angles between successive line segments.
 * For each convex sub-polygon:
 *     Find the largest and smallest longitude and latitude associated with all of the
 *         segments in the sub-polygon.
 *     Sweep through all points within the rectangle determined by the longitude/latitude
 *         bounding box and check each point to determine whether it is within the area
 *         in question.  This can be done by following the line segments around the polygon
 *         and checking that the point is always on the same side of each segment.  (The
 *         sign of the line segment identifier(s) determines which side the point is
 *         supposed to be on.  See the DLG documentation for details.)
 *
 * Although there is a lot of handwaving in the above description, it should be obvious
 * that this algorithm would be incredibly slow.  One could obviously come up with some
 * ways to speed it up, since it is designed for simplicity of description rather than
 * efficiency of operation, but it is not immediately obvious how to make it really fast.
 * Nor is it immediately obvious (at least to me) how to come up with a different algorithm
 * that would be both robust and fast.  Also, the current version appears to work pretty
 * well, with occasional inevitable glitches.  Thus, for the time being, we are stuck with
 * the code that follows.
 */
void
fill_small_area(long x1, long y1, long color)
{
	/*
	 * Check that we have not wandered outside of the area
	 * covered by the data from this DLG file.
	 */
	if ((x1 < dlg_x_low) || (x1 > dlg_x_high) || (y1 < dlg_y_low) || (y1 > dlg_y_high))  {
		return;
	}

	/*
	 * Fill in the given pixel, and recusively fill in the pixels to the
	 * left, right, top, and bottom.
	 */
	*(image_out + (y1 + TOP_BORDER) * x_prime + x1 + LEFT_BORDER) = color;

	if (*(image_out + (y1 - 1 + TOP_BORDER) * x_prime + x1 + LEFT_BORDER) != color)  {
		fill_small_area(x1, y1 - 1, color);
	}
	if (*(image_out + (y1 + 1 + TOP_BORDER) * x_prime + x1 + LEFT_BORDER) != color)  {
		fill_small_area(x1, y1 + 1, color);
	}
	if (*(image_out + (y1 + TOP_BORDER) * x_prime + x1 - 1 + LEFT_BORDER) != color)  {
		fill_small_area(x1 - 1, y1, color);
	}
	if (*(image_out + (y1 + TOP_BORDER) * x_prime + x1 + 1 + LEFT_BORDER) != color)  {
		fill_small_area(x1 + 1, y1, color);
	}
}
void
fill_area(double px1, double py1, long color)
{
	double latitude1, longitude1;
	long xx1, yy1;

	/* Find the latitude and longitude of the representative point and convert them into index values. */
	latitude1 = lat_reference + lat_offset + c1 * (py1 - grid_y_reference) + c2 * (px1 - grid_x_reference);
	longitude1 = long_reference + long_offset + c3 * (py1 - grid_y_reference) + c4 * (px1 - grid_x_reference);
	
	xx1 = -1 + round((longitude1 - longitude_low) * (double)x / (longitude_high - longitude_low));
	yy1 = y - 1 - round((latitude1 - latitude_low) * (double)y / (latitude_high - latitude_low));
	if ((xx1 < -1) || (xx1 >= x) || (yy1 < -1) || (yy1 >= y))  {
/*		fprintf(stderr, "fill_area() was given a starting point outside the map area:  (%d %d) (%f %f)\n", xx1, yy1, latitude1, longitude1); */
		return;
	}

	if ((xx1 < dlg_x_low) || (xx1 > dlg_x_high) || (yy1 < dlg_y_low) || (yy1 > dlg_y_high))  {
		fprintf(stderr, "fill_area() was passed a bad starting point:  (%d %d) (%f %f)\n\tlimits are: %d %d   %d %d\n",
			xx1, yy1, latitude1, longitude1, dlg_x_low, dlg_x_high, dlg_y_low, dlg_y_high);
		return;
	}

/*
 * Some debugging code to figure out where the representative point
 * for each area falls on the image.
 */
/*
{
static h = 0;
long la, lo;
double long_prime = fabs(longitude1) - 0.5;
la = latitude1;
lo = long_prime;
la = la * 10000 + ((int)((latitude1 - la) * 60.0)) * 100 + (int)((latitude1 - la - ((int)((latitude1 - la) * 60.0)) / 60.0) * 3600.0 + 0.5);
lo = lo * 10000 + ((int)((long_prime - lo) * 60.0)) * 100 + (int)((long_prime - lo - ((int)((long_prime - lo) * 60.0)) / 60.0) * 3600.0 + 0.5);

fprintf(stderr, "lat=%f long=%f     %d %d\n", la, lo, xx1, yy1);
fprintf(stdout, "Area %2.2d                                                                             island   Blaine                                             30005%6.6dN%7.7dW                     %f %f                         \n", h, la, lo, px1, py1);
h++;

*(image_out + (yy1 + TOP_BORDER) * x_prime + xx1 + LEFT_BORDER) = B_GREEN;
return;
}
*/

	/*
	 * Some small areas are so small that the lines around their borders have
	 * already filled them in.  If the representative point is already set to
	 * the target color, then we assume we are in such an area.  In such cases,
	 * we immediately return, because otherwise (if we happen to be sitting
	 * right on the boundary) we will begin filling in the area outside the
	 * boundary and potentially fill large swaths of the image.  The risk of
	 * simply returning (rather than doing a more thorough investigation of
	 * what is going on) is that the boundary lines may not have actually
	 * filled the area in, but rather
	 * that the representative point just happens to fall very near
	 * (or on) the boundary.  There is not much we can do about this potential
	 * problem, unless we re-write the whole area-filling algorithm
	 * (not necessarily a bad idea).  However, in practice, things seem
	 * to generally work out okay for the data sets I have tried.
	 */
	if (*(image_out + (yy1 + TOP_BORDER) * x_prime + xx1 + LEFT_BORDER) == color)  {
		return;
	}

	/* Recursively call fill_small_area() to do most of the work. */
	fill_small_area(xx1, yy1, color);
}