/*--------------------------------------------------------------------
 *	$Id: grdimage.c,v 1.54 2006/04/13 08:34:41 pwessel Exp $
 *
 *	Copyright (c) 1991-2006 by P. Wessel and W. H. F. Smith
 *	See COPYING file for copying and redistribution conditions.
 *
 *	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; version 2 of the License.
 *
 *	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.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * grdimage will read a grdfile and image the area using the PostScript
 * image command. If non-linear scaling is chosen, we first resample the data
 * onto the new grid before calling the image command.  The output image
 * will be 1-, 8-, or 24-bit depending on colors used.  As an option, grdimage
 * can read three grdfiles with red, green, blue components in the 0-255
 * range and use those colors directly (no cpt required).
 *
 * Author:	Paul Wessel
 * Date:	17-SEP-2001
 * Ver:		4
 */

#include "gmt.h"
#include "pslib.h"

struct GRDIMAGE_CTRL {
	struct C {	/* -C<cptfile> */
		BOOLEAN active;
		char *file;
	} C;
	struct E {	/* -Ei|<dpi> */
		BOOLEAN active;
		BOOLEAN device_dpi;
		int dpi;
	} E;
	struct G {	/* -G[f|b]<rgb> */
		BOOLEAN active;
		int f_rgb[3];
		int b_rgb[3];
	} G;
	struct I {	/* -I<intensfile> */
		BOOLEAN active;
		char *file;
	} I;
	struct M {	/* -M */
		BOOLEAN active;
	} M;
	struct Q {	/* -Q */
		BOOLEAN active;
	} Q;
	struct S {	/* -S<radius> */
		BOOLEAN active;
		double radius;
	} S;
	struct T {	/* -T[s][o[<pen>] */
		BOOLEAN active;
		BOOLEAN skip;
		BOOLEAN outline;
		struct GMT_PEN pen;
	} T;
};

int main (int argc, char **argv)
{
	BOOLEAN error = FALSE, do_rgb = FALSE, need_to_project, normal_x, normal_y;

	unsigned char *bitimage_8 = NULL, *bitimage_24 = NULL;

	char *grdfile[3], *c_method[2] = {"colorimage", "colortiles",};

	int i, j, k, kk, rgb[3], nm, nm2 = 0, byte, nx_f, ny_f, grid_type = 0;
	int nx, ny, nx_proj = 0, ny_proj = 0, node, n_grids = 0, index;
	int PS_colormask_off = 0, PS_colormask, PS_interpolate;

	float *tmp1[3], *tmp2 = NULL, *map[3], *intensity = NULL;
	
	double  dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0;
	double west, east, south, north, data_west, data_east, data_south, data_north;

	struct GRD_HEADER g_head[3], r_head[3], i_head, j_head;
	struct GRDIMAGE_CTRL *Ctrl;

	void GMT_limits (struct GRD_HEADER *r, struct GRD_HEADER *g);
	void GMT_set_proj_limits (struct GRD_HEADER *r, struct GRD_HEADER *g);
	void *New_Grdimage_Ctrl (), Free_Grdimage_Ctrl (struct GRDIMAGE_CTRL *C);
	
	argc = GMT_begin (argc, argv);

	Ctrl = (struct GRDIMAGE_CTRL *)New_Grdimage_Ctrl ();	/* Allocate and initialize a new control structure */
	
	grdfile[0] = grdfile[1] = grdfile[2] = CNULL;
	west = east = south = north = 0.0;

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* Common parameters */

				case 'B':
				case 'J':
				case 'K':
				case 'O':
				case 'P':
				case 'R':
				case 'U':
				case 'V':
				case 'X':
				case 'x':
				case 'Y':
				case 'y':
				case 'c':
				case 'f':
				case '\0':
					error += GMT_parse_common_options (argv[i], &west, &east, &south, &north);
					break;

				/* Supplemental parameters */

				case 'C':
					Ctrl->C.active = TRUE;
					Ctrl->C.file = strdup (&argv[i][2]);
					break;
				case 'E':
					Ctrl->E.active = TRUE;
					if (argv[i][2] == 'i')	/* Interpolate image to device resolution */
						Ctrl->E.device_dpi = TRUE;
					else
						Ctrl->E.dpi = atoi (&argv[i][2]);
					break;
				case 'G':	/* 1-bit fore or background color for transparent masks */
					Ctrl->G.active = TRUE;
					switch (argv[i][2]) {
						case 'F':
						case 'f':
							if (GMT_getrgb (&argv[i][3], Ctrl->G.f_rgb)) {
								GMT_rgb_syntax ('G');
								error++;
							}
							else
								Ctrl->G.b_rgb[0] = -1;
							break;
						case 'B':
						case 'b':
							if (GMT_getrgb (&argv[i][3], Ctrl->G.b_rgb)) {
								GMT_rgb_syntax ('G');
								error++;
							}
							else
								Ctrl->G.f_rgb[0] = -1;
							break;
						default:	/* Same as -Gf */
							if (GMT_getrgb (&argv[i][2], Ctrl->G.f_rgb)) {
								GMT_rgb_syntax ('G');
								error++;
							}
							else
								Ctrl->G.b_rgb[0] = -1;
							break;
					}
					break;
				case 'I':
					Ctrl->I.active = TRUE;
					Ctrl->I.file = strdup (&argv[i][2]);
					break;
				case 'M':
					Ctrl->M.active = TRUE;
					break;
				case '0':	/* Backwards compatibility */
					gmtdefs.color_image = 0;
					break;
				case '1':	/* Backwards compatibility */
					gmtdefs.color_image = 1;
					break;
				case 'Q':
					Ctrl->Q.active = TRUE;
					break;
				case 'S':
					Ctrl->S.active = TRUE;
					Ctrl->S.radius = atof (&argv[i][2]);
					break;
				case 'T':
					Ctrl->T.active = TRUE;
					k = 2;
					if (argv[i][2] == 's') Ctrl->T.skip = TRUE, k = 3;
					if (argv[i][k] == 'o') {	/* Want tile outline also */
						Ctrl->T.outline = TRUE;
						k++;
						if (argv[i][k] && GMT_getpen (&argv[i][k], &Ctrl->T.pen)) {
							GMT_pen_syntax ('T');
							error++;
						}
					}
					break;

				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else {
			if (n_grids == 3) {
				fprintf (stderr, "%s: ERROR, give 1 or 3 grid files\n", GMT_program);
				exit (EXIT_FAILURE);
			}
			grdfile[n_grids++] = argv[i];
		}
	}

	if (argc == 1 || GMT_give_synopsis_and_exit) {
		fprintf (stderr,"grdimage %s - Plot grdfiles in 2-D\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdimage <grd_z|grd_r grd_g grd_b> %s [%s] [-C<cpt_file>] [-Ei|<dpi>] [-G[f|b]<rgb>]\n", GMT_J_OPT, GMT_B_OPT);
		fprintf (stderr, "\t[-I<intensity_file>] [-K] [-M] [-O] [-P] [-Q] [%s] [-S<radius>] [-T[s][o[<pen>]]]\n", GMT_Rgeo_OPT);
		fprintf (stderr, "\t[%s] [-V] [%s] [%s] [%s]\n\n", GMT_U_OPT, GMT_X_OPT, GMT_Y_OPT, GMT_c_OPT);

		if (GMT_give_synopsis_and_exit) exit (EXIT_FAILURE);

		fprintf (stderr,"\t<grd_z> is data set to be plotted.  Its z-values are in user units and will be\n");
		fprintf (stderr,"\t  converted to rgb colors via the cpt file.  Alternatively, give three separate\n");
		fprintf (stderr,"\t  grid files that contain the red, green, and blue components in the 0-255 range.\n");
		GMT_explain_option ('j');
		fprintf (stderr,"\n\tOPTIONS:\n");
		GMT_explain_option ('b');
		fprintf (stderr,"\t-C color palette file to convert z to rgb\n");
		fprintf(stderr, "\t-E sets dpi for the projected grid which must be constructed\n");
		fprintf(stderr, "\t   if -Jx or -Jm is not selected [Default gives same size as input grid]\n");
		fprintf(stderr, "\t   Give i to do the interpolation in PostScript at device resolution.\n");
		fprintf (stderr,"\t-G sets transparency color for images that\n");
		fprintf (stderr,"\t   otherwise would result in 1-bit images\n");
		fprintf (stderr,"\t-I use illumination.  Append name of intensity grd file\n");
		GMT_explain_option ('K');
		fprintf (stderr, "\t-M force monochrome image\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf(stderr, "\t-Q use PS Level 3 colormasking for nodes with z = NaN.\n");
		GMT_explain_option ('R');
		fprintf(stderr, "\t-S sets the search radius in projected units [Default avoids aliasing]\n");
		fprintf(stderr, "\t-T will image the data without interpolation by painting polygonal tiles\n");
		fprintf(stderr, "\t   Append s to skip tiles for nodes with z = NaN [Default paints all tiles]\n");
		fprintf(stderr, "\t   Append o[<pen>] to draw tile outline [Default uses no outline]\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		GMT_explain_option ('X');
		GMT_explain_option ('c');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	do_rgb = (n_grids == 3);

	for (i = 0; i < n_grids; i++) {
		if (!grdfile[i]) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
			error++;
		}
	}
	if (!Ctrl->C.file && !do_rgb) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify color palette table\n", GMT_program);
		error++;
	}
	if (Ctrl->I.active && !Ctrl->I.file) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify intensity file\n", GMT_program);
		error++;
	}
	if (Ctrl->E.active && !Ctrl->E.device_dpi && Ctrl->E.dpi <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  dpi must be positive\n", GMT_program);
		error++;
	}
	if (Ctrl->G.f_rgb[0] < 0 && Ctrl->G.b_rgb[0] < 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Only one of fore/back-ground can be transparent for 1-bit images\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */

	/* Get color palette file */

	if (!do_rgb) GMT_read_cpt (Ctrl->C.file);
	PS_interpolate = (Ctrl->E.device_dpi) ? -1 : +1;
	PS_colormask   = (Ctrl->Q.active) ? -1 : +1;;


	if (gmtdefs.verbose) fprintf (stderr, "%s: Allocates memory and read data file\n", GMT_program);

	for (i = 0; i < n_grids; i++) {
		GMT_grd_init (&g_head[i], argc, argv, FALSE);
		GMT_grd_init (&r_head[i], argc, argv, FALSE);

		if (GMT_read_grd_info (grdfile[i], &g_head[i])) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grdfile[i]);
			exit (EXIT_FAILURE);
		}
	}
	if (do_rgb) {	/* Must ensure all three grids are coregistered */
		if (!(g_head[0].x_min == g_head[1].x_min && g_head[0].x_min == g_head[2].x_min)) error++;
		if (!(g_head[0].y_min == g_head[1].y_min && g_head[0].y_min == g_head[2].y_min)) error++;
		if (!(g_head[0].x_inc == g_head[1].x_inc && g_head[0].x_inc == g_head[2].x_inc)) error++;
		if (!(g_head[0].nx == g_head[1].nx && g_head[0].nx == g_head[2].nx)) error++;
		if (!(g_head[0].ny == g_head[1].ny && g_head[0].ny == g_head[2].ny)) error++;
		if (!(g_head[0].node_offset == g_head[1].node_offset && g_head[0].node_offset == g_head[2].node_offset)) error++;
		if (error) {
			fprintf (stderr, "%s: The r, g, and b grids are not congruent\n", GMT_program);
			exit (EXIT_FAILURE);
		}
	}

	/* Determine what wesn to pass to map_setup */

	if (!project_info.region_supplied) {
		west = g_head[0].x_min;
		east = g_head[0].x_max;
		south = g_head[0].y_min;
		north = g_head[0].y_max;
	}
	
	GMT_map_setup (west, east, south, north);

	/* Determine the wesn to be used to read the grdfile; or exit if file is outside -R */

	if (GMT_grd_setregion (&g_head[0], &data_west, &data_east, &data_south, &data_north)) {	/* No grid to plot; just do empty map and exit */
		GMT_plotinit (argc, argv);
		GMT_map_basemap ();
		GMT_plotend ();
		GMT_end (argc, argv);
	}

	nx_f = g_head[0].nx;
	ny_f = g_head[0].ny;

	/* Read data */

	nx = GMT_get_n (data_west, data_east, g_head[0].x_inc, g_head[0].node_offset);
	ny = GMT_get_n (data_south, data_north, g_head[0].y_inc, g_head[0].node_offset);
	nm = nx * ny;

	for (i = 0; i < n_grids; i++) {
		tmp1[i] = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
		if (GMT_read_grd (grdfile[i], &g_head[i], tmp1[i], data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grdfile[i]);
			exit (EXIT_FAILURE);
		}
	}

	/* If given, get intensity file or compute intensities */

	if (Ctrl->I.active) {	/* Illumination wanted */

		if (gmtdefs.verbose) fprintf (stderr, "%s: Allocates memory and read intensity file\n", GMT_program);

		GMT_grd_init (&i_head, argc, argv, FALSE);
		GMT_grd_init (&j_head, argc, argv, FALSE);

		if (GMT_read_grd_info (Ctrl->I.file, &i_head)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, Ctrl->I.file);
			exit (EXIT_FAILURE);
		}

		if (i_head.nx != nx_f || i_head.ny != ny_f) {
			fprintf (stderr, "%s: Intensity file has improper dimensions!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		tmp2 = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);

		if (GMT_read_grd (Ctrl->I.file, &i_head, tmp2, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, Ctrl->I.file);
			exit (EXIT_FAILURE);
		}
	}

	for (i = 0; i < n_grids; i++) GMT_set_proj_limits (&r_head[i], &g_head[i]);
	need_to_project = (project_info.projection != GMT_LINEAR || project_info.xyz_projection[0] == GMT_LOG10 || project_info.xyz_projection[0] == GMT_POW || project_info.xyz_projection[1] == GMT_LOG10 || project_info.xyz_projection[1] == GMT_POW);

	if (!Ctrl->T.active && (Ctrl->E.dpi > 0 || need_to_project)) {	/* Need to resample the grd file */

		if (gmtdefs.verbose) fprintf (stderr, "%s: project grdfiles\n", GMT_program);

		if (Ctrl->E.dpi == 0) {	/* Use input # of nodes as # of projected nodes */
			nx_proj = g_head[0].nx;
			ny_proj = g_head[0].ny;
		}
		for (i = 0; i < n_grids; i++) {
			grid_type = (Ctrl->E.dpi > 0) ? 1 : g_head[i].node_offset;	/* Force pixel if dpi is set */
			GMT_grdproject_init (&r_head[i], 0.0, 0.0, nx_proj, ny_proj, Ctrl->E.dpi, grid_type);
			GMT_init_search_radius (&Ctrl->S.radius, &r_head[i], &g_head[i], FALSE);
			nm2 = r_head[i].nx * r_head[i].ny;
			map[i] = (float *) GMT_memory (VNULL, (size_t)nm2, sizeof (float), "grdproject");
			GMT_grd_forward (tmp1[i], &g_head[i], map[i], &r_head[i], Ctrl->S.radius);
			GMT_free ((void *)tmp1[i]);
		}
		if (Ctrl->I.active) {
			j_head.x_min = r_head[0].x_min;	j_head.x_max = r_head[0].x_max;
			j_head.y_min = r_head[0].y_min;	j_head.y_max = r_head[0].y_max;

			if (Ctrl->E.dpi == 0) {	/* Use input # of nodes as # of projected nodes */
				nx_proj = i_head.nx;
				ny_proj = i_head.ny;
			}
			GMT_grdproject_init (&j_head, 0.0, 0.0, nx_proj, ny_proj, Ctrl->E.dpi, grid_type);
			intensity = (float *) GMT_memory (VNULL, (size_t)nm2, sizeof (float), "grdproject");
			GMT_grd_forward (tmp2, &i_head, intensity, &j_head, Ctrl->S.radius);
			GMT_free ((void *)tmp2);
		}
		nm = nm2;
		/* project_info.xyz_pos[0] = project_info.xyz_pos[1] = 1; */	/* Since we have applied the projection we do no longer have 'negative' axes */
	}
	else {	/* Simply copy g_head[0] info to r_head[0] */
		for (i = 0; i < n_grids; i++) {
			map[i] = tmp1[i];
			r_head[i].nx = g_head[i].nx;		r_head[i].ny = g_head[i].ny;
			r_head[i].x_inc = g_head[i].x_inc;	r_head[i].y_inc = g_head[i].y_inc;
		}
		if (Ctrl->I.active) {
			j_head.nx = i_head.nx;		j_head.ny = i_head.ny;
			j_head.x_inc = i_head.x_inc;	j_head.y_inc = i_head.y_inc;
			intensity = tmp2;
		}
		grid_type = g_head[0].node_offset;
	}

	GMT_plotinit (argc, argv);

	GMT_map_clip_on (GMT_no_rgb, 3);

	if (Ctrl->T.active) {	/* Plot image as polygonal pieces and then exit */
		int np, index = 0, plot_n, n_use, start;
		double *grd_x, grd_y, *xx, *yy, *xp, *yp, dx2, dy2, y_bot, y_top;
		struct GMT_FILL *f;
		BOOLEAN skip;

		dx2 = 0.5 * g_head[0].x_inc;
		dy2 = 0.5 * g_head[0].y_inc;

		if (gmtdefs.verbose) fprintf (stderr, "%s: Tiling without interpolation\n", GMT_program);
		grd_x = (double *) GMT_memory (VNULL, (size_t)g_head[0].nx, sizeof (double), GMT_program);

		if (Ctrl->T.outline) GMT_setpen (&Ctrl->T.pen);
		for (i = 0; i < g_head[0].nx; i++) {
			grd_x[i] = GMT_i_to_x (i, g_head[0].x_min, g_head[0].x_max, g_head[0].x_inc, g_head[0].xy_off, g_head[0].nx);
			if (GMT_io.in_col_type[0] == GMT_IS_LON) {
				if (grd_x[i] < project_info.w && (grd_x[i] + 360.0) <= project_info.e) grd_x[i] += 360.0;
				if (grd_x[i] > project_info.e && (grd_x[i] - 360.0) >= project_info.w) grd_x[i] -= 360.0;
			}
		}
		for (j = node = 0; j < g_head[0].ny; j++) {
			grd_y = GMT_j_to_y (j, g_head[0].y_min, g_head[0].y_max, g_head[0].y_inc, g_head[0].xy_off, g_head[0].ny);
			y_bot = grd_y - dy2;
			y_top = grd_y + dy2;
			if (project_info.degree[1]) {
				if (y_bot < -90.0) y_bot = -90.0;
				if (y_top > +90.0) y_top = +90.0;
			}
			for (i = 0; i < g_head[0].nx; i++, node++) {	/* Compute rgb for each pixel */
				if (do_rgb) {
					for (k = 0, skip = FALSE; k < 3 && !skip; k++) {
						if (GMT_is_fnan (map[k][node])) {	/* If one is NaN they are all assumed to be NaN */
							if (Ctrl->T.skip)	/* Here we want to skip this node entirely */
								skip = TRUE;
							else {			/* Here we simply set the color to the NaN color */
								k = 3;
								memcpy ((void *)rgb, (void *)GMT_bfn[GMT_NAN].rgb, 3 * sizeof (int));
								index = -1;	/* Ensures no illumination done later */
							}
						}
						else {				/* Set color, let index = 0 so illuminate test will work */
							rgb[k] = irint (map[k][node]);	if (rgb[k] < 0) rgb[k] = 0; else if (rgb[k] > 255) rgb[k] = 255;	/* Clip */
							index = 0;
						}
					}
					if (skip) continue;
				}
				else {
					if (GMT_is_fnan (map[0][node]) && Ctrl->T.skip) continue;
					index = GMT_get_rgb24 (map[0][node], rgb);
					if (GMT_cpt_skip) continue;	/* Skip this z-slice */
				}
				if (Ctrl->I.active && index != -1) GMT_illuminate (intensity[node], rgb);
				if (Ctrl->M.active) rgb[0] = rgb[1] = rgb[2] = GMT_YIQ (rgb);

				np = GMT_graticule_path (&xx, &yy, 1, grd_x[i] - dx2, grd_x[i] + dx2, y_bot, y_top);

				plot_n = GMT_clip_to_map (xx, yy, np, &xp, &yp);
				GMT_free ((void *)xx);
				GMT_free ((void *)yy);
				if (plot_n == 0) continue;	/* Outside */

				if ((*GMT_will_it_wrap) (xp, yp, plot_n, &start)) {	/* Polygon wraps */

					/* First truncate agains left border */

					GMT_n_plot = (*GMT_truncate) (xp, yp, plot_n, start, -1);
					n_use = GMT_compact_line (GMT_x_plot, GMT_y_plot, GMT_n_plot, FALSE, 0);
					if (do_rgb)
						ps_polygon (GMT_x_plot, GMT_y_plot, n_use, rgb, Ctrl->T.outline);
					else if (index >= 0 && (f = GMT_lut[index].fill))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, Ctrl->T.outline);
					else if (index < 0 && (f = GMT_bfn[index+3].fill))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, Ctrl->T.outline);
					else
						ps_polygon (GMT_x_plot, GMT_y_plot, n_use, rgb, Ctrl->T.outline);

					/* Then truncate agains right border */

					GMT_n_plot = (*GMT_truncate) (xp, yp, plot_n, start, +1);
					n_use = GMT_compact_line (GMT_x_plot, GMT_y_plot, GMT_n_plot, FALSE, 0);
					if (do_rgb)
						ps_polygon (GMT_x_plot, GMT_y_plot, n_use, rgb, Ctrl->T.outline);
					else if (index >= 0 && (f = GMT_lut[index].fill))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, Ctrl->T.outline);
					else if (index < 0 && (f = GMT_bfn[index+3].fill))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, Ctrl->T.outline);
					else
						ps_polygon (GMT_x_plot, GMT_y_plot, n_use, rgb, Ctrl->T.outline);

				}
				else {
					if (do_rgb)
						ps_polygon (xp, yp, plot_n, rgb, Ctrl->T.outline);
					else if (index >= 0 && (f = GMT_lut[index].fill))
						GMT_fill (xp, yp, plot_n, f, Ctrl->T.outline);
					else if (index < 0 && (f = GMT_bfn[index+3].fill))
						GMT_fill (xp, yp, plot_n, f, Ctrl->T.outline);
					else
						ps_polygon (xp, yp, plot_n, rgb, Ctrl->T.outline);
				}
				GMT_free ((void *)xp);
				GMT_free ((void *)yp);
			}
		}

		GMT_free ((void *)grd_x);
		for (i = 0; i < n_grids; i++) GMT_free ((void *)map[i]);
		if (Ctrl->I.active) GMT_free ((void *)intensity);

		GMT_map_clip_off();

		GMT_map_basemap ();

		GMT_plotend ();

		GMT_end (argc, argv);
	}

	if (gmtdefs.verbose) {
		if (GMT_cpt_pattern) fprintf (stderr, "%s: Warning: Patterns in cpt file only apply to -T\n", GMT_program);
		fprintf (stderr, "%s: Evaluate pixel colors\n", GMT_program);
	}

	if (Ctrl->M.active || GMT_gray)
		bitimage_8 = (unsigned char *) GMT_memory (VNULL, (size_t)nm, sizeof (char), GMT_program);
	else {
		if (PS_colormask == -1) PS_colormask_off = 3;
		bitimage_24 = (unsigned char *) GMT_memory (VNULL, (size_t)(3 * nm + PS_colormask_off), sizeof (char), GMT_program);
		if (PS_colormask == -1) {
			bitimage_24[0] = (unsigned char)GMT_bfn[GMT_NAN].rgb[0];
			bitimage_24[1] = (unsigned char)GMT_bfn[GMT_NAN].rgb[1];
			bitimage_24[2] = (unsigned char)GMT_bfn[GMT_NAN].rgb[2];
		}
	}
	normal_x = !(project_info.projection == GMT_LINEAR && !project_info.xyz_pos[0]);
	normal_y = !(project_info.projection == GMT_LINEAR && !project_info.xyz_pos[1]);
	for (j = 0, byte = PS_colormask_off; j < r_head[0].ny; j++) {
		kk = r_head[0].nx * (normal_y ? j : r_head[0].ny - j - 1);
		for (i = 0; i < r_head[0].nx; i++) {	/* Compute rgb for each pixel */
			node = kk + (normal_x ? i : r_head[0].nx - i - 1);
			if (do_rgb) {
				for (k = 0; k < 3; k++) {
					if (GMT_is_fnan (map[k][node])) {	/* If one is NaN they are all assumed to be NaN */
						k = 3;
						memcpy ((void *)rgb, (void *)GMT_bfn[GMT_NAN].rgb, 3 * sizeof (int));
						index = -1;	/* Ensures no illumination done later */
					}
					else {				/* Set color, let index = 0 so illuminate test will work */
						rgb[k] = irint (map[k][node]);	if (rgb[k] < 0) rgb[k] = 0; else if (rgb[k] > 255) rgb[k] = 255;	/* Clip */
						index = 0;
					}
				}
			}
			else
				index = GMT_get_rgb24 (map[0][node], rgb);
			if (Ctrl->I.active && index != -1) GMT_illuminate (intensity[node], rgb);

			if (GMT_gray)	/* Color table only has grays, pick r */
				bitimage_8[byte++] = (unsigned char) rgb[0];
			else if (Ctrl->M.active)	/* Convert rgb to gray using the GMT_YIQ transformation */
				bitimage_8[byte++] = (unsigned char) GMT_YIQ (rgb);
			else {
				bitimage_24[byte++] = (unsigned char) rgb[0];
				bitimage_24[byte++] = (unsigned char) rgb[1];
				bitimage_24[byte++] = (unsigned char) rgb[2];
			}
		}
	}

	for (i = 0; i < n_grids; i++) GMT_free ((void *)map[i]);
	if (Ctrl->I.active) GMT_free ((void *)intensity);

	/* Get actual size of each pixel */

	dx = GMT_get_inc (r_head[0].x_min, r_head[0].x_max, r_head[0].nx, r_head[0].node_offset);
	dy = GMT_get_inc (r_head[0].y_min, r_head[0].y_max, r_head[0].ny, r_head[0].node_offset);

	/* Set lower left position of image on map */

	x0 = r_head[0].x_min;	y0 = r_head[0].y_min;
	if (grid_type == 0) {	/* Grid registration, move 1/2 pixel down/left */
		x0 -= 0.5 * dx;
		y0 -= 0.5 * dy;
	}

	x_side = dx * r_head[0].nx;
	y_side = dy * r_head[0].ny;

	if (gmtdefs.verbose) fprintf (stderr, "%s: Creating PostScript image ", GMT_program);

	if (GMT_gray) for (k = 0, GMT_b_and_w = TRUE; GMT_b_and_w && k < nm; k++) if (!(bitimage_8[k] == 0 || bitimage_8[k] == 255)) GMT_b_and_w = FALSE;

	if (GMT_b_and_w) {	/* Can get away with 1 bit image */
		int nx8, shift, byte, b_or_w, k8, nx_pixels;
		unsigned char *bit;

		if (gmtdefs.verbose) fprintf (stderr, "[1-bit B/W image]\n");

		nx8 = (int)ceil (r_head[0].nx / 8.0);
		nx_pixels = nx8 * 8;
		bit = (unsigned char *) GMT_memory (VNULL, (size_t)(nx8 * r_head[0].ny), sizeof (char), GMT_program);

		for (j = k = k8 = 0; j < r_head[0].ny; j++) {
			shift = byte = 0;
			for (i = 0; i < r_head[0].nx; i++, k++) {
				b_or_w = (bitimage_8[k] == 255);
				byte |= b_or_w;
				shift++;
				if (shift == 8) {	/* Time to dump out byte */
					bit[k8++] = (unsigned char) byte;
					byte = shift = 0;
				}
				else
					byte <<= 1;
			}
			if (shift) {
				byte |= 1;
				shift++;
				while (shift < 8) {
					byte <<= 1;
					byte |= 1;
					shift++;
				}
				bit[k8++] = (unsigned char) byte;
			}
		}
		GMT_free ((void *)bitimage_8);

		x_side = nx_pixels * dx;
		ps_bitimage (x0, y0, x_side, y_side, bit, nx_pixels, r_head[0].ny, FALSE, Ctrl->G.f_rgb, Ctrl->G.b_rgb);
		GMT_free ((void *)bit);
	}
	else if (GMT_gray || Ctrl->M.active) {
		if (gmtdefs.verbose) fprintf (stderr, "[8-bit grayshade image]\n");
		ps_colorimage (x0, y0, x_side, y_side, bitimage_8, r_head[0].nx, r_head[0].ny, 8 *PS_interpolate);
		GMT_free ((void *)bitimage_8);
	}
	else {
		if (gmtdefs.verbose) fprintf (stderr, "24-bit [%s]\n", c_method[gmtdefs.color_image]);
		GMT_color_image (x0, y0, x_side, y_side, bitimage_24, r_head[0].nx * PS_colormask, r_head[0].ny, 24 * PS_interpolate);
		GMT_free ((void *)bitimage_24);
	}

	GMT_map_clip_off();

	GMT_map_basemap ();

	GMT_plotend ();

	Free_Grdimage_Ctrl (Ctrl);	/* Deallocate control structure */

	GMT_end (argc, argv);

	return(0);
}

void GMT_set_proj_limits (struct GRD_HEADER *r, struct GRD_HEADER *g)
{
	/* Sets the projected extent of the grid given the map projection
	 * The extreme x/y coordinates are returned in r, and dx/dy, and
	 * nx/ny are set accordingly.  Not that some of these may change
	 * if GMT_grdproject_init is called at a later stage */

	int i, j;
	BOOLEAN all_lats, all_lons;
	double lon, lat, x, y;

	r->nx = g->nx;
	r->ny = g->ny;

	all_lats = (fabs (g->y_max - 90.0) < GMT_CONV_LIMIT && fabs (g->y_min + 90.0) < GMT_CONV_LIMIT);
	all_lons = (fabs (g->x_max - g->x_min - 360.0) < GMT_CONV_LIMIT);

	if (GMT_IS_AZIMUTHAL && all_lons && all_lats) {	/* Whole globe, get square box */
		r->x_min = project_info.xmin;	r->x_max = project_info.xmax;
		r->y_min = project_info.ymin;	r->y_max = project_info.ymax;
		return;
	}

	/* Must search for extent along perimeter */

	r->x_min = r->y_min = +DBL_MAX;
	r->x_max = r->y_max = -DBL_MAX;

	for (i = j = 0; i < g->nx; i++, j++) {	/* South and north */
		lon = g->x_min + i * g->x_inc;
		GMT_geo_to_xy (lon, g->y_min, &x, &y);
		r->x_min = MIN (r->x_min, x);	r->x_max = MAX (r->x_max, x);
		r->y_min = MIN (r->y_min, y);	r->y_max = MAX (r->y_max, y);
		GMT_geo_to_xy (lon, g->y_max, &x, &y);
		r->x_min = MIN (r->x_min, x);	r->x_max = MAX (r->x_max, x);
		r->y_min = MIN (r->y_min, y);	r->y_max = MAX (r->y_max, y);
	}
	for (i = 0; i < g->ny; j++, i++) {	/* East and west */
		lat = g->y_min + i * g->y_inc;
		GMT_geo_to_xy (g->x_min, lat, &x, &y);
		r->x_min = MIN (r->x_min, x);	r->x_max = MAX (r->x_max, x);
		r->y_min = MIN (r->y_min, y);	r->y_max = MAX (r->y_max, y);
		GMT_geo_to_xy (g->x_max, lat, &x, &y);
		r->x_min = MIN (r->x_min, x);	r->x_max = MAX (r->x_max, x);
		r->y_min = MIN (r->y_min, y);	r->y_max = MAX (r->y_max, y);
	}

	if (GMT_IS_AZIMUTHAL && all_lons) {	/* Full 360, use min/max for x */
		r->x_min = project_info.xmin;	r->x_max = project_info.xmax;
	}
	if (GMT_IS_AZIMUTHAL && all_lats) {	/* Full -90/+90, use min/max for y */
		r->y_min = project_info.ymin;	r->y_max = project_info.ymax;
	}
}

void *New_Grdimage_Ctrl () {	/* Allocate and initialize a new control structure */
	struct GRDIMAGE_CTRL *C;
	
	C = (struct GRDIMAGE_CTRL *) GMT_memory (VNULL, 1, sizeof (struct GRDIMAGE_CTRL), "New_Grdimage_Ctrl");
	
	/* Initialize values whose defaults are not 0/FALSE/NULL */
	
	GMT_init_pen (&C->T.pen, GMT_PENWIDTH);
	C->G.f_rgb[0] = C->G.f_rgb[1] = C->G.f_rgb[2] = 0;
	C->G.b_rgb[0] = C->G.b_rgb[1] = C->G.b_rgb[2] = 255;
	
	return ((void *)C);
}

void Free_Grdimage_Ctrl (struct GRDIMAGE_CTRL *C) {	/* Deallocate control structure */
	if (C->C.file) GMT_free ((void *)C->C.file);	
	if (C->I.file) GMT_free ((void *)C->I.file);	
	GMT_free ((void *)C);	
}
