/*--------------------------------------------------------------------
 *	$Id: grdcontour.c,v 1.57 2004/09/10 23:16:15 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 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
 *--------------------------------------------------------------------*/
/*
 * grdcontour reads a 2-D grd file and contours it. This algorithm handles
 * cases where the old contourm failed, e.g. when a contour line passes
 * exactly through a data point. This is achieved by adding a tiny
 * number to those data points that would have a contour line going
 * through them. This is cheating, but the human eye cannot tell the
 * difference, so who cares?
 * A multitude of options exist and are explained in the usage message.
 * The 2-D grd file format is outlined in the routines read_grd/write_grd.
 *
 * Author:	Paul Wessel
 * Date:	6-JAN-1991
 * Version:	2.0	Based on old v1.x by Paul Wessel & Walter Smith
 * Revised:	21-JUL-1998 for GMT 3.1
 * Revised:	03-MAR-1999 for GMT 3.2
 * Revised:	25-FEB-2000 for GMT 3.3.4
 * Revised:	25-JUN-2000 for GMT 3.3.5
 * Version:	3.4 Added Kaj Jancke's modification for transparent label boxes
 * Version:	4.0 15-JUL-2004 Implemented label clipping in PostScript. Expeanded -G option
 *
 */

#include "gmt.h"

#define MIN_LENGTH 0.01

float *grd;
int *edge;
double *x, *y;	/* Arrays holding the contour xy values */

struct SAVE {
	double *x, *y;
	double cval;
	int n;
	struct GMT_PEN pen;
	BOOLEAN do_it, high;
} *save;

void sort_and_plot_ticks (struct SAVE *save, int n, struct GRD_HEADER *h, float *grd, double tick_gap, double tick_length, BOOLEAN tick_low, BOOLEAN tick_high, BOOLEAN tick_label, char *labels);
void GMT_grd_minmax (float *a, struct GRD_HEADER *h, double xyz[2][3]);

main (int argc, char **argv)
{
	int i, j, k, n, nn, nm, c, n_edges, n_alloc;
	int section, n_contours, smooth_factor = 0, side, id, bad;
	int off, nx, ny, n_save = 0, got, ncut = 0, rgb[3];
	
	BOOLEAN error = FALSE, first = TRUE, dump = FALSE, interior, begin;
	BOOLEAN do_ticks = FALSE, use_cpt_colors = FALSE, use_cpt_annot = TRUE;
	BOOLEAN tick_high = FALSE, tick_low = FALSE, cpt_given = FALSE, tick_label = FALSE;
	
	char *grdfile = 0, dfile[BUFSIZ], line[BUFSIZ], *cont_type, *cont_do_tick;
	char *cpt_file, *labels, txt_a[32], txt_b[32], cont_label[128], format[128];
	
	double c_int = 0.0, c_low, c_high, aval, a_int = 0.0, tick_gap = 0.2;
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, tick_length = 0.04, *xp, *yp;
	double *contour, cval, mult = 1.0, shift = 0.0, min, max, small, label_angle = 0.0, xyz[2][3];
	double small_x, small_y, data_west, data_east, data_south, data_north, tmp, *cont_angle;
	
	float *grd_original;
	
	FILE *fpc = NULL;
	
	struct GRD_HEADER header;
	
	struct GMT_PEN pen[2];
	struct GMT_CONTOUR G;

	argc = GMT_begin (argc, argv);
	
	GMT_init_pen (&pen[0], GMT_PENWIDTH);
	GMT_init_pen (&pen[1], 3.0 * GMT_PENWIDTH);
	c_low = -DBL_MAX;
	c_high = DBL_MAX;
	dfile[0] = 0;
	GMT_contlabel_init (&G);
	if (gmtdefs.measure_unit == GMT_CM) {
		tick_gap = 0.5 / 2.54;
		tick_length = 0.1 / 2.54;
	}
	
	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 'b':
				case 'c':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
			
				case 'A':
					/* Format can be one of two:
					 * 3.4.x: -A[-|aint][f<fontsize>][a<angle>][/<r/g/b>][t|o]
					 * 4.x:   -A[-|aint][+a<angle>][+c<dx>[/<dy>]][+f<font>][+g<fill>][+j<just>][+l<label>][+o|O|t][+s<size>][+p<pen>][+u<unit>]
					 */

					bad = GMT_contlabel_specs (&argv[i][2], &G);
					if (argv[i][2] == '-')
						use_cpt_annot = FALSE;
					else {
						n = sscanf (&argv[i][2], "%lf", &a_int);
						G.annot = TRUE;
					}
					if (bad) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -A option.  Correct syntax:\n", GMT_program);
						fprintf (stderr, "\t-A[-][aint][+a<angle>][+c<dx>[/<dy>]][+f<font>][+g[<fill>]][+j<just>][+o][+p[<pen>]][+s<size>][+u<unit>][+v]\n");
						error += bad;
					}
					break;
				case 'C':
					if ((fpc = fopen (&argv[i][2], "r")) == NULL)
						c_int = atof (&argv[i][2]);
					else {
						c_int = 1.0;
						cpt_given = (int) strstr (&argv[i][2], ".cpt");
						cpt_file = &argv[i][2];
					}
					break;
				case 'D':
					dump = TRUE;
					strcpy (dfile, &argv[i][2]);
					break;
				case 'E':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'G':
					error += GMT_contlabel_info ('G', &argv[i][2], &G);
					break;
				case 'L':
					sscanf (&argv[i][2], "%lf/%lf", &c_low, &c_high);
					break;
				case 'M':	/* with -D, create one multiple line segments */
					GMT_multisegment (&argv[i][2]);
					GMT_io.multi_segments = 2;
					break;
				case 'N':	/* Backwards compatibility - now done in -A instead */
					if (argv[i][2])
						strcpy (G.unit, &argv[i][2]);
					else
						strcpy (G.unit, "z_units");
					break;
				case 'S':
					smooth_factor = atoi (&argv[i][2]);
					break;
				case 'Q':
					ncut = atoi (&argv[i][2]);
					break;
				case 'T':
					do_ticks = tick_high = tick_low = TRUE;	/* Default if just -T is given */
					if (argv[i][2]) {	/* But here we gave more options */
						if (argv[i][2] == '+')			/* Only tick local highs */
							tick_low = FALSE, j = 1;
						else if (argv[i][2] == '-')		/* Only tick local lows */
							tick_high = FALSE, j = 1;
						else
							j = 0;
						n = 0;
						if (strchr (&argv[i][2+j], '/')) {	/* Gave gap/length */
							n = sscanf (&argv[i][2+j], "%[^/]/%[^:]", txt_a, txt_b);
							if (n == 2) {
								tick_gap = GMT_convert_units (txt_a, GMT_INCH);
								tick_length = GMT_convert_units (txt_b, GMT_INCH);
							}
						}
						for (j = 2; argv[i][j] && argv[i][j] != ':'; j++);
						if (argv[i][j]) {	/* Gave high/low markers */
							j++;
							tick_label = TRUE;
							labels = &argv[i][j];
						}
						if (n == 1 || tick_gap <= 0.0 || tick_length == 0.0) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -T option.  Correct syntax:\n", GMT_program);
							fprintf (stderr, "\t-T[+|-][<tick_gap>[c|i|m|p]/<tick_length>[c|i|m|p]][:LH], <tick_gap> must be > 0\n");
							error = TRUE;
						}
					}
					break;
				case 'W':
					k = 2;
					if (argv[i][k] == '+') use_cpt_colors = TRUE, k++;
					j = (argv[i][k] == 'a' || argv[i][k] == 'c') ? k+1 : k;
					if (j == k) {	/* Set both */
						if (GMT_getpen (&argv[i][j], &pen[0])) {
							GMT_pen_syntax ('W');
							error++;
						}
						else pen[1] = pen[0];
					}
					else {
						id = (argv[i][k] == 'a') ? 1 : 0;
						if (GMT_getpen (&argv[i][j], &pen[id])) {
							GMT_pen_syntax ('W');
							error++;
						}
					}
					break;
				case 'Z':
					if (argv[i][2] && argv[i][2] != 'p') n = sscanf (&argv[i][2], "%lf/%lf", &mult, &shift);
					GMT_z_periodic = (argv[i][strlen(argv[i])-1] == 'p');	/* Phase data */
					break;
				default:
					error++;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			grdfile = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdcontour %s - Contouring of 2-D gridded data sets\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdcontour <grdfile> -C<cont_int> -J<params>\n");
		fprintf (stderr, "\t[-A[-|<annot_int>][<labelinfo>] [-B<tickinfo>] [-D<dumpfile>] [-Eaz/el] [-G[d|f|l|n|x]<info>]\n");
		fprintf (stderr, "\t[-K] [-L<Low/high>] [-M[<flag>]] [-O] [-P] [-Q<cut>] [-R<west/east/south/north>] [-S<smooth>]\n");
		fprintf (stderr, "\t[-T[+|-][<gap>[c|i|m|p]/<length>[c|i|m|p]][:LH]] [-U[<label>]] [-V] [-W[+]<type><pen>]\n");
		fprintf (stderr, "\t[-X<x_shift>] [-Y<y_shift>] [-Z[<fact>[/shift>]][p]] [-bo[s][<n>]] [-c<ncopies>]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<grdfile> is 2-D netCDF grdfile to be contoured\n");
		fprintf (stderr, "\t-C Contours to be drawn can be specified in one of three ways:\n");
		fprintf (stderr, "\t   1. Fixed contour interval\n");
		fprintf (stderr, "\t   2. Name of file with contour levels in col 1 and C(ont) or A(not) in col 2\n");
		fprintf (stderr, "\t      [and optionally an individual annotation angle in col 3.]\n");
		fprintf (stderr, "\t   3. Name of cpt-file\n");
		fprintf (stderr, "\t   If -T is used, only contours with upper case C or A is ticked\n");
		fprintf (stderr, "\t     [cpt-file contours are set to C unless last column has flags; Use -A to force all to A]\n");
		GMT_explain_option ('j');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A Annotation label information. [Default is no annoted contours].\n");
		fprintf (stderr, "\t   Give annotation interval OR - to disable all contour annotations implied in -C\n");
		fprintf (stderr, "\t   <labelinfo> controls the specifics of the labels.  Append what you need:\n");
		GMT_label_syntax (5, 0);
		GMT_explain_option ('b');
		fprintf (stderr, "\t-D to Dump contour lines to individual files (but see -M)\n");
		fprintf (stderr, "\t   Append file prefix [contour].  Files will be called <dumpfile>_<cont>_#[_i].xyz|b\n");
		fprintf (stderr, "\t   where <cont> is the contour value and # is a segment counter.\n");
		fprintf (stderr, "\t   _i is inserted for interior (closed) contours, with xyz (ascii) or b (binary) as extension.\n");
		fprintf (stderr, "\t   However, if -D- is given then files are C#_e or C#_i plus extension, where # is a running number.\n");
		fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
		fprintf (stderr, "\t-G Controls placement of labels along contours.  Choose among five algorithms:\n");
		GMT_cont_syntax (3, 0);
		GMT_explain_option ('K');
		fprintf (stderr, "\t-L only contour inside this range\n");
		fprintf (stderr, "\t-M Used with -D.   Create a single multiple segment file where contours are separated by a record\n");
		fprintf (stderr, "\t   whose first character is <flag> ['>'].  This header also has the contour level value\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr, "\t-Q Do not draw contours with less than <cut> points [Draw all contours]\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\t   [Default is extent of grdfile]\n");
		fprintf (stderr, "\t-S will Smooth contours by splining and resampling\n");
		fprintf (stderr, "\t   at approximately gridsize/<smooth> intervals\n");
		fprintf (stderr, "\t-T will embellish innermost, closed contours with ticks pointing in the downward direction\n");
		fprintf (stderr, "\t   User may specify to tick only highs (-T+) or lows (-T-) [-T means both]\n");
		fprintf (stderr, "\t   Append spacing/ticklength (apppend units) to change defaults [%g/%g %s]\n",
			tick_gap*GMT_u2u[GMT_INCH][gmtdefs.measure_unit], tick_length*GMT_u2u[GMT_INCH][gmtdefs.measure_unit], GMT_unit_names[gmtdefs.measure_unit]);
		fprintf (stderr, "\t   Append :LH to plot the characters L and H in the center of closed contours\n");
		fprintf (stderr, "\t   for local Lows and Highs (e.g, give :-+ to plot - and + signs)\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets pen attributes. <type> can be 'a' for annotated contours and\n");
		fprintf (stderr, "\t   'c' for regular contours.  The default settings are\n");
		fprintf (stderr, "\t   Contour pen:  width = %gp, color = (%d/%d/%d), texture = solid\n", pen[0].width, pen[0].rgb[0], pen[0].rgb[1], pen[0].rgb[2]);
		fprintf (stderr, "\t   Annotate pen: width = %gp, color = (%d/%d/%d), texture = solid\n", pen[1].width, pen[1].rgb[0], pen[1].rgb[1], pen[1].rgb[2]);
		fprintf (stderr, "\t   Use + to draw colored contours based on the cpt file\n");
		GMT_explain_option ('X');
		fprintf (stderr, "\t-Z to subtract <shift> and multiply data by <fact> before contouring [1/0].\n");
		fprintf (stderr, "\t   Append p for z-data that is periodic in 360 (i.e., phase data)\n");
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('c');
		exit (EXIT_FAILURE);
	}
	
	if (a_int > 0.0 && (!fpc && c_int == 0.0)) c_int = a_int;
	if (!fpc && c_int <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Must specify contour interval, file name with levels, or cpt-file\n", GMT_program);
		error++;
	}
	if (c_low >= c_high) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -L option: lower limit >= upper!\n", GMT_program);
		error++;
	}
	if (smooth_factor < 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Smooth_factor must be > 0\n", GMT_program);
		error++;
	}
	if (ncut < 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Q option:  Value must be >= 0\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 (!grdfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (G.label_dist_spacing <= 0.0 || G.half_width <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option.  Correct syntax:\n", GMT_program);
		fprintf (stderr, "\t-G<annot_dist>/<npoints>, both values must be > 0\n");
		error++;
	}
	if (mult == 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  factor must be nonzero\n", GMT_program);
		error++;
	}
	if (use_cpt_colors && !cpt_given) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -W option:  + only valid if -C sets a cpt file\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

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


	if (dump && dfile[0] == 0) strcpy (dfile,"contour");
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Allocate memory and read data file\n", GMT_program);
	
	if (!strcmp (grdfile,  "=")) {
		fprintf (stderr, "%s: Piping of grdfile not supported!\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	
	if (GMT_read_grd_info (grdfile, &header)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}
	off = (header.node_offset) ? 0 : 1;
	if (!strcmp (G.unit, "z_units")) strcpy (G.unit, header.z_units);
	
	/* Determine what wesn to pass to map_setup */

	if (!project_info.region_supplied) {
		west = header.x_min;
		east = header.x_max;
		south = header.y_min;
		north = header.y_max;
	}

	GMT_map_setup (west, east, south, north);

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

	if (GMT_grd_setregion (&header, &data_west, &data_east, &data_south, &data_north)) {	/* No grid to plot; just do empty map and exit */
		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, gmtdefs.encoding.name, GMT_epsinfo (argv[0]));
		GMT_echo_command (argc, argv);
		if (gmtdefs.unix_time) GMT_timestamp (argc, argv);
		GMT_map_basemap ();
		ps_plotend (gmtdefs.last_page);
		GMT_end (argc, argv);
	}
	

	/* Read data */

	nx = irint ( (data_east - data_west) / header.x_inc) + off;
	ny = irint ( (data_north - data_south) / header.y_inc) + off;
	nm = nx * ny;
	grd_original = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
	if (GMT_read_grd (grdfile, &header, grd_original, 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);
	}

	n_edges = header.ny * (int )ceil (header.nx / 16.0);
	edge = (int *) GMT_memory (VNULL, (size_t)n_edges, sizeof (int), GMT_program);
	
	if (!(mult == 1.0 && shift == 0.0)) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Subtracting %g and multiplying grid by %g\n", GMT_program, shift, mult);
		for (i = 0; i < nm; i++) grd_original[i] = (float)((grd_original[i] - shift) * mult);
		header.z_min = (header.z_min - shift) * mult;
		header.z_max = (header.z_max - shift) * mult;
		if (mult < 0.0) d_swap (header.z_min, header.z_max);
	}
	if (c_low > header.z_min) header.z_min = c_low;
	if (c_high < header.z_max) header.z_max = c_high;
	
	small = c_int * 1.0e-6;
	if (a_int == 0.0) a_int = c_int;
	
	small_x = 0.01 * header.x_inc;	small_y = 0.01 * header.y_inc;
	
	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, gmtdefs.encoding.name, GMT_epsinfo (argv[0]));
	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);

	GMT_map_clip_on (GMT_no_rgb, 3);
		
	if (G.annot) {	/* Wants annotated contours */
		aval = floor (header.z_min / a_int) * a_int;
		if (aval < header.z_min) aval += a_int;
	}
	else
		aval = header.z_max + 1.0;
	
	n_alloc = GMT_CHUNK;
	contour = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	cont_type = (char *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (char), GMT_program);
	cont_do_tick = (char *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (char), GMT_program);
	cont_angle = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	
	if (cpt_given) {	/* Presumably got a cpt-file */
		GMT_read_cpt (cpt_file);
		n_contours = GMT_n_colors + 1;
		if (n_contours > n_alloc) {
			contour = (double *) GMT_memory ((void *)contour, (size_t)n_contours, sizeof (double), GMT_program);
			cont_type = (char *) GMT_memory ((void *)cont_type, (size_t)n_contours, sizeof (char), GMT_program);
			cont_do_tick = (char *) GMT_memory ((void *)cont_do_tick, (size_t)n_contours, sizeof (char), GMT_program);
			cont_angle = (double *) GMT_memory ((void *)cont_angle, (size_t)n_contours, sizeof (double), GMT_program);
		}
		for (i = c = 0; i < GMT_n_colors; i++) {
			if (GMT_lut[i].skip) continue;
			contour[c] = GMT_lut[i].z_low;
			if (!use_cpt_annot)
				cont_type[c] = 'C';
			else if (GMT_lut[i].annot)
				cont_type[c] = 'A';
			else
				cont_type[c] = (G.annot) ? 'A' : 'C';
			cont_angle[c] = (G.angle_type == 2) ? label_angle : GMT_d_NaN;
			cont_do_tick[c] = do_ticks;
			c++;
		}
		contour[c] = GMT_lut[GMT_n_colors-1].z_high;
		if (!use_cpt_annot)
			cont_type[c] = 'C';
		else if (GMT_lut[GMT_n_colors-1].annot & 2)
			cont_type[c] = 'A';
		else
			cont_type[c] = (G.annot) ? 'A' : 'C';
		cont_angle[c] = (G.angle_type == 2) ? label_angle : GMT_d_NaN;
		cont_do_tick[c] = do_ticks;
		n_contours = c + 1;
	}
	else if (fpc != NULL) {	/* read contour info from file */
		n_contours = 0;
		while (fgets (line, BUFSIZ, fpc)) {
			if (line[0] == '#') continue;
			got = sscanf (line, "%lf %c %lf", &contour[n_contours], &cont_type[n_contours], &tmp);
			if (cont_type[n_contours] == 0) cont_type[n_contours] = 'C';
			cont_do_tick[n_contours] = (do_ticks && ((cont_type[n_contours] == 'C') || (cont_type[n_contours] == 'A'))) ? 1 : 0;
			cont_angle[n_contours] = (got == 3) ? tmp : GMT_d_NaN;
			n_contours++;
			if (n_contours == n_alloc) {
				n_alloc += GMT_CHUNK;
				contour = (double *) GMT_memory ((void *)contour, (size_t)n_alloc, sizeof (double), GMT_program);
				cont_type = (char *) GMT_memory ((void *)cont_type, (size_t)n_alloc, sizeof (char), GMT_program);
				cont_do_tick = (char *) GMT_memory ((void *)cont_do_tick, (size_t)n_alloc, sizeof (char), GMT_program);
				cont_angle = (double *) GMT_memory ((void *)cont_angle, (size_t)n_alloc, sizeof (double), GMT_program);
			}
		}
		fclose (fpc);
	}
	else {	/* Set up contour intervals automatically from c_int and a_int */
		min = floor (header.z_min / c_int) * c_int; if (!GMT_z_periodic && min < header.z_min) min += c_int;	
		max = ceil (header.z_max / c_int) * c_int; if (max > header.z_max) max -= c_int;
		for (c = irint (min/c_int), n_contours = 0; c <= irint (max/c_int); c++, n_contours++) {
			if (n_contours == n_alloc) {
				n_alloc += GMT_CHUNK;
				contour = (double *) GMT_memory ((void *)contour, (size_t)n_alloc, sizeof (double), GMT_program);
				cont_type = (char *) GMT_memory ((void *)cont_type, (size_t)n_alloc, sizeof (char), GMT_program);
				cont_angle = (double *) GMT_memory ((void *)cont_angle, (size_t)n_alloc, sizeof (double), GMT_program);
			}
			contour[n_contours] = c * c_int;
			if (G.annot && (contour[n_contours] - aval) > SMALL) aval += a_int;
			cont_type[n_contours] = (fabs (contour[n_contours] - aval) < SMALL) ? 'A' : 'C';
			cont_angle[n_contours] = (G.angle_type == 2) ? label_angle : GMT_d_NaN;
			cont_do_tick[n_contours] = do_ticks;
		}
	}

	if (GMT_z_periodic && fabs (contour[n_contours-1] - contour[0] - 360.0) < SMALL) {	/* Get rid of redundant contour */
		n_contours--;
	}

	contour = (double *) GMT_memory ((void *)contour, (size_t)n_contours, sizeof (double), GMT_program);
	cont_type = (char *) GMT_memory ((void *)cont_type, (size_t)n_contours, sizeof (char), GMT_program);
	cont_do_tick = (char *) GMT_memory ((void *)cont_do_tick, (size_t)n_contours, sizeof (char), GMT_program);
	cont_angle = (double *) GMT_memory ((void *)cont_angle, (size_t)n_contours, sizeof (double), GMT_program);
	
	if (do_ticks) save = (struct SAVE *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (struct SAVE), GMT_program);
	GMT_grd_minmax (grd_original, &header, xyz);
	if (GMT_contlabel_prep (&G, xyz, 3)) exit (EXIT_FAILURE);	/* Prep for crossing lines, if any */
	
	/* Because we are doing single-precision, we cannot subtract incrementally but must start with the
	 * original grid values and subtract the current contour value. */
	 
	grd = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);

	for (c = 0; c < n_contours; c++) {	/* For each contour value cval */
		
		/* Reset markers and set up new zero-contour*/
		
		cval = contour[c];
		if (gmtdefs.verbose) fprintf (stderr, "%s: Tracing the %g contour\n", GMT_program, cval);
		
		/* New approach to avoid round-off */
		
		for (i = 0; i < nm; i++) {
			grd[i] = grd_original[i] - (float)cval;		/* If there are NaNs they will remain NaNs */
			if (grd[i] == 0.0) grd[i] += (float)small;	/* There will be no actual zero-values, just -ve and +ve values */
		}
		
		section = 0;
		first = FALSE;
		id = (cont_type[c] == 'A' || cont_type[c] == 'a') ? 1 : 0;

		G.line_pen = pen[id];	/* Load current pen into contour structure */
		if (use_cpt_colors) {	/* Override pen & label color according to cpt file */
			GMT_get_rgb24 (cval, rgb);
			memcpy ((void *)&G.line_pen.rgb, (void *)rgb, (size_t)(3 * sizeof (int)));
			if (!G.got_font_rgb && G.curved_text) memcpy ((void *)&G.font_rgb, (void *)rgb, (size_t)(3 * sizeof (int)));
		}
		
		side = 0;
		begin = TRUE;
		while (side < 5) {
			while ((n = GMT_contours (grd, &header, smooth_factor, gmtdefs.interpolant, &side, edge, begin, &x, &y)) > 0) {
		
				if (fabs (x[0] - x[n-1]) < small_x && fabs (y[0] - y[n-1]) < small_y) {
					interior = TRUE;
					x[n-1] = x[0];	y[n-1] = y[0];	/* Ensure exact closure */
				}
				else
					interior = FALSE;
					
				if (n >= ncut) {
					if (dump) GMT_dump_contour (x, y, n, cval, section, interior, dfile);
					if ((nn = GMT_clip_to_map (x, y, n, &xp, &yp))) {	/* Lines inside the region */
						/* From here on, xp/yp are map inches */
						if (cont_type[c] == 'A' || cont_type[c] == 'a') {	/* Annotated contours */
							GMT_get_format (cval, G.unit, CNULL, format);
							sprintf (cont_label, format, cval);
						}
						if (cont_do_tick[c] && interior) {
							save[n_save].x = (double *) GMT_memory (VNULL, (size_t)nn, sizeof (double), GMT_program);
							save[n_save].y = (double *) GMT_memory (VNULL, (size_t)nn, sizeof (double), GMT_program);
							memcpy ((void *)save[n_save].x, (void *)xp, (size_t)(nn * sizeof (double)));
							memcpy ((void *)save[n_save].y, (void *)yp, (size_t)(nn * sizeof (double)));
							save[n_save].n = nn;
							memcpy ((void *)&save[n_save].pen, (void *)&pen[id], sizeof (struct GMT_PEN));
							save[n_save].do_it = TRUE;
							save[n_save].cval = cval;
							n_save++;
							if (n_save == n_alloc) {
								n_alloc += GMT_CHUNK;
								save = (struct SAVE *) GMT_memory ((void *)save, (size_t)n_alloc, sizeof (struct SAVE), GMT_program);
							}
						}
						GMT_hold_contour (&xp, &yp, nn, cont_label, cont_type[c], cont_angle[c], interior, &G);
						GMT_free ((void *)xp);
						GMT_free ((void *)yp);
					}
					section++;
				}
				begin = FALSE;
				GMT_free ((void *)x);
				GMT_free ((void *)y);
			}
			begin = FALSE;
		}
	}
	
	if (pen[0].texture || pen[1].texture) ps_setdash (CNULL, 0);

	if (do_ticks && n_save) {
		save = (struct SAVE *) GMT_memory ((void *)save, (size_t)n_save, sizeof (struct SAVE), GMT_program);
		sort_and_plot_ticks (save, n_save, &header, grd_original, tick_gap, tick_length, tick_low, tick_high, tick_label, labels);
		for (i = 0; i < n_save; i++) {
			GMT_free ((void *)save[i].x);
			GMT_free ((void *)save[i].y);
		}
		GMT_free ((void *)save);
	}
	
	GMT_contlabel_plot (&G);

	GMT_map_clip_off ();

	GMT_map_basemap ();
	
	ps_setpaint (gmtdefs.background_rgb);
	
	if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);

	ps_plotend (gmtdefs.last_page);
	
	GMT_free ((void *)grd);
	GMT_free ((void *)grd_original);
	GMT_free ((void *)edge);
	GMT_free ((void *)contour);
	GMT_free ((void *)cont_type);
	GMT_free ((void *)cont_angle);
	GMT_free ((void *)cont_do_tick);

	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program);
	
	GMT_end (argc, argv);
}

void sort_and_plot_ticks (struct SAVE *save, int n, struct GRD_HEADER *h, float *grd, double tick_gap, double tick_length, BOOLEAN tick_low, BOOLEAN tick_high, BOOLEAN tick_label, char *labels)
{
	int np, i, j, inside, ix, iy, done, n_ticks;
	double x0, y0, add, sign, dx, dy, x_back, y_back, x_front, y_front;
	double xmin, xmax, ymin, ymax, x_mean, y_mean, inc, dist, a, small, this_lon, this_lat, sa, ca, *s;
	char txt[2][2];
	
	/* The x/y coordinates in SAVE are now all projected to map inches */
	
	small = 0.1 * project_info.xmax / h->nx;
	if (tick_label) {
		txt[0][0] = labels[0];	txt[0][1] = 0;
		txt[1][0] = labels[1];	txt[1][1] = 0;
	}
	
	for (i = 0; i < n; i++) {	/* Mark polygons that have other polygons inside them */
		np = save[i].n;
		for (j = 0; save[i].do_it && j < n; j++) {
			inside = GMT_non_zero_winding (save[j].x[0], save[j].y[0], save[i].x, save[i].y, np);
			if (inside == 2) save[i].do_it = FALSE;
		}
	}
	
	/* Here, only the polygons that are innermost (containing the local max/min, will have do_it = TRUE */
	
	for (i = 0; i < n; i++) {
		if (!save[i].do_it) continue;
		np = save[i].n;
		s = (double *) GMT_memory (VNULL, (size_t)np, sizeof (double), GMT_program);
		
		/* Get rectangular surrounding box for this polygon */
		
		xmin = xmax = save[i].x[0];
		ymin = ymax = save[i].y[0];
		for (j = 1; j < np; j++) {
			xmin = MIN (xmin, save[i].x[j]);
			xmax = MAX (xmax, save[i].x[j]);
			ymin = MIN (ymin, save[i].y[j]);
			ymax = MAX (ymax, save[i].y[j]);
		}
		
		/* Find a point inside the polygon to test for high or low */
		
		GMT_xy_to_geo (&this_lon, &this_lat, 0.5 * (xmin + xmax), ymax);
		ix = irint ((this_lon - h->x_min) / h->x_inc) + 1;
		iy = irint ((h->y_max - this_lat) / h->y_inc) + 1;
		this_lon = h->x_min + ix * h->x_inc;
		done = FALSE;
		while (!done && iy < h->ny) {
			this_lat = h->y_max - iy * h->y_inc;
			GMT_geo_to_xy (this_lon, this_lat, &x0, &y0);
			inside = GMT_non_zero_winding (x0, y0, save[i].x, save[i].y, np);
			if (inside == 2)
				done = TRUE;
			else
				iy++;
		}
		if (iy == h->ny) iy--;	/* Check to make sure we don't exceed array bounds in next line. However, iy may thus be wrong */ 
		save[i].high = (grd[iy*h->nx+ix] > save[i].cval);
		
		if (save[i].high && !tick_high) continue;	/* Dont tick highs */
		if (!save[i].high && !tick_low) continue;	/* Dont tick lows */
		
		if (tick_label) {
			x_mean = save[i].x[0];
			y_mean = save[i].y[0];
		}
		for (j = 1, s[0] = 0.0; j < np; j++) {
			s[j] = s[j-1] + hypot (save[i].x[j]-save[i].x[j-1], save[i].y[j]-save[i].y[j-1]);
			if (tick_label) {
				x_mean += save[i].x[j];
				y_mean += save[i].y[j];
			}
		}
		if (s[np-1] < MIN_LENGTH) {
			GMT_free ((void *)s);
			continue;
		}
		dx = save[i].x[1] - save[i].x[0];
		dy = save[i].y[1] - save[i].y[0];
		a = atan2 (dy, dx) + M_PI_2;
		sincos (a, &sa, &ca);
		x0 = 0.5 * (save[i].x[0] + save[i].x[1]) + small * ca;
		y0 = 0.5 * (save[i].y[0] + save[i].y[1]) + small * sa;
		inside = GMT_non_zero_winding (x0, y0, save[i].x, save[i].y, np);
		sign = (inside == 2) ? 1.0 : -1.0;
		
		n_ticks = (int)(s[np-1] / tick_gap);
		if (n_ticks == 0) {
			GMT_free ((void *)s);
			continue;
		}
		inc = s[np-1] / n_ticks;
		if (tick_label) {
			x_mean /= np;
			y_mean /= np;
			GMT_text3D (x_mean, y_mean, project_info.z_level, gmtdefs.annot_font_size[0], gmtdefs.annot_font[0], txt[save[i].high], 0.0, 6, 0);
		}
		
		x_back = save[i].x[np-1];
		y_back = save[i].y[np-1];
		dist = 0.0;
		j = 0;
		add = sign * ((save[i].high) ? -1.0 : 1.0);
		GMT_setpen (&save[i].pen);
		while (j < np-1) {
			x_front = save[i].x[j+1];
			y_front = save[i].y[j+1];
			if (s[j] >= dist) {	/* Time for tick */
				dx = x_front - x_back;
				dy = y_front - y_back;
				a = atan2 (dy, dx) + add * M_PI_2;
				sincos (a, &sa, &ca);
				if (project_info.three_D) {
					GMT_xy_do_z_to_xy (save[i].x[j], save[i].y[j], project_info.z_level, &x0, &y0);
					ps_plot (x0, y0, 3);
					GMT_xy_do_z_to_xy (save[i].x[j] + tick_length * ca, save[i].y[j] + tick_length * sa, project_info.z_level, &x0, &y0);
					ps_plot (x0, y0, -2);
				}
				else {
					ps_plot (save[i].x[j], save[i].y[j], 3);
					ps_plotr (tick_length * ca, tick_length * sa, -2);
				}
				dist += inc;
			}
			x_back = x_front;
			y_back = y_front;
			j++;
		}
		GMT_free ((void *)s);
	}
}

void GMT_grd_minmax (float *a, struct GRD_HEADER *h, double xyz[2][3])
{
	int i, i_min, i_max;
	float z_min, z_max;
	double half;
	
	z_min = FLT_MAX;	z_max = -FLT_MAX;
	i_min = i_max = 0;
	half = (h->node_offset) ? 0.5 : 0.0;
	for (i = 0; i < h->nx * h->ny; i++) {
		if (GMT_is_fnan (a[i])) continue;
		if (a[i] < z_min) {
			z_min = a[i];
			i_min = i;
		}
		if (a[i] > z_max) {
			z_max = a[i];
			i_max = i;
		}
	}
	xyz[0][0] = h->x_min + (i_min%h->nx + half) * h->x_inc;
	xyz[0][1] = h->y_max - (i_min/h->nx + half) * h->y_inc;
	xyz[1][0] = h->x_min + (i_max%h->nx + half) * h->x_inc;
	xyz[1][1] = h->y_max - (i_max/h->nx + half) * h->y_inc;
	xyz[0][2] = z_min;	xyz[1][2] = z_max;
}
