/*--------------------------------------------------------------------
 *    The GMT-system:	@(#)grdview.c	2.111  11/29/99
 *
 *	Copyright (c) 1991-1999 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: www.soest.hawaii.edu/gmt
 *--------------------------------------------------------------------*/
/*
 * grdview will read a topofile and produce a 3-D perspective plot
 * of the surface z = f(x,y) using PostScript. The surface can
 * be represented as:
 *	1) A Mesh plot
 *	2) A shaded (or colored) surface w/wo contourlines and w/wo
 *	   illumination by artificial sun(s).
 *
 * grdview calls contours to find the line segments that make up the
 * contour lines. This allows the user to specify that the contours
 * should be smoothed before plotting. This will make the resulting
 * image smoother, especially if nx and ny are relatively small.
 *
 * Author:	Paul Wessel
 * Date:	8-DEC-1988
 * Modified:	5-JUL-1994	New version 3.0
 * Modified:	15-JUL-1998	New version 3.1
 * Modified:	15-MAR-1999	New version 3.2
 * Modified:	07-OCT-1999	3.3.2 - Bug in image fixed, plus added correct
 *				image mask clipping
 *
 */
 
#include "gmt.h"
#include "gmt_boundcond.h"
#include "gmt_bcr.h"

/* Declarations needed for binning of smooth contours */

struct BIN {
	struct CONT *first_cont;
} *binij;

struct CONT {
	struct POINT *first_point;
	struct CONT *next_cont;
	double value;
};

struct POINT {
	double x, y;
	struct POINT *next_point;
};

/* Global variables */

float *grd, *intensity, *topo;

int *edge;
int bin_inc[4], ij_inc[4];
double x_inc[4], y_inc[4], *x, *y, *z, *v, *xx, *yy;

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

void grdview_init_setup(struct GRD_HEADER *header, float *topo, int two, BOOLEAN draw_plane, double plane_level);
int pixel_inside(int ip, int jp, int *ix, int *iy, int bin), quick_idist(int x1, int y1, int x2, int y2);

main (int argc, char **argv)
{
	int 	i, j, ij, n_edges, k, k1, n, max, i_bin, j_bin, i_bin_old, j_bin_old;
	int	smooth = 1, side, way, nm, nx_f, ny_f, off, nx, ny;
	int	rgb_facade[3], bin, two, mx, my, tiling = 0;
	int	c, rgb[3], i_start, i_stop, j_start, j_stop, i_inc, j_inc;
	int	n_saddle, prev_side, entry, check, n_added, dpi_i = 100;
	int	sw, se, nw, ne, index, id, n4, nk;
	
	BOOLEAN mesh = FALSE, draw_plane = FALSE, facade = FALSE, get_contours, bad, set_z = FALSE;
	BOOLEAN error = FALSE, outline = FALSE, surface = FALSE, pen_not_set, bilinear = FALSE;
	BOOLEAN first, begin, draw_contours = FALSE, moved[4], intens = FALSE, drape = FALSE;
	BOOLEAN saddle, cross_ok, go, image = FALSE, monochrome = FALSE, subset = FALSE;
	
	double cval, x_left, x_right, y_top, y_bottom, small, z_ave, this_intensity, *xval, *yval;
	double plane_level = 0.0, dx2, dy2, del_x, del_y, take_out, get_intensity(float *intensity, int k, int nx);
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, new_z_level = 0.0, del_off;
	double data_west, data_east, data_south, data_north, delx, dely, z_val, next_up, xmesh[4], ymesh[4];
	double x_pixel_size, y_pixel_size, *x_imask, *y_imask;
	
	float was[4];
	
	char *topofile, *intensfile, *drapefile, *cpt_file;
	
	struct CONT *get_cont_struct(int bin, double value), *start_cont, *this_cont, *last_cont;
	struct POINT *get_point(double x, double y), *this_point, *last_point;
	
	struct GRD_HEADER header, t_head, i_head, d_head;
	
	struct GMT_PEN pen[2];
	struct GMT_FILL fill;

        struct GMT_EDGEINFO edgeinfo;

	for (i = 0; i < 3; i++) rgb_facade[i] = 0;
	
	cpt_file = topofile = intensfile = drapefile = CNULL;
	
	argc = GMT_begin (argc, argv);
	
	GMT_init_pen (&pen[0], 3.0 * GMT_PENWIDTH);	/* Contour pen */
	GMT_init_pen (&pen[1], GMT_PENWIDTH);	/* Mesh pen */
	GMT_init_fill (&fill, 255 - gmtdefs.basemap_frame_rgb[0], 255 - gmtdefs.basemap_frame_rgb[0], 255 - gmtdefs.basemap_frame_rgb[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':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'G':
					drapefile = &argv[i][2];
					drape = TRUE;
					break;
				case 'I':
					intensfile = &argv[i][2];
					intens = TRUE;
					break;
				case 'L':
					if (argv[i][2])
						error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]);
					else
						bilinear = TRUE;
					break;
				case 'N':
					if (argv[i][2]) {
						n = sscanf (&argv[i][2], "%lf/%d/%d/%d", &plane_level, &rgb_facade[0], &rgb_facade[1], &rgb_facade[2]);
						draw_plane = TRUE;
						if (n == 4) facade = TRUE;
						if ((n > 1 && n != 4) || (n == 4 && GMT_check_rgb (rgb_facade))) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR option -N:  RGB components must all be in 0-255 range\n", GMT_program);
							error = TRUE;
						}
					}
					else {
						fprintf (stderr, "%s: GMT SYNTAX ERROR option -N:  Must specify plane level\n", GMT_program);
						error = TRUE;
					}
					break;
				case 'Q':
					switch (argv[i][2]) {
						case 'm':	/* Mesh plot */
							mesh = TRUE;
							if (argv[i][3] == '/' && GMT_getfill (&argv[i][4], &fill)) {
								fprintf (stderr, "%s: GMT SYNTAX ERROR -Qm option: RGB components must all be in 0-255 range\n", GMT_program);
								error = TRUE;
							}
							break;
						case 's':	/* Color wo/ contours */
							surface = TRUE;
							if (argv[i][3] == 'm') outline = TRUE;
							break;
						case 'i':	/* image w/ clipmask */
						case 'I':	/* Backward compatibility, gives -Qi */
							image = TRUE;
							if (argv[i][3]) dpi_i = atoi (&argv[i][3]);
							break;
						default:
							fprintf (stderr, "%s: GMT SYNTAX ERROR:  Unrecognized qualifier (%c) for option -%c\n", GMT_program, argv[i][2], argv[i][1]);
							error = TRUE;
							break;
					}
					monochrome = (argv[i][strlen(argv[i])-1] == 'g');
					break;
				case 'S':
					smooth = atoi (&argv[i][2]);
					break;
				case 'T':
					tiling = 1;
					if (argv[i][2] == 's') tiling = 2;
					break;
				case 'W':	/* Contour and/or mesh pens */
					j = (argv[i][2] == 'm' || argv[i][2] == 'c') ? 3 : 2;
                                        if (j == 2) {   /* Set both */
						if (GMT_getpen (&argv[i][2], &pen[0])) {
							GMT_pen_syntax ('W');
							error++;
						}
						else
							pen[1] = pen[0];
						draw_contours = TRUE;
					}
					else {
						id = (argv[i][2] == 'm') ? 1 : 0;
						if (GMT_getpen (&argv[i][j], &pen[id])) {
							GMT_pen_syntax ('W');
							error++;
						}
						if (id == 0) draw_contours = TRUE;

					}
					break;
				case 'Z':
					if (argv[i][2]) {
						new_z_level = atof (&argv[i][2]);
						set_z = TRUE;
					}
					break;
					
				/* Illegal options */
			
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			topofile = argv[i];
	}

	if (!(mesh || surface || image)) mesh = TRUE;	/* Default is mesh plot */
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdview %s - Plot topofiles in 3-D\n\n", GMT_VERSION);
		fprintf (stderr,"usage: grdview <topofile> -J<params> [-B<tickinfo>] [-C<cpt_file>] [-E<v_az/v_el>]\n");
		fprintf (stderr,"	[-G<drapefile>] [-I<intensfile>] [-K] [-L[<flags>]] [-N<level>[/<r/g/b>]] [-O] [-P]\n");
		fprintf (stderr,"	[-Q<type>] [-R<region>] [-S<smooth>] [-T[s]] [-U] [-V] [-W<type><pen>]\n");
		fprintf (stderr,"	[-X<x_shift>] [-Y<y_shift>] [-Z<z_level>] [-c<ncopies>]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr,"	<topofile> is data set to be plotted\n");
		GMT_explain_option ('j');
		fprintf (stderr,"\n\tOPTIONS:\n");
		GMT_explain_option ('b');
		fprintf (stderr,"	-C Color palette file\n");
		fprintf (stderr,"	-E Draw perspective from viewpoint at azimuth <v_az>, elevation <v_el> (degrees).\n");
		fprintf (stderr,"	   Default is looking straight down [180/90].\n");
		fprintf (stderr,"	-G <drapefile> rather than <topofile> is the data set to color-code.\n");
		fprintf (stderr,"	   Use <topofile> as the relief and \'drape\' the image on top.\n");
		fprintf (stderr,"	   Note that -Jz and -N always refers to the <topofile>\n");
		fprintf (stderr,"	-I gives name of intensity file and selects illumination\n");
		fprintf (stderr, "      -L sets boundary conditions when resampling the grid.  <flags> can be either\n");
		fprintf (stderr, "         g for geographic boundary conditions\n");
		fprintf (stderr, "         or one or both of\n");
		fprintf (stderr, "         x for periodic boundary conditions on x\n");
		fprintf (stderr, "         y for periodic boundary conditions on y\n");
		fprintf (stderr,"	   If no <flags> are set, use bilinear rather than bicubic [Default] resampling \n");
		GMT_explain_option ('Z');
		GMT_explain_option ('K');
		fprintf (stderr,"	-N<level> will draw a plane at z = level.  Append color [/r/g/b] to paint\n");
		fprintf (stderr,"	   the facade between the plane and the data perimeter.\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr,"	-Q sets plot reQuest. Choose one of the following:\n");
		fprintf (stderr,"	   -Qm for Mesh plot [Default].  Append color for mesh paint [%d/%d/%d]\n", fill.rgb[0], fill.rgb[1], fill.rgb[2]);
		fprintf (stderr,"	   -Qs[m] for colored or shaded Surface.  Append m to draw meshlines on the surface.\n");
		fprintf (stderr,"	   -Qi for scanline converting polygons to rasterimage.  Append effective dpi [100].\n");
		fprintf (stderr,"	  To force a monochrome image using the YIQ transformation, append g\n");
		GMT_explain_option ('R');
		fprintf (stderr,"	-S will smooth contours first (see grdcontour for <smooth> value info).\n");
		fprintf(stderr, "	-T will image the data without interpolation by painting polygonal tiles\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		fprintf (stderr,"	-W sets pen attributes. <type> can be c for contours or m for mesh\n");
		fprintf (stderr,"	   -Wc draw scontours on top of surface or mesh.  [Default is no contours]\n");
		fprintf (stderr, "	       Optionally append pen attributes [width = %lgp, color = (%d/%d/%d), solid line].\n", 
			pen[0].width, pen[0].rgb[0], pen[0].rgb[1], pen[0].rgb[2]);
		fprintf (stderr,"	   -Wm sets attributes for mesh lines [[width = %lgp, color = (%d/%d/%d), solid line].\n", 
			pen[1].width, pen[1].rgb[0], pen[1].rgb[1], pen[1].rgb[2]);
		fprintf (stderr,"	       Requires -Qm or -Qsm to take effect.\n");
		GMT_explain_option ('X');
		fprintf (stderr, "      -Z For 3-D plots: Set the z-level of map [0]\n");
		GMT_explain_option ('c');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
		
	if (!topofile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (drape && !drapefile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -G:  Must specify drape file\n", GMT_program);
		error++;
	}
	if (intens && !intensfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -I:  Must specify intensity file\n", GMT_program);
		error++;
	}
	if ((surface || image || draw_contours) && !cpt_file) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify color palette table\n", GMT_program);
		error++;
	}
	if (image && dpi_i <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Qi option:  Must specify positive dpi\n", GMT_program);
		error++;
	}
	if (smooth <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  smooth value must be positive\n", GMT_program);
		error++;
	}
	if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */
	
	if (cpt_file) {
		GMT_read_cpt (cpt_file);
		if (GMT_b_and_w) monochrome = TRUE;
	}
	
	get_contours = ( (mesh && draw_contours) || (surface && GMT_n_colors > 1) );
	
	two = (drape && get_contours) ? 2 : 0;	/* Must read topofile with 2 boundary columns/rows */
	
	if (!strcmp (topofile, "=")) {
		fprintf (stderr, "%s: Piping of topofile not supported!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (GMT_read_grd_info (topofile, &t_head)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, topofile);
		exit (EXIT_FAILURE);
	}
	
	if (intens && GMT_read_grd_info (intensfile, &i_head)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, intensfile);
		exit (EXIT_FAILURE);
	}
	
	if (drape && GMT_read_grd_info (drapefile, &d_head)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, drapefile);
		exit (EXIT_FAILURE);
	}
	
	off = (t_head.node_offset) ? 0 : 1;
	del_off = (t_head.node_offset) ? 0.5 : 0.0;
	
	GMT_boundcond_init (&edgeinfo);
	
	/* Determine what wesn to pass to map_setup */

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

	if (project_info.z_bottom == 0.0 && project_info.z_top == 0.0) {
		project_info.z_bottom = t_head.z_min;
		project_info.z_top = t_head.z_max;
		if (draw_plane && plane_level < project_info.z_bottom) project_info.z_bottom = plane_level;
		if (draw_plane && plane_level > project_info.z_top) project_info.z_top = plane_level;
	}
	
	GMT_map_setup (west, east, south, north);

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

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

	/* Read data */

	if (gmtdefs.verbose) fprintf (stderr, "%s: Processing shape file\n", GMT_program);
	
	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = two;
	nx = irint ( (data_east - data_west) / t_head.x_inc) + off;
	ny = irint ( (data_north - data_south) / t_head.y_inc) + off;
	mx = nx + 2 * two;
	my = ny + 2 * two;
	nm = nx * ny;
	topo = (float *) GMT_memory (VNULL, (size_t)(mx * my), sizeof (float), GMT_program);
	if (GMT_read_grd (topofile, &t_head, topo, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, topofile);
		exit (EXIT_FAILURE);
	}

	nx_f = t_head.nx;
	ny_f = t_head.ny;

	delx = (t_head.node_offset) ? 0.5 * t_head.x_inc :0.0;
	dely = (t_head.node_offset) ? 0.5 * t_head.y_inc :0.0;
	
	xval = (double *) GMT_memory (VNULL, (size_t)t_head.nx, sizeof (double), GMT_program);
	yval = (double *) GMT_memory (VNULL, (size_t)t_head.ny, sizeof (double), GMT_program);
	for (i = 0; i < t_head.nx; i++) xval[i] = t_head.x_min + i * t_head.x_inc + delx;
	for (j = 0; j < t_head.ny; j++) yval[j] = t_head.y_max - j * t_head.y_inc - dely;
		
	GMT_boundcond_param_prep (&t_head, &edgeinfo);

	if (drape) {	/* draping wanted */
	
		if (gmtdefs.verbose) fprintf (stderr, "%s: Processing drape file\n", GMT_program);
	
		grd = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
		if (GMT_read_grd (drapefile, &d_head, grd, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, drapefile);
			exit (EXIT_FAILURE);
		}
		if (d_head.nx != nx_f || d_head.ny != ny_f) {
			fprintf (stderr, "%s: Drape file has improper dimensions!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		header = d_head;
	}
	else {
		grd = topo;
		header = t_head;
	}
	
	if (!project_info.xyz_pos[2])	/* Negative z-scale, must flip */
		d_swap (project_info.z_bottom, project_info.z_top);

	ij_inc[0] = 0;		ij_inc[1] = 1;	ij_inc[2] = 1 - mx;	ij_inc[3] = -mx;
	bin_inc[0] = 0;		bin_inc[1] = 1;	bin_inc[2] = 1 - header.nx;	bin_inc[3] = -header.nx;
	nw = two * mx + two;
	ne = nw + t_head.nx - 1;
	sw = (t_head.ny + two - 1) * mx + two;
	se = sw + t_head.nx - 1;
	
	grdview_init_setup (&t_head, topo, two, draw_plane, plane_level);	/* Find projected min/max in y-direction */
			
	i_start = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 0 : header.nx - 2;
	i_stop  = (z_project.quadrant == 1 || z_project.quadrant == 2) ? header.nx - 1 : -1;
	i_inc   = (z_project.quadrant == 1 || z_project.quadrant == 2) ? 1 : -1;
	j_start = (z_project.quadrant == 1 || z_project.quadrant == 4) ? header.ny - 1 : 1;
	j_stop  = (z_project.quadrant == 1 || z_project.quadrant == 4) ? 0 : header.ny;
	j_inc   = (z_project.quadrant == 1 || z_project.quadrant == 4) ? -1 : 1;
	x_inc[0] = x_inc[3] = 0.0;	x_inc[1] = x_inc[2] = header.x_inc;
	y_inc[0] = y_inc[1] = 0.0;	y_inc[2] = y_inc[3] = header.y_inc;

	if (get_contours) {	/* Need to find contours */
		if (gmtdefs.verbose) fprintf (stderr, "%s: Find contours\n", GMT_program);
		n_edges = header.ny * (int )ceil (header.nx / 16.0);
		edge = (int *) GMT_memory (VNULL, (size_t)n_edges, sizeof (int), GMT_program);
		binij = (struct BIN *) GMT_memory (VNULL, (size_t)nm, sizeof (struct BIN), GMT_program);
		small = SMALL * (header.z_max - header.z_min);
		if (gmtdefs.verbose) fprintf (stderr, "%s: Trace and bin contours...\n", GMT_program);
		first = TRUE;
		for (c = 0; c < GMT_n_colors+1; c++) {	/* For each color change */
		
			/* Reset markers and set up new zero-contour*/
		
			cval = (c == GMT_n_colors) ? GMT_lut[c-1].z_high : GMT_lut[c].z_low;

			if (cval < header.z_min || cval > header.z_max) continue;

			if (gmtdefs.verbose) fprintf (stderr, "%s: Now tracing contour interval %8lg\r", GMT_program, cval);
			take_out = (first) ? cval : cval - GMT_lut[c-1].z_low;
			first = FALSE;
			for (i = 0; i < nm; i++) {
				if (!GMT_is_fnan (grd[i])) grd[i] -= (float)take_out;
				if (grd[i] == 0.0) grd[i] += (float)small;
			}
		
			side = 0;
			begin = TRUE;
			while (side < 5) {
				while ((n = GMT_contours (grd, &header, smooth, gmtdefs.interpolant, &side, edge, begin, &x, &y)) > 0) {
		
					begin = FALSE;
		
					i_bin_old = j_bin_old = -1;
					for (i = 1; i < n; i++) {
						i_bin = (int)floor (((0.5 * (x[i-1] + x[i]) - header.x_min) / header.x_inc) - del_off);
						j_bin = (int)floor (((header.y_max - 0.5 * (y[i-1] + y[i])) / header.y_inc) - del_off) + 1;
						if (i_bin != i_bin_old || j_bin != j_bin_old) {	/* Entering new bin */
							bin = j_bin * header.nx + i_bin;
							this_cont = get_cont_struct (bin, cval);
							this_cont->value = cval;
							this_cont->first_point = get_point (x[i-1], y[i-1]);
							this_point = this_cont->first_point;
							i_bin_old = i_bin;
							j_bin_old = j_bin;
						}
						this_point->next_point = get_point (x[i], y[i]);
						this_point = this_point->next_point;
					}
					GMT_free ((void *)x);
					GMT_free ((void *)y);
				}
				begin = FALSE;
			}
		}
	
		/* Remove temporary variables */
		
		GMT_free ((void *)edge);
		
		/* Go back to beginning and reread since grd has been destroyed */
		
		if (drape) {
			if (GMT_read_grd_info (drapefile, &d_head)) {
				fprintf (stderr, "%s: Error opening file %s\n", GMT_program, drapefile);
				exit (EXIT_FAILURE);
			}
			if (GMT_read_grd (drapefile, &d_head, grd, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error reading file %s\n", GMT_program, drapefile);
				exit (EXIT_FAILURE);
			}
			header = d_head;
		}
		else {
			if (GMT_read_grd_info (topofile, &t_head)) {
				fprintf (stderr, "%s: Error opening file %s\n", GMT_program, topofile);
				exit (EXIT_FAILURE);
			}
			if (GMT_read_grd (topofile, &t_head, topo, data_west, data_east, data_south, data_north, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error reading file %s\n", GMT_program, topofile);
				exit (EXIT_FAILURE);
			}
			grd = topo;
			header = t_head;
		}
		if (gmtdefs.verbose) fprintf (stderr, "\n");
	}
	
	if (intens) {	/* Illumination wanted */
	
		if (gmtdefs.verbose) fprintf (stderr, "%s: Processing illumination file\n", GMT_program);
	
		intensity = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
		if (GMT_read_grd (intensfile, &i_head, intensity, 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);
		}
		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);
		}
	}
	
	if (two) {	/* Initialize bcr stuff */
		GMT_bcr_init (&t_head, GMT_pad, bilinear);

		/* Set boundary conditions  */

		GMT_boundcond_set (&t_head, &edgeinfo, GMT_pad, topo);
	
	}

	for (j = bin = ij = 0; j < header.ny; j++) {	/* Nodes outside -R is set to NaN and not used */
		ij = (two) ? (j + 2) * mx + two : bin;
		for (i = 0; i < header.nx; i++, bin++, ij++) if (GMT_map_outside (xval[i], yval[j])) topo[ij] = GMT_f_NaN;
	}
	
	max = 2 * (smooth * (((header.nx > header.ny) ? header.nx : header.ny) + 2)) + 1;
	x = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program);
	y = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program);
	z = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program);
	v = (double *) GMT_memory (VNULL, (size_t)max, sizeof (double), GMT_program);
	
	dx2 = 0.5 * header.x_inc;		dy2 = 0.5 * header.y_inc;
			
	if (gmtdefs.verbose) fprintf (stderr, "%s: Start creating PostScript plot\n", GMT_program);

	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]));

	ps_setformat (3);
	GMT_echo_command (argc, argv);
	if (gmtdefs.unix_time) GMT_timestamp (argc, argv);

	if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);

	if (frame_info.plot) GMT_map_basemap ();	/* Plot basemap */

	if (project_info.z_pars[0] == 0.0) GMT_map_clip_on (GMT_no_rgb, 3);

	if (set_z) project_info.z_level = new_z_level;
	
	max = 2 * (smooth * (((header.nx > header.ny) ? header.nx : header.ny) + 2)) + 1;
	xx = (double *) GMT_memory (VNULL, (size_t)max, sizeof(double), GMT_program);
	yy = (double *) GMT_memory (VNULL, (size_t)max, sizeof(double), GMT_program);
	if (draw_plane) {
		ps_comment ("Plot the plane at desired level");
		if (!z_project.draw[0])	{	/* Southern side */
			for (i = 0; i < header.nx; i++) GMT_geoz_to_xy (header.x_min + i * header.x_inc + delx, header.y_min + dely, plane_level, &xx[i], &yy[i]);
			ps_line (xx, yy, header.nx, 3, TRUE, TRUE);
		}
		if (!z_project.draw[2])	{	/* Northern side */
			for (i = 0; i < header.nx; i++) GMT_geoz_to_xy (header.x_min + i * header.x_inc + delx, header.y_max - dely, plane_level, &xx[i], &yy[i]);
			ps_line (xx, yy, header.nx, 3, TRUE, TRUE);
		}
		if (!z_project.draw[3])	{	/* Western side */
			for (j = 0; j < header.ny; j++) GMT_geoz_to_xy (header.x_min + delx, header.y_max - j * header.y_inc - dely, plane_level, &xx[j], &yy[j]);
			ps_line (xx, yy, header.ny, 3, TRUE, TRUE);
		}
		if (!z_project.draw[1])	{	/* Eastern side */
			for (j = 0; j < header.ny; j++) GMT_geoz_to_xy (header.x_max - delx, header.y_max - j * header.y_inc - dely, plane_level, &xx[j], &yy[j]);
			ps_line (xx, yy, header.ny, 3, TRUE, TRUE);
		}
		
		GMT_geoz_to_xy (header.x_min + delx, header.y_min + dely, plane_level, &xx[0], &yy[0]);
		GMT_geoz_to_xy (header.x_max - delx, header.y_min + dely, plane_level, &xx[1], &yy[1]);
		GMT_geoz_to_xy (header.x_max - delx, header.y_max - dely, plane_level, &xx[2], &yy[2]);
		GMT_geoz_to_xy (header.x_min + delx, header.y_max - dely, plane_level, &xx[3], &yy[3]);
		if (!GMT_is_fnan (topo[nw])) {
			GMT_geoz_to_xy (header.x_min + delx, header.y_max - dely, (double)(topo[nw]), &x_left, &y_top);
			ps_plot (x_left, y_top, 3);	ps_plot (xx[3], yy[3], 2);
		}
		if (!GMT_is_fnan (topo[ne])) {
			GMT_geoz_to_xy (header.x_max - delx, header.y_max - dely, (double)(topo[ne]), &x_right, &y_top);
			ps_plot (x_right, y_top, 3);	ps_plot (xx[2], yy[2], 2);
		}
		if (!GMT_is_fnan (topo[se])) {
			GMT_geoz_to_xy (header.x_max - delx, header.y_min + dely, (double)(topo[se]), &x_right, &y_bottom);
			ps_plot (x_right, y_bottom, 3);	ps_plot (xx[1], yy[1], 2);
		}
		if (!GMT_is_fnan (topo[sw])) {
			GMT_geoz_to_xy (header.x_min + delx, header.y_min + dely, (double)(topo[sw]), &x_left, &y_bottom);
			ps_plot (x_left, y_bottom, 3);	ps_plot (xx[0], yy[0], 2);
		}
		ps_plot (0.0, 0.0, 3);	/* To make sure path is stroked */
	}
	
	if (tiling) {	/* Plot image as polygonal pieces */
		int p;
		double *xx, *yy, d, xpos, ypos;
		
		d = (off) ? 0.0 : 0.5;
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Tiling without interpolation\n", GMT_program);
		
		for (j = k = 0; j < header.ny; j++) {
			for (i = 0; i < header.nx; i++, k++) {	/* Compute rgb for each pixel */
				if (GMT_is_fnan (topo[k]) && tiling == 2) continue;
				GMT_get_rgb24 (topo[k], rgb);
				if (intens) GMT_illuminate (intensity[k], rgb);
				
				n = GMT_graticule_path (&xx, &yy, 1, xval[i] - dx2, xval[i] + dx2, yval[j] - dy2, yval[j] + dy2);
				
				/* Do map projection and plot */
				
				for (p = 0; p < n; p++) {
					GMT_geoz_to_xy (xx[p], yy[p], topo[k], &xpos, &ypos);
					xx[p] = xpos;	yy[p] = ypos;
				}
				ps_polygon (xx, yy, n, rgb, FALSE);
				GMT_free ((void *)xx);
				GMT_free ((void *)yy);
			}
		}
	}
	if (image) {	/* compute image */
		int nx_i, ny_i, kk, ip, jp, min_i, max_i, min_j, max_j, dist, node, nm_i, layers, *ix, *iy;
		int *top_jp, *bottom_jp, p;
		BOOLEAN done;
		double xp, yp, sum_w, w, sum_i, x_width, y_width;
		double sum_r, sum_g, sum_b, intval;
		unsigned char *bitimage_24, *bitimage_8;
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: get and store projected vertices\n", GMT_program);
		ps_comment ("Plot 3-D surface using scanline conversion of polygons to raster image");

		x_width = z_project.xmax - z_project.xmin;	/* Size of image in inches */
		y_width = z_project.ymax - z_project.ymin;
		nx_i = irint (x_width * dpi_i);	/* Size of image in pixels */
		ny_i = irint (y_width * dpi_i);
		
		ix = (int *) GMT_memory (VNULL, (size_t)nm, sizeof (int), GMT_program);
		iy = (int *) GMT_memory (VNULL, (size_t)nm, sizeof (int), GMT_program);
		
		for (j = bin = ij = 0; j < header.ny; j++) {	/* Get projected coordinates converted to pixel locations */
			ij = (two) ? (j + 2) * mx + two : bin;
			for (i = 0; i < header.nx; i++, bin++, ij++) {
				if (GMT_is_fnan (topo[ij])) {	/* Outside -R or NaNs not used */
					ix[bin] = iy[bin] = -1;
				}
				else {
					GMT_geoz_to_xy (xval[i], yval[j], (double)topo[ij], &xp, &yp);
					ix[bin] = (int)floor ((xp - z_project.xmin) * dpi_i);
					iy[bin] = (int)floor ((yp - z_project.ymin) * dpi_i);
				}
			}
		}
		GMT_free ((void *) xval);
		GMT_free ((void *) yval);
		
		/* Allocate image array and set background to PAGE_COLOR */
		
		if (monochrome) {
			int gray;
			
			nm_i = nx_i * ny_i;
			layers = 1;
			bitimage_8 = (unsigned char *) GMT_memory (VNULL, (size_t)nm_i, sizeof (char), GMT_program);
			gray = YIQ (gmtdefs.page_rgb);
			memset ((void *)bitimage_8, gray, (size_t)nm_i);
		}
		else {
			nm_i = nx_i * ny_i * 3;
			layers = 3;
			bitimage_24 = (unsigned char *) GMT_memory (VNULL, (size_t)nm_i, sizeof (char), GMT_program);
			kk = 0;
			while (kk < nm_i) {
				bitimage_24[kk++] = (unsigned char)gmtdefs.page_rgb[0];
				bitimage_24[kk++] = (unsigned char)gmtdefs.page_rgb[1];
				bitimage_24[kk++] = (unsigned char)gmtdefs.page_rgb[2];
			}
		}

		/* Set up arrays for staircase clippath and initialize them */
		top_jp = (int *) GMT_memory (VNULL, (size_t)nx_i, sizeof (int), GMT_program);
		bottom_jp = (int *) GMT_memory (VNULL, (size_t)nx_i, sizeof (int), GMT_program);
		for (ip = 0; ip < nx_i; ip++) bottom_jp[ip] = ny_i;
			
		/* Plot from back to front */
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Start rasterization\n", GMT_program);
		for (j = j_start; j != j_stop; j += j_inc) {

			if (gmtdefs.verbose) fprintf (stderr, "%s: Scan line conversion at j-line %.4d\r", GMT_program, j);

			for (i = i_start; i != i_stop; i += i_inc) {
				bin = j * header.nx + i;
				ij = (two) ? (j + 2) * header.nx + i + 2 : bin;
				for (k = bad = 0; !bad && k < 4; k++) bad = (ix[bin+bin_inc[k]] < 0 || iy[bin+bin_inc[k]] < 0);
				if (bad) continue;

				min_i = max_i = ix[bin];
				min_j = max_j = iy[bin];
				for (k = 1; k < 4; k++) {
					p = bin+bin_inc[k];
					if (ix[p] < min_i) min_i = ix[p];
					if (ix[p] > max_i) max_i = ix[p];
					if (iy[p] < min_j) min_j = iy[p];
					if (iy[p] > max_j) max_j = iy[p];
				}
				for (jp = min_j; jp <= max_j; jp++) {
					if (jp < 0 || jp >= ny_i) continue;
					for (ip = min_i; ip <= max_i; ip++) {
						if (ip < 0 || ip >= nx_i) continue;
						if (!pixel_inside (ip, jp, ix, iy, bin)) continue;
						
						/* Update clip mask */
						if (jp > top_jp[ip]) top_jp[ip] = jp; 
						if (jp < bottom_jp[ip]) bottom_jp[ip] = jp;
						
						kk = layers * ((ny_i-jp-1) * nx_i + ip);
						sum_r = sum_g = sum_b = sum_w = sum_i = 0.0;
						done = FALSE;
						for (k = 0; !done && k < 4; k++) {
							node = bin + bin_inc[k];
							GMT_get_rgb24 (grd[node], rgb);
							dist = quick_idist (ip, jp, ix[node], iy[node]);
							if (dist == 0) {	/* Only need this corner value */
								done = TRUE;
								if (intens) intval = intensity[node];
							}
							else {
								w = 1.0 / (double)dist;
								sum_r += rgb[0] * w;
								sum_g += rgb[1] * w;
								sum_b += rgb[2] * w;
								if (intens) sum_i += intensity[node] * w;
								sum_w += w;
							}
						}
						if (!done) {	/* Must get weighted value */
							sum_w = 1.0 / sum_w;
							rgb[0] = irint (sum_r * sum_w);
							rgb[1] = irint (sum_g * sum_w);
							rgb[2] = irint (sum_b * sum_w);
							if (intens) intval = sum_i * sum_w;
						}
						if (intens) GMT_illuminate (intval, rgb);
						if (monochrome) /* YIQ transformation */
							bitimage_8[kk] = (unsigned char) YIQ (rgb);
						else {
							bitimage_24[kk++] = (unsigned char) rgb[0];
							bitimage_24[kk++] = (unsigned char) rgb[1];
							bitimage_24[kk] = (unsigned char) rgb[2];
						}
					}
				}
			}
		}
		if (gmtdefs.verbose) fprintf (stderr, "\n");
		
		/* Must implement the clip path for the image */
			
		x_pixel_size = x_width / (double)nx_i;
		y_pixel_size = y_width / (double)ny_i;
		n4 = 4 * nx_i;
		x_imask = (double *) GMT_memory (VNULL, (size_t)n4, sizeof (double), GMT_program);
		y_imask = (double *) GMT_memory (VNULL, (size_t)n4, sizeof (double), GMT_program);
		nk = n4 - 1;
		
		for (ip = k = 0; ip < nx_i; ip++, k+= 2) {
			k1 = k + 1;
			x_imask[k]  = x_imask[nk-k]  = z_project.xmin + ip * x_pixel_size;
			x_imask[k1] = x_imask[nk-k1] = x_imask[k] + x_pixel_size;
			if (top_jp[ip] < bottom_jp[ip]) {	/* No pixels set in this column */
				y_imask[k] = y_imask[k1] = y_imask[nk-k] = y_imask[nk-k1] = z_project.ymin;
			}
			else {	/* Set top of upper pixel and bottom of lower pixel */
				y_imask[k] = y_imask[k1] = z_project.ymin + (top_jp[ip] + 1) * y_pixel_size;
				y_imask[nk-k] = y_imask[nk-k1] = z_project.ymin + bottom_jp[ip] * y_pixel_size;
			}
		}
		ps_clipon (x_imask, y_imask, 4 * nx_i, GMT_no_rgb, 3);
		GMT_free ((void *)x_imask);
		GMT_free ((void *)y_imask);
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Creating PostScript image ", GMT_program);
		if (monochrome) {
			if (gmtdefs.verbose) fprintf (stderr, "[B/W image]\n");
			ps_image (z_project.xmin, z_project.ymin, x_width, y_width, bitimage_8, nx_i, ny_i, 8);
			GMT_free ((void *)bitimage_8);
		}
		else {
			if (gmtdefs.verbose) fprintf (stderr, "[%s]\n", c_method[gmtdefs.color_image]);
			GMT_color_image (z_project.xmin, z_project.ymin, x_width, y_width, bitimage_24, nx_i, ny_i);
			GMT_free ((void *)bitimage_24);
		}
		ps_clipoff ();
		GMT_free ((void *)top_jp);
		GMT_free ((void *)bottom_jp);
		
		GMT_free ((void *)ix);
		GMT_free ((void *)iy);
	}

	if (mesh) {
		ps_comment ("Start of mesh plot");
		GMT_setpen (&pen[1]);
		for (j = j_start; j != j_stop; j += j_inc) {
			y_bottom = yval[j];
			y_top = y_bottom + abs (j_inc) * header.y_inc;
			for (i = i_start; i != i_stop; i += i_inc) {
				bin = j * header.nx + i;
				ij = (two) ? (j + 2) * mx + i + 2 : bin;
				for (k = bad = 0; !bad && k < 4; k++) bad += GMT_is_fnan (topo[ij+ij_inc[k]]);
				if (bad) continue;
				x_left = xval[i];
				x_right = x_left + abs (i_inc) * header.x_inc;
				GMT_geoz_to_xy (x_left, y_bottom, (double)(topo[ij+ij_inc[0]]), &xx[0], &yy[0]);
				GMT_geoz_to_xy (x_right, y_bottom, (double)(topo[ij+ij_inc[1]]), &xx[1], &yy[1]);
				GMT_geoz_to_xy (x_right, y_top, (double)(topo[ij+ij_inc[2]]), &xx[2], &yy[2]);
				GMT_geoz_to_xy (x_left, y_top, (double)(topo[ij+ij_inc[3]]), &xx[3], &yy[3]);
				ps_patch (xx, yy, 4, fill.rgb, TRUE);
				if (draw_contours) {
					pen_not_set = TRUE;
					if (binij[bin].first_cont == NULL) continue;
					for (this_cont = binij[bin].first_cont->next_cont; this_cont; this_cont = this_cont->next_cont) {
						for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) {
							z_val = (drape) ? GMT_get_bcr_z (&t_head, (double)this_point->x, (double)this_point->y, topo, &edgeinfo) : this_cont->value;
							if (GMT_is_dnan (z_val)) continue;
							GMT_geoz_to_xy ((double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]);
							k++;
						}
						if (pen_not_set) {
							GMT_setpen (&pen[0]);
							pen_not_set = FALSE;
						}
						ps_line (xx, yy, k, 3, FALSE, TRUE);
					}
					if (!pen_not_set) GMT_setpen (&pen[1]);
				}					
			}
		}
		GMT_free ((void *) xval);
		GMT_free ((void *) yval);
	}
	else if (surface) {
		ps_comment ("Start of filled surface");
		if (outline) GMT_setpen (&pen[1]);
		
		for (j = j_start; j != j_stop; j += j_inc) {
			y_bottom = yval[j];
			y_top = y_bottom + header.y_inc;
			for (i = i_start; i != i_stop; i += i_inc) {
				bin = j * header.nx + i;
				ij = (two) ? (j + 2) * mx + i + 2 : bin;
				for (k = bad = 0; !bad && k < 4; k++) bad += GMT_is_fnan (topo[ij+ij_inc[k]]);
				if (bad) continue;
				x_left = xval[i];
				x_right = x_left + header.x_inc;
				if (intens) {
					this_intensity = get_intensity (intensity, bin, header.nx);
					if (GMT_is_dnan (this_intensity)) continue;

				}
				/* Get mesh polygon */

				GMT_geoz_to_xy (x_left, y_bottom, (double)(topo[ij+ij_inc[0]]), &xmesh[0], &ymesh[0]);
				GMT_geoz_to_xy (x_right, y_bottom, (double)(topo[ij+ij_inc[1]]), &xmesh[1], &ymesh[1]);
				GMT_geoz_to_xy (x_right, y_top, (double)(topo[ij+ij_inc[2]]), &xmesh[2], &ymesh[2]);
				GMT_geoz_to_xy (x_left, y_top, (double)(topo[ij+ij_inc[3]]), &xmesh[3], &ymesh[3]);

				if (get_contours && binij[bin].first_cont) {	/* Contours go thru here */
					start_cont = binij[bin].first_cont->next_cont;
					
					/* Paint the part of the tile that has a value below the contour value */
					
					/* Set the right color */
					
					if (GMT_continuous) {
						z_ave = 0.0;
						for (n = 0; n < 4; n++) z_ave += MIN (grd[bin+bin_inc[n]], start_cont->value);
						z_ave *= 0.25;
					}
					else
						z_ave = start_cont->value - small;
		
					index = GMT_get_rgb24 (z_ave, rgb);
					if (index >= 0 && GMT_lut[index].skip) continue;
					if (intens) GMT_illuminate (this_intensity, rgb);
					ps_patch (xmesh, ymesh, 4, rgb, FALSE);
					
					/* Now loop over all the contours that crossed this tile */
					
					saddle = FALSE;
					cross_ok = TRUE;
					n_saddle = 0;
					
					for (this_cont = start_cont; this_cont; this_cont = this_cont->next_cont) {
					
						next_up = (this_cont->next_cont) ? this_cont->next_cont->value : DBL_MAX;
						/* First get all the x/y pairs for this contour */
						
						for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) {
							x[k] = this_point->x;
							y[k] = this_point->y;
							z[k] = (drape) ? GMT_get_bcr_z (&t_head, x[k], y[k], topo, &edgeinfo) : this_cont->value;
							if (GMT_is_dnan (z[k])) continue;
							v[k] = this_cont->value;
							k++;
						}
						
						/* Check for saddle situation */
						
						if (!saddle && this_cont->next_cont && this_cont->next_cont->value == this_cont->value) saddle = TRUE;
						if (saddle) n_saddle++;
						
						/* First make sure no grd-value equals the contour value */
						
						
						for (side = 0; side < 4; side++) {
							if (grd[bin+bin_inc[side]] == this_cont->value) {
								was[side] = grd[bin+bin_inc[side]];
								grd[bin+bin_inc[side]] += (float)small;
								moved[side] = TRUE;
							}
							else
								moved[side] = FALSE;
						}
						
						/* Then figure out on what side the contour exits */
						
						del_x = x[k-1] - x_left;
						if (del_x > dx2) del_x = header.x_inc - del_x;
						del_y = y[k-1] - y_bottom;
						if (del_y > dy2) del_y = header.y_inc - del_y;
						if (del_x < del_y) /* Cutting N-S gridlines */
							side = ((x[k-1]-x_left) > dx2) ? 1 : 3;
						else /* Cutting E-W gridlines */
							side = ((y[k-1]-y_bottom) > dy2) ? 2 : 0;
							
						/* If saddle, figure out on what side the contour enters */
						
						if (saddle && n_saddle == 1) {
							del_x = x[0] - x_left;
							if (del_x > dx2) del_x = header.x_inc - del_x;
							del_y = y[0] - y_bottom;
							if (del_y > dy2) del_y = header.y_inc - del_y;
							if (del_x < del_y) /* Cutting N-S gridlines */
								entry = ((x[0]-x_left) > dx2) ? 1 : 3;
							else /* Cutting E-W gridlines */
								entry = ((y[0]-y_bottom) > dy2) ? 2 : 0;
							
							if ((entry == 0 && side == 3) || (entry == 3 && side == 0))
								check = 0;
							else if ((entry == 1 && side == 2) || (entry == 2 && side == 1))
								check = 2;
							else if ((entry == 0 && side == 1) || (entry == 1 && side == 0))
								check = 1;
							else
								check = 3;
							
							cross_ok = (grd[bin+bin_inc[check]] < this_cont->value);
								
						}
							
						/* Now determine which way we have to go around the tile
						   outline to close the polygon */
						   
						way = (grd[bin+bin_inc[side]] >= this_cont->value) ? -1 : 1;
						
						/* Add the corner points that are needed to close the polygon */
						
						prev_side = side;
						if (way == 1) side++;
						for (n = n_added = 0; n < 4; n++) {
							side %= 4;
							if (cross_ok)
								go = TRUE;
							else if (n_added == 0)
								go = TRUE;
							else
								go = (side != (prev_side + 2)%4);
							if (go && grd[bin+bin_inc[side]] > this_cont->value) {
								x[k] = x_left + x_inc[side];
								y[k] = y_bottom + y_inc[side];
								z[k] = topo[ij+ij_inc[side]];
								v[k] = grd[bin+bin_inc[side]];
								k++;
								n_added++;
								prev_side = side;
							}
							side += way;
							if (side < 0) side += 4;
						}
						
						/* Compute the xy from the xyz triplets */
						
						for (n = 0; n < k; n++) GMT_geoz_to_xy (x[n], y[n], z[n], &xx[n], &yy[n]);
						
						if (GMT_continuous) {
							z_ave = 0.0;
							for (n = 0; n < k; n++) z_ave += MIN (v[n], next_up);
							z_ave /= k;
						}
						else
							z_ave = this_cont->value;
						
						/* Set the right color */
						
						GMT_get_rgb24 (z_ave+small, rgb);
						if (intens) GMT_illuminate (this_intensity, rgb);

						if (k > 2) ps_patch (xx, yy, k, rgb, FALSE);	/* Contours drawn separately below if desired */
						
						/* Reset nodes that was moved off contour-value */
						
						for (side = 0; side < 4; side++) if (moved[side]) grd[bin+bin_inc[side]] = was[side];
						
						if (n_saddle == 2) {
							n_saddle = 0;
							saddle = FALSE;
							cross_ok = TRUE;
						}
					}
					
					/* Draw contour lines if desired */
					
					pen_not_set = TRUE;
					for (this_cont = start_cont; draw_contours && this_cont; this_cont = this_cont->next_cont) {
						for (k = 0, this_point = this_cont->first_point; this_point; this_point = this_point->next_point) {
							z_val = (drape) ? GMT_get_bcr_z (&t_head, (double)this_point->x, (double)this_point->y, topo, &edgeinfo) : this_cont->value;
							if (GMT_is_dnan (z_val)) continue;

							GMT_geoz_to_xy ((double)this_point->x, (double)this_point->y, z_val, &xx[k], &yy[k]);
							k++;
						}
						if (pen_not_set) {
							GMT_setpen (&pen[0]);
							pen_not_set = FALSE;
						}
						ps_line (xx, yy, k, 3, FALSE, TRUE);
					}
					if (!pen_not_set) GMT_setpen (&pen[1]);
					if (outline) ps_patch (xmesh, ymesh, 4, GMT_no_rgb, TRUE);
				}
				else {	/* No Contours, paint tile (any of the corner values will give correct color */
					
					if (GMT_continuous) {
						for (n = 0, z_ave = 0.0; n < 4; n++) z_ave += grd[bin+bin_inc[n]];
						z_ave *= 0.25;
					}
					else
						z_ave = grd[bin];
					index = GMT_get_rgb24 (z_ave, rgb);
					if (index >= 0 && GMT_lut[index].skip) continue;
					if (intens) GMT_illuminate (this_intensity, rgb);
					ps_patch (xmesh, ymesh, 4, rgb, outline);
				}
			}
		}
		GMT_free ((void *) xval);
		GMT_free ((void *) yval);
	}			
	
	if (pen[1].texture || pen[0].texture) ps_setdash (CNULL, 0);

	if (project_info.z_pars[0] == 0.0) GMT_map_clip_off();

	if (facade) {	/* Cover the two front sides */
		ps_comment ("Paiting the frontal facade");
		GMT_setpen (&pen[0]);
		if (!z_project.draw[0])	{	/* Southern side */
			for (i = n = 0, ij = sw; i < header.nx; i++, ij++) {
				if (GMT_is_fnan (topo[ij])) continue;
				GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_min+dely, (double)(topo[ij]), &xx[n], &yy[n]);
				n++;
			}
			for (i = header.nx - 1; i >= 0; i--, n++) GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_min+dely, plane_level, &xx[n], &yy[n]);
			ps_polygon (xx, yy, n, rgb_facade, TRUE);
		}
		if (!z_project.draw[1]) {	/*	Eastern side */
			for (j = n = 0, ij = ne; j < header.ny; j++, ij += mx) {
				if (GMT_is_fnan (topo[ij])) continue;
				GMT_geoz_to_xy (header.x_max-delx, header.y_max-j*header.y_inc-dely, (double)(topo[ij]), &xx[n], &yy[n]);
				n++;
			}
			for (j = header.ny - 1; j >= 0; j--, n++) GMT_geoz_to_xy (header.x_max-delx, header.y_max-j*header.y_inc-dely, plane_level, &xx[n], &yy[n]);
			ps_polygon (xx, yy, n, rgb_facade, TRUE);
		}
		if (!z_project.draw[2])	{	/* Northern side */
			for (i = n = 0, ij = nw; i < header.nx; i++, ij++) {
				if (GMT_is_fnan (topo[ij])) continue;
				GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_max-dely, (double)(topo[i]), &xx[n], &yy[n]);
				n++;
			}
			for (i = header.nx - 1; i >= 0; i--, n++) GMT_geoz_to_xy (header.x_min+i*header.x_inc+delx, header.y_max-dely, plane_level, &xx[n], &yy[n]);
			ps_polygon (xx, yy, n, rgb_facade, TRUE);
		}
		if (!z_project.draw[3]) {	/*	Western side */
			for (j = n = 0, ij = nw; j < header.ny; j++, ij += mx) {
				if (GMT_is_fnan (topo[ij])) continue;
				GMT_geoz_to_xy (header.x_min+delx, header.y_max-j*header.y_inc-dely, (double)(topo[ij]), &xx[n], &yy[n]);
				n++;
			}
			for (j = header.ny - 1; j >= 0; j--, n++) GMT_geoz_to_xy (header.x_min+delx, header.y_max-j*header.y_inc-dely, plane_level, &xx[n], &yy[n]);
			ps_polygon (xx, yy, n, rgb_facade, TRUE);
		}
	}
	
	if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
	
	ps_plotend (gmtdefs.last_page);
	
	/* Free memory */
	
	if (get_contours) {
		for (ij = 0; ij < nm; ij++) {
			if (!binij[ij].first_cont) continue;
			last_cont = binij[ij].first_cont;
			for (this_cont = binij[ij].first_cont->next_cont; this_cont; this_cont = this_cont->next_cont) {
				if (this_cont->first_point) {
					last_point = this_cont->first_point;
					for (this_point = this_cont->first_point->next_point; this_point; this_point = this_point->next_point) {
						GMT_free ((void *)last_point);
						last_point = this_point;
					}
					GMT_free ((void *)last_point);
				}
				GMT_free ((void *)last_cont);
				last_cont = this_cont;
			}
			GMT_free ((void *)last_cont);
		}
		GMT_free ((void *)binij);
	}
	
	GMT_free ((void *)xx);
	GMT_free ((void *)yy);
	GMT_free ((void *)x);
	GMT_free ((void *)y);
	GMT_free ((void *)z);
	GMT_free ((void *)v);
	GMT_free ((void *)topo);
	if (intens) GMT_free ((void *)intensity);
	if (drape) GMT_free ((void *)grd);
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program);
	
	GMT_end (argc, argv);
}

struct CONT *get_cont_struct (int bin, double value)
{
	struct CONT *cont, *new_cont;
	
	if (!binij[bin].first_cont) binij[bin].first_cont = (struct CONT *) GMT_memory (VNULL, (size_t)1, sizeof (struct CONT), GMT_program);

	for (cont = binij[bin].first_cont; cont->next_cont && cont->next_cont->value <= value; cont = cont->next_cont);
	
	new_cont = (struct CONT *) GMT_memory (VNULL, (size_t)1, sizeof (struct CONT), GMT_program);
	if (cont->next_cont) {	/* Put it in the link */
		new_cont->next_cont = cont->next_cont;
		cont->next_cont = new_cont;
	}
	else	/* End of list */
		cont->next_cont = new_cont;
	return (new_cont);
}

struct POINT *get_point (double x, double y)
{
	struct POINT *point;
	
	point = (struct POINT *) GMT_memory (VNULL, (size_t)1, sizeof (struct POINT), GMT_program);
	point->x = x;
	point->y = y;
	return (point);
}

void grdview_init_setup (struct GRD_HEADER *header, float *topo, int two, BOOLEAN draw_plane, double plane_level)
{
	int i, j, ij;
	double xtmp, ytmp, tmp, delx, dely;
	
	delx = (header->node_offset) ? 0.5 * header->x_inc :0.0;
	dely = (header->node_offset) ? 0.5 * header->y_inc :0.0;
	/* Find projected min/max in y-direction */
	
	z_project.ymax = z_project.ymin;	/* Reset from whatever it was */
	
	for (j = 0; j < header->ny; j++) {
		ij = (j + two) * header->nx + two;
		for (i = 0; i < header->nx; i++, ij++) {
			if (GMT_is_fnan (topo[ij])) continue;
			GMT_geoz_to_xy (header->x_min + i * header->x_inc, header->y_max - j * header->y_inc, (double)topo[ij], &xtmp, &ytmp);
			z_project.ymin = MIN (z_project.ymin, ytmp);
			z_project.ymax = MAX (z_project.ymax, ytmp);
		}
	}
	if (draw_plane) {	/* plane or facade may exceed the found min/max */
		for (i = 0; i < header->nx; i++) {
			tmp = header->x_min + i * header->x_inc + delx;
			GMT_geoz_to_xy (tmp, header->y_min + dely, plane_level, &xtmp, &ytmp);
			z_project.ymin = MIN (z_project.ymin, ytmp);
			z_project.ymax = MAX (z_project.ymax, ytmp);
			GMT_geoz_to_xy (tmp, header->y_max-dely, plane_level, &xtmp, &ytmp);
			z_project.ymin = MIN (z_project.ymin, ytmp);
			z_project.ymax = MAX (z_project.ymax, ytmp);
		}
		for (j = 0; j < header->ny; j++) {
			tmp = header->y_max - j * header->y_inc - dely;
			GMT_geoz_to_xy (header->x_min+delx, tmp, plane_level, &xtmp, &ytmp);
			z_project.ymin = MIN (z_project.ymin, ytmp);
			z_project.ymax = MAX (z_project.ymax, ytmp);
			GMT_geoz_to_xy (header->x_max-delx, tmp, plane_level, &xtmp, &ytmp);
			z_project.ymin = MIN (z_project.ymin, ytmp);
			z_project.ymax = MAX (z_project.ymax, ytmp);
		}
	}
}

double get_intensity (float *intensity, int k, int nx)
{
	/* Finds the agerage intensity for this polygon */
	return (0.25 * (intensity[k] + intensity[k+1] + intensity[k-nx] + intensity[k-nx+1]));
}

int pixel_inside (int ip, int jp, int *ix, int *iy, int bin)
{
	int i, what;
	double x[6], y[6];
	
	for (i = 0; i < 4; i++) {
		x[i] = ix[bin+bin_inc[i]];
		y[i] = iy[bin+bin_inc[i]];
	}
	x[4] = x[0];	y[4] = y[0];
	what = GMT_non_zero_winding ((double)ip, (double)jp, x, y, 5);
	return (what);
}

int quick_idist (int x1, int y1, int x2, int y2)
{
	if ((x2 -= x1) < 0) x2 = -x2;
	if ((y2 -= y1) < 0) y2 = -y2;
	return (x2 + y2 - (((x2 > y2) ? y2 : x2) >> 1));
}