/*--------------------------------------------------------------------
 *	$Id: grdimage.c,v 1.5.4.3 2002/01/17 19:04:41 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 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.
 *
 * Author:	Paul Wessel
 * Date:	20-SEP-1999
 * Ver:		3.0	based on old 2.x
 * Ver:		3.1	based on old 3.0
 * Ver:		3.3.2	Added -S as in grdproject
 * Ver:		3.3.3	Added -T[s] for non-interpolated tiling
 * Ver:		3.3.6	-T needed polygon clipping and wrapping (like psxy)
 * Ver:		3.4.1
 */
 
#include "gmt.h"

float *tmp1, *tmp2, *map, *intensity;
unsigned char *bitimage_8, *bitimage_24;

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

void GMT_limits (struct GRD_HEADER *r, struct GRD_HEADER *g);
void GMT_set_proj_limits (struct GRD_HEADER *r, struct GRD_HEADER *g);

main (int argc, char **argv)
{
	int i, j, k, kk, rgb[3], nm, nm2, byte, off, nx_f, ny_f, grid_type, polarity;
	int nx, ny, dpi = 0, nx_proj = 0, ny_proj = 0, node, f_rgb[3], tiling = 0, index;
	
	BOOLEAN error = FALSE, intens = FALSE, monochrome = FALSE;
	BOOLEAN mapped = FALSE, set_dpi = FALSE, subset = FALSE, transparent = FALSE;
	
	double  dx, dy, x_side, y_side, x0 = 0.0, y0 = 0.0, max_radius = 0.0;
	double west, east, south, north, data_west, data_east, data_south, data_north;
	
	char *grdfile, *intensfile, *cpt_file;
	
	struct GRD_HEADER g_head, r_head, i_head, j_head;
	
	argc = GMT_begin (argc, argv);

	GMT_grd_init (&g_head, argc, argv, FALSE);
	GMT_grd_init (&r_head, argc, argv, FALSE);
	
	grdfile = intensfile = cpt_file = 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 '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
			
				case 'C':
					cpt_file = &argv[i][2];
					break;
				case 'E':
					dpi = atoi (&argv[i][2]);
					set_dpi = TRUE;
					break;
				case 'G':	/* 1-bit fore or background color for transparent masks */
					transparent = TRUE;
					polarity = 1;
					j = 2;
					if (argv[i][2] == 'F' || argv[i][2] == 'f') j++;
					if (argv[i][2] == 'B' || argv[i][2] == 'b') polarity = 0, j++;
					if (GMT_getrgb (&argv[i][j], f_rgb)) {
						GMT_rgb_syntax ('G');
						error++;
					}
					break;
				case 'I':
					intens = TRUE;
					intensfile = &argv[i][2];
					break;
				case 'M':
					monochrome = TRUE;
					break;
				case '0':
					gmtdefs.color_image = 0;
					break;
				case 'S':
					max_radius = atof (&argv[i][2]);
					break;
				case 'T':
					tiling = 1;
					if (argv[i][2] == 's') tiling = 2;
					break;
				case '1':
					gmtdefs.color_image = 1;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			grdfile = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdimage %s - Plot grdfiles in 2-D\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdimage <grdfile> -C<cpt_file> -J<params> [-B<tickinfo>] [-E<dpi>] [-G[f|b]<rgb>]\n");
		fprintf (stderr, "\t[-I<intensity_file>] [-K] [-M] [-O] [-P] [-R<w/e/s/n>] [-S<radius>] [-T[s]] [-U] [-V]\n");
		fprintf (stderr, "\t[-X<x_shift>] [-Y<y_shift>] [-c<ncopies>]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr,"\t<grdfile> is data set to be plotted\n");
		fprintf (stderr,"\t-C color palette file\n");
		GMT_explain_option ('j');
		fprintf (stderr,"\n\tOPTIONS:\n");
		GMT_explain_option ('b');
		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-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');
		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");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		GMT_explain_option ('X');
		GMT_explain_option ('c');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!grdfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (!cpt_file) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify color palette table\n", GMT_program);
		error++;
	}
	if (intens && !intensfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify intensity file\n", GMT_program);
		error++;
	}
	if (set_dpi && dpi <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  dpi must be positive\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */
	
	/* Get color palette file */
	
	GMT_read_cpt (cpt_file);

	/* Check limits and get data file */
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Allocates memory and read data file\n", GMT_program);
	
	if (GMT_read_grd_info (grdfile, &g_head)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}
	off = (g_head.node_offset) ? 0 : 1;

	/* Determine what wesn to pass to map_setup */

	if (!project_info.region_supplied) {
		west = g_head.x_min;
		east = g_head.x_max;
		south = g_head.y_min;
		north = g_head.y_max;
	}
	else if (!(west == g_head.x_min && east == g_head.x_max && south == g_head.y_min && north == g_head.y_max))
		subset = TRUE;

	GMT_map_setup (west, east, south, north);

	/* Determine the wesn to be used to read the grdfile */

	GMT_grd_setregion (&g_head, &data_west, &data_east, &data_south, &data_north);

	nx_f = g_head.nx;
	ny_f = g_head.ny;
	
	/* Read data */

	nx = irint ( (data_east - data_west) / g_head.x_inc) + off;
	ny = irint ( (data_north - data_south) / g_head.y_inc) + off;
	nm = nx * ny;
	tmp1 = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
	if (GMT_read_grd (grdfile, &g_head, tmp1, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}

	/* If given, get intensity file or compute intensities */
	
	if (intens) {	/* 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 (intensfile, &i_head)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, intensfile);
			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 (intensfile, &i_head, tmp2, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, intensfile);
			exit (EXIT_FAILURE);
		}
	}
	
	GMT_set_proj_limits (&r_head, &g_head);

	/* if (dpi > 0 || MAPPING) { */
	if (!tiling && (dpi > 0 || project_info.projection > 5)) {	/* Need to resample the grd file */
		
		mapped = TRUE;
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: project grdfiles\n", GMT_program);
	
		if (dpi == 0) {	/* Use input # of nodes as # of projected nodes */
			nx_proj = g_head.nx;
			ny_proj = g_head.ny;
		}
		grid_type = (dpi > 0) ? 1 : g_head.node_offset;	/* Force pixel if dpi is set */
		GMT_grdproject_init (&r_head, 0.0, 0.0, nx_proj, ny_proj, dpi, grid_type);
		GMT_init_search_radius (&max_radius, &r_head, &g_head, FALSE);
		nm2 = r_head.nx * r_head.ny;
		map = (float *) GMT_memory (VNULL, (size_t)nm2, sizeof (float), "grdproject");

		GMT_grd_forward (tmp1, &g_head, map, &r_head, max_radius);
		
		GMT_free ((void *)tmp1);
		if (intens) {
			j_head.x_min = r_head.x_min;	j_head.x_max = r_head.x_max;
			j_head.y_min = r_head.y_min;	j_head.y_max = r_head.y_max;
		
			if (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, dpi, grid_type);
			intensity = (float *) GMT_memory (VNULL, (size_t)nm2, sizeof (float), "grdproject");
			GMT_grd_forward (tmp2, &i_head, intensity, &j_head, max_radius);
			GMT_free ((void *)tmp2);
		}
		nm = nm2;
	}
	else {	/* Simply copy g_head info to r_head */
		map = tmp1;

		r_head.nx = g_head.nx;		r_head.ny = g_head.ny;
		r_head.x_inc = g_head.x_inc;	r_head.y_inc = g_head.y_inc;
		if (intens) {
			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.node_offset;
	}
	
	ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
		gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
		gmtdefs.dpi, GMT_INCH , gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo (argv[0]));
	GMT_echo_command (argc, argv);
	if (gmtdefs.unix_time) GMT_timestamp (argc, argv);
	
	GMT_map_clip_on (GMT_no_rgb, 3);

	if (tiling) {	/* Plot image as polygonal pieces and then exit */
		int np, index, plot_n, n_use, start;
		double *grd_x, grd_y, *xx, *yy, *xp, *yp, d, dx2, dy2, y_bot, y_top;
		struct GMT_FILL *f;
		
		d = (off) ? 0.0 : 0.5;
		dx2 = 0.5 * g_head.x_inc;
		dy2 = 0.5 * g_head.y_inc;
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Tiling without interpolation\n", GMT_program);
		grd_x = (double *) GMT_memory (VNULL, (size_t)g_head.nx, sizeof (double), GMT_program);
		
		for (i = 0; i < g_head.nx; i++) grd_x[i] = g_head.x_min + (i + d) * g_head.x_inc;
		for (j = node = 0; j < g_head.ny; j++) {
			grd_y = g_head.y_max - (j + d) * g_head.y_inc;
			y_bot = grd_y - dy2;
			y_top = grd_y + dy2;
			if (MAPPING) {
				if (y_bot < -90.0) y_bot = -90.0;
				if (y_top > +90.0) y_top = +90.0;
			}
			for (i = 0; i < g_head.nx; i++, node++) {	/* Compute rgb for each pixel */
				if (GMT_is_fnan (map[node]) && tiling == 2) continue;
				index = GMT_get_rgb24 (map[node], rgb);
				if (GMT_cpt_skip) continue;	/* Skip this z-slice */
				if (intens && index != -1) GMT_illuminate (intensity[node], rgb);
				if (monochrome) rgb[0] = rgb[1] = rgb[2] = 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 (index >= 0 && (f = GMT_lut[index].fill))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, FALSE);
					else if (index < 0 && (f = GMT_bfn.fill[index+3]))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, FALSE);
					else
						ps_polygon (GMT_x_plot, GMT_y_plot, n_use, rgb, FALSE);
													
					/* 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 (index >= 0 && (f = GMT_lut[index].fill))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, FALSE);
					else if (index < 0 && (f = GMT_bfn.fill[index+3]))
						GMT_fill (GMT_x_plot, GMT_y_plot, n_use, f, FALSE);
					else
						ps_polygon (GMT_x_plot, GMT_y_plot, n_use, rgb, FALSE);
					
				}
				else {
					if (index >= 0 && (f = GMT_lut[index].fill))
						GMT_fill (xp, yp, plot_n, f, FALSE);
					else if (index < 0 && (f = GMT_bfn.fill[index+3]))
						GMT_fill (xp, yp, plot_n, f, FALSE);
					else
						ps_polygon (xp, yp, plot_n, rgb, FALSE);
				}
				GMT_free ((void *)xp);
				GMT_free ((void *)yp);
			}
		}
		
		GMT_free ((void *)grd_x);
	
		GMT_map_clip_off();
	
		if (frame_info.plot) {
			ps_setpaint (gmtdefs.basemap_frame_rgb);
			GMT_map_basemap ();
		}
		ps_plotend (gmtdefs.last_page);
	
		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 (monochrome || GMT_gray)
		bitimage_8 = (unsigned char *) GMT_memory (VNULL, (size_t)nm, sizeof (char), GMT_program);
	else
		bitimage_24 = (unsigned char *) GMT_memory (VNULL, (size_t)(3 * nm), sizeof (char), GMT_program);
	
	for (j = byte = 0; j < r_head.ny; j++) {
		kk = r_head.nx * (project_info.xyz_pos[1] ? j : r_head.ny - j - 1);
		for (i = 0; i < r_head.nx; i++) {	/* Compute rgb for each pixel */
			node = kk + (project_info.xyz_pos[0] ? i : r_head.nx - i - 1);
			index = GMT_get_rgb24 (map[node], rgb);
			if (intens && 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 (monochrome)	/* Convert rgb to gray using the YIQ transformation */
				bitimage_8[byte++] = (unsigned char) YIQ (rgb);
			else {
				bitimage_24[byte++] = (unsigned char) rgb[0];
				bitimage_24[byte++] = (unsigned char) rgb[1];
				bitimage_24[byte++] = (unsigned char) rgb[2];
			}
		}
	}
	
	GMT_free ((void *)map);
	if (intens) GMT_free ((void *)intensity);
	
	/* Get actual size of each pixel */

	dx = (r_head.x_max - r_head.x_min) / (r_head.nx - off);
	dy = (r_head.y_max - r_head.y_min) / (r_head.ny - off);

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

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

	if (subset && GRID_CLIP_OK) GMT_grid_clip_on (&g_head, GMT_no_rgb, 3);

	x_side = dx * r_head.nx;
	y_side = dy * r_head.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;
		
		nx8 = (int)ceil (r_head.nx / 8.0);
		nx_pixels = nx8 * 8;
		bit = (unsigned char *) GMT_memory (VNULL, (size_t)(nx8 * r_head.ny), sizeof (char), GMT_program);
		
		for (j = k = k8 = 0; j < r_head.ny; j++) {
			shift = byte = 0;
			for (i = 0; i < r_head.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;
		if (transparent) {
			if (gmtdefs.verbose) fprintf (stderr, "[1-bit image mask painted %d/%d/%d]\n", f_rgb[0], f_rgb[1], f_rgb[2]);
			ps_imagemask (x0, y0, x_side, y_side, bit, nx_pixels, r_head.ny, polarity, f_rgb);
		}
		else {
			if (gmtdefs.verbose) fprintf (stderr, "[1-bit B/W image]\n");
			ps_image (x0, y0, x_side, y_side, bit, nx_pixels, r_head.ny, 1);
		}
		GMT_free ((void *)bit);
	}
	else if (GMT_gray || monochrome) {
		if (gmtdefs.verbose) fprintf (stderr, "[8-bit grayshade image]\n");
		ps_image (x0, y0, x_side, y_side, bitimage_8, r_head.nx, r_head.ny, 8);
		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.nx, r_head.ny);
		GMT_free ((void *)bitimage_24);
	}
		
	if (subset && GRID_CLIP_OK) GMT_grid_clip_off();

	GMT_map_clip_off();
	
	if (frame_info.plot) {
		ps_setpaint (gmtdefs.basemap_frame_rgb);
		GMT_map_basemap ();
	}
	ps_plotend (gmtdefs.last_page);
	
	GMT_end (argc, argv);
}

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

	/* Then truncate, if neccesary, values to fit inside plotbox */
	
	r->x_min = MAX (r->x_min, project_info.xmin);	r->x_max = MIN (r->x_max, project_info.xmax);
	r->y_min = MAX (r->y_min, project_info.ymin);	r->y_max = MIN (r->y_max, project_info.ymax);

	r->x_inc = (r->x_max - r->x_min) / (r->nx - r->node_offset);
	r->y_inc = (r->y_max - r->y_min) / (r->ny - r->node_offset);
}