/*--------------------------------------------------------------------
 *	$Id: grdcontour.c,v 1.3.4.4 2002/02/27 17:58:55 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
 *--------------------------------------------------------------------*/
/*
 * grdcontour reads a 2-D grd file and contours it. This algoritm 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:	3.4.1
 *
 */

#include "gmt.h"

#define MIN_LENGTH 0.01

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

struct LABEL {
	double x, y;
	double angle;
	char label[32];
	struct LABEL *next_label, *prev_label;
} *anchor, *old_label;

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

void dump_contour (double *xx, double *yy, int nn, double cval, int id, BOOLEAN interior, char *file);
void draw_contour (double *xx, double *yy, int nn, double cval, char ctype, double cangle, int closed, double gap, int half_width, char *unit, int box, int label_font_size);
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 plot_labels (int size, int box, int rgb[]);

int int_cont_count = 0, ext_cont_count = 0;

main (int argc, char **argv)
{
	int i, j, k, n, nm, c, n_edges, label_font_size = 9, b_rgb[3], n_alloc;
	int section, n_contours, smooth_factor = 0, side, id, box = 1, bad;
	int off, nx, ny, half = 5, n_save = 0, got, ncut = 0, rgb[3];
	
	BOOLEAN error = FALSE, first = TRUE, dump = FALSE, anot = FALSE, interior, begin;
	BOOLEAN fix_angle = FALSE, do_ticks = FALSE, add_unit = FALSE, use_cpt_colors = FALSE, use_cpt_anot = 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, c_unit[32], txt_a[32], txt_b[32];
	
	double c_int = 0.0, c_low, c_high, take_out, aval, a_int = 0.0, tick_gap = 0.2;
	double anot_dist = 4.0, west = 0.0, east = 0.0, south = 0.0, north = 0.0, tick_length = 0.04;
	double *contour, cval, mult = 1.0, shift = 0.0, min, max, small, label_angle = 0.0;
	double small_x, small_y, data_west, data_east, data_south, data_north, tmp, *cont_angle;
	
	FILE *fpc = NULL;
	
	struct GRD_HEADER header;
	
	struct GMT_PEN pen[2];

	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] = c_unit[0] = 0;
	
	if (gmtdefs.measure_unit == GMT_CM) {
		anot_dist = 10.0 / 2.54;
		tick_gap = 0.5 / 2.54;
		tick_length = 0.1 / 2.54;
	}

	for (i = 0; i < 3; i++) b_rgb[i] = gmtdefs.page_rgb[i];	/* Default box color is page color */
	
	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 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'A':
					if (argv[i][2] == '-') use_cpt_anot = FALSE;
					n = sscanf (&argv[i][2], "%lf", &a_int);
					
					for (j = 2, bad = 0; argv[i][j] && argv[i][j] != 'f'; j++);
					if (argv[i][j])	{ /* Found font size option */
						label_font_size = atoi (&argv[i][j+1]);
						if (label_font_size <= 0) bad++;
					}
						
					for (j = 2; argv[i][j] && argv[i][j] != 'a'; j++);
					if (argv[i][j])	{ /* Found fixed angle option */
						label_angle = atof (&argv[i][j+1]);
						fix_angle = TRUE;
						if (label_angle <= -90.0 || label_angle > 180.0) bad++;
					}
						
					for (j = 2; argv[i][j] && argv[i][j] != '/'; j++);
					if (argv[i][j] && GMT_getrgb (&argv[i][j+1], b_rgb)) bad++;
					anot = TRUE;
					if (strchr (argv[i], 'o')) box = 2;
					if (strchr (argv[i], 't')) box = 0;
					if (bad) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -A option.  Correct syntax:\n", GMT_program);
						fprintf (stderr, "\t-A[-][aint][f<fontsize>][a<fixedangle>][/<r/g/b>][o]\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':
					sscanf (&argv[i][2], "%[^/]/%d", txt_a, &half);
					anot_dist = GMT_convert_units (txt_a, GMT_INCH);
					half /= 2;
					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':
					add_unit = TRUE;
					if (argv[i][2]) strcpy (c_unit, &argv[i][2]);
					break;
				case 'S':
					smooth_factor = atoi (&argv[i][2]);
					break;
				case 'Q':
					ncut = atoi (&argv[i][2]);
					break;
				case 'T':
					do_ticks = TRUE;
					if (argv[i][2]) {
						if (argv[i][2] == '+')
							tick_high = TRUE, j = 1;
						else if (argv[i][2] == '-')
							tick_low = TRUE, j = 1;
						else
							tick_high = tick_low = TRUE, j = 0;
						if (argv[i][2+j] != ':') {
							n = sscanf (&argv[i][2+j], "%[^/]/%s", 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]) {
							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;
						}
					}
					else
						tick_high = tick_low = TRUE;	/* Do default */
					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> [-A[-][<anot_int>][f<fontsize>][/r/g/b][a<angle>][o]]\n");
		fprintf (stderr, "\t[-B<tickinfo>] [-D<dumpfile>] [-Eaz/el] [-G<gap>[c|i|m|p]/<npoints>] [-K] [-L<Low/high>]\n");
		fprintf (stderr, "\t[-M[<flag>]] [-N[unit]] [-O] [-P] [-Q<cut>] [-R<west/east/south/north>] [-S<smooth>] [-T[+|-][<gap>[c|i|m|p]/<length>[c|i|m|p]][:LH]] [-U[<label>]] [-V]\n");
		fprintf (stderr, "\t[-W[+]<type><pen>] [-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   -Fixed contour interval\n");
		fprintf (stderr, "\t   -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   -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 information. [Default is no annotation].\n");
		fprintf (stderr, "\t   Append - to disable all contour anotations\n");
		fprintf (stderr, "\t   Append f followed by desired font size in points [Default is 9].\n");
		fprintf (stderr, "\t   Append /r/g/b to change color of text box [Default is %d/%d/%d]\n", b_rgb[0], b_rgb[1], b_rgb[2]);
		fprintf (stderr, "\t   Append o to draw outline of text box [Default is no outline]\n");
		fprintf (stderr, "\t   Append t for transparent text box [Default is opaque]\n");
		fprintf (stderr, "\t   Append a<angle> to force anotations at this fixed angle [Default follows contour]\n");
		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 countour 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 Gap between anotations in %s and number of (x,y) values used to find L2 label angle [Default = %lg/10]\n", GMT_unit_names[gmtdefs.measure_unit], anot_dist);
		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");
		fprintf (stderr, "\t-N add z-unit from grdfile to contour labels; append unit to override name.\n");
		fprintf (stderr, "\t   Start unit with - for no space between anotation and unit\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 (in %s) to change defaults [%lg/%lg]\n", GMT_unit_names[gmtdefs.measure_unit], tick_gap, tick_length);
		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 anotated contours and\n");
		fprintf (stderr, "\t   'c' for regular contours.  The default settings are\n");
		fprintf (stderr, "\t   Contour pen:  width = %lgp, 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 = %lgp, 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 subract <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 ('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 (anot_dist <= 0.0 || half <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option.  Correct syntax:\n", GMT_program);
		fprintf (stderr, "\t-G<anot_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 .gmtcommands */


	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 (add_unit && !c_unit[0]) strcpy (c_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 */

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

	/* 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 = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program);
	if (GMT_read_grd (grdfile, &header, grd, 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);
	}

	anchor = old_label = (struct LABEL *) GMT_memory (VNULL, (size_t)1, sizeof (struct LABEL), 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);
	
	if (!(mult == 1.0 && shift == 0.0)) {
		if (gmtdefs.verbose) fprintf (stderr, "%s: Subtracting %lg and multiplying grid by %lg\n", GMT_program, shift, mult);
		for (i = 0; i < nm; i++) grd[i] = (float)((grd[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, 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);
	
	if (anot) {	/* Wants anotated 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_anot)
				cont_type[c] = 'C';
			else if (GMT_lut[i].anot)
				cont_type[c] = 'A';
			else
				cont_type[c] = (anot) ? 'A' : 'C';
			cont_angle[c] = (fix_angle) ? 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_anot)
			cont_type[c] = 'C';
		else if (GMT_lut[GMT_n_colors-1].anot & 2)
			cont_type[c] = 'A';
		else
			cont_type[c] = (anot) ? 'A' : 'C';
		cont_angle[c] = (fix_angle) ? 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 (anot && (contour[n_contours] - aval) > SMALL) aval += a_int;
			cont_type[n_contours] = (fabs (contour[n_contours] - aval) < SMALL) ? 'A' : 'C';
			cont_angle[n_contours] = (fix_angle) ? 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 redundent 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);
	
	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 %lg contour\n", GMT_program, cval);
		take_out = (first) ? cval : contour[c] - contour[c-1];
		for (i = 0; i < nm; i++) {
			grd[i] -= (float)take_out;
			if (grd[i] == 0.0) grd[i] += (float)small;
		}
		
		section = 0;
		first = FALSE;
		id = (cont_type[c] == 'A' || cont_type[c] == 'a') ? 1 : 0;

		GMT_setpen (&pen[id]);
		if (use_cpt_colors) {	/* Set pen color according to cpt file */
			GMT_get_rgb24 (cval, rgb);
			ps_setpaint (rgb);
		}
		
		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];
				}
				else
					interior = FALSE;
				if (dump && n >= ncut) dump_contour (x, y, n, cval, section, interior, dfile);
				if (n >= ncut) draw_contour (x, y, n, cval, cont_type[c], cont_angle[c], interior, anot_dist, half, c_unit, box, label_font_size);
				if (cont_do_tick[c] && interior && n >= ncut) {
					save[n_save].x = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
					save[n_save].y = (double *) GMT_memory (VNULL, (size_t)n, sizeof (double), GMT_program);
					memcpy ((void *)save[n_save].x, (void *)x, (size_t)(n * sizeof (double)));
					memcpy ((void *)save[n_save].y, (void *)y, (size_t)(n * sizeof (double)));
					save[n_save].n = n;
					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);
					}
				}
				if (n >= ncut) 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);
		GMT_read_grd_info (grdfile, &header);
		GMT_read_grd (grdfile, &header, grd, data_west, data_east, data_south, data_north, GMT_pad, FALSE);
		if (!(mult == 1.0 && shift == 0.0)) {
			for (i = 0; i < nm; i++) grd[i] = (float)((grd[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);
		}
		sort_and_plot_ticks (save, n_save, &header, grd, 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);
	}
	
	plot_labels (label_font_size, box, b_rgb);
	
	if (frame_info.plot) {
		ps_setpaint (gmtdefs.basemap_frame_rgb);
		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 *)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 dump_contour (double *xx, double *yy, int nn, double cval, int id, BOOLEAN interior, char *file)
{
	int i;
	char fname[BUFSIZ], format[80], suffix[4];
	double out[3];
	FILE *fp;
	
	if (nn < 2) return;
	
	out[2] = cval;
	(GMT_io.binary[1]) ? strcpy (suffix, "b") : strcpy (suffix, "xyz");
	sprintf (format, "%s\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
	if (!GMT_io.binary[1] && GMT_io.multi_segments) {
		if (GMT_io.multi_segments == 2) {	/* Must create file the first time around */
			fp = GMT_fopen (file, "w");
			GMT_io.multi_segments = TRUE;
		}
		else	/* Later we append to it */
			fp = GMT_fopen (file, "a+");
		sprintf (GMT_io.segment_header, "%c %lg contour\n", GMT_io.EOF_flag, cval);
		GMT_write_segmentheader (fp, 3);
	}
	else {
		if (interior) {
			if (file[0] == '-')	/* Running numbers only */
				sprintf (fname, "C%d_i.%s", int_cont_count++, suffix);
			else
				sprintf (fname, "%s_%lg_%d_i.%s", file, cval, id, suffix);
		}
		else {
			if (file[0] == '-')	/* Running numbers only */
				sprintf (fname, "C%d_e.%s", ext_cont_count++, suffix);
			else
				sprintf (fname, "%s_%lg_%d.%s", file, cval, id, suffix);
		}
		fp = GMT_fopen (fname, GMT_io.w_mode);
	}
	for (i = 0; i < nn; i++) {
		out[0] = xx[i];	out[1] = yy[i];
		GMT_output (fp, 3, out);
	}
	GMT_fclose (fp);
}

void draw_contour (double *xx, double *yy, int nn, double cval, char ctype, double cangle, int closed, double gap, int half_width, char *unit, int box, int label_font_size)
{
	int i, j, start = 0;
	double dist, angle, dx, dy, width, sum_x2, sum_xy, one = 1.0;
	double label_height, label_width, x, y;
	char label[128], format[128];
	struct LABEL *new_label;
	
	if (nn < 2) return;
		
	sprintf (format, "%lg contour", cval);
	ps_comment (format);
	
	GMT_n_plot = GMT_geo_to_xy_line (xx, yy, nn);
	if (box > 0) GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot);	/* Draw entire contour for opaque label boxes */
	
	if (ctype == 'A' || ctype == 'a') {	/* Annotated contours */
		GMT_get_format (cval, unit, format);
		sprintf (label, format, cval);
		label_height = label_width = label_font_size * GMT_u2u[GMT_PT][GMT_INCH] * 0.58;  /* Constants from GMT_textbox3d()KJ*/
		label_width *= 0.81 * strlen(label) * 0.78;
		label_width *= 0.66;  /* was still too wide KJ*/
		dist = (closed) ? gap - one : 0.0;	/* Label closed contours longer than 1 inch */
		for (i = 1; i < GMT_n_plot; i++) {
			if (GMT_pen[i] == 3) continue;

			dx = GMT_x_plot[i] - GMT_x_plot[i-1];
			if (fabs (dx) > (width = GMT_half_map_width (GMT_y_plot[i-1]))) {
				width *= 2.0;
				dx = copysign (width - fabs (dx), -dx);
				if (GMT_x_plot[i] < width)
					GMT_x_plot[i-1] -= width;
				else
					GMT_x_plot[i-1] += width;
			}
			dy = GMT_y_plot[i] - GMT_y_plot[i-1];
			dist += hypot (dx, dy);
			if (dist > gap) {	/* Time for label */
				new_label = (struct LABEL *) GMT_memory (VNULL, (size_t)1, sizeof (struct LABEL), GMT_program);
				new_label->x = 0.5 * (GMT_x_plot[i-1] + GMT_x_plot[i]);
				new_label->y = 0.5 * (GMT_y_plot[i-1] + GMT_y_plot[i]);
				strcpy (new_label->label, label);
				sum_x2 = sum_xy = 0.0;
				for (j = i - half_width; j <= i + half_width; j++) {	/* L2 fit for slope */
					if (j < 0 || j >= GMT_n_plot) continue;
					dx = GMT_x_plot[j] - new_label->x;
					dy = GMT_y_plot[j] - new_label->y;
					sum_x2 += dx * dx;
					sum_xy += dx * dy;
				}
				if (GMT_is_dnan (cangle)) {	/* Must calculate label angle */
					angle = d_atan2 (sum_xy, sum_x2) * R2D;
					if (angle < 0.0) angle += 360.0;
					if (angle > 90.0 && angle < 270) angle -= 180.0;
				}
				else
					angle = cangle;
				new_label->angle = angle;
				new_label->prev_label = old_label;
				dist = 0.0;
				old_label->next_label = new_label;
				old_label = new_label;

				if (box == 0) {	/* Transpartent label box */
					/* new segments of the contour line from label to label KJ*/
					j = 1;
					x = GMT_x_plot[i];
					y = GMT_y_plot[i];
	   				while ((GMT_x_plot[i-j]> (x - label_width)) && (GMT_x_plot[i-j] < (x + label_width)) &&
	                        		(GMT_y_plot[i-j] > (y - label_height)) && (GMT_y_plot[i-j] < (y+label_height))) j++;
					GMT_plot_line (&GMT_x_plot[start], &GMT_y_plot[start], &GMT_pen[start], i-j-start); 
					start = i + 1;
	   				while ( (GMT_x_plot[start] > (x - label_width)) && (GMT_x_plot[start] < (x + label_width)) &&
	                      	 		(GMT_y_plot[start] > (y - label_height)) && (GMT_y_plot[start] < (y + label_height))) start++;
					GMT_pen[start] = 3;	/* Since we may have moved up a bit */
				}
			}
			GMT_x_plot[i-1] = GMT_x_plot[i];	GMT_y_plot[i-1] = GMT_y_plot[i];
		}
		if (box == 0) GMT_plot_line (&GMT_x_plot[start], &GMT_y_plot[start], &GMT_pen[start], GMT_n_plot-start); /*KJ the last segment, supposed not to end with a label*/
	}
	else if (box == 0) {   /* just one line , no holes for labels */
		GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot);
	}
}

void plot_labels (int size, int box, int rgb[])
{
	/* box = 1: white box only, box = 2: white box + draw outline */
        /* box = 0: no box (append 't' instead of 'o' in -A parameter) */
	double dx, dy;
	struct LABEL *this, *old;
	
	dx = 0.5 * gmtdefs.anot_offset;
	dy = 0.05 * gmtdefs.anot_offset;
	box--;
	ps_comment ("Contour annotations:");
	
	ps_setpaint (gmtdefs.basemap_frame_rgb);
	
	for (old = anchor; box >= 0 && old->next_label; old = old->next_label) {	/* First draw boxes if not 3-D*/
		this = old->next_label;
		GMT_textbox3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, box, dx, dy, rgb);
	}
	for (old = anchor; old->next_label; old = old->next_label) { /* Then labels */
		this = old->next_label;
		GMT_text3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, 0);
	}
		
	ps_setpaint (gmtdefs.background_rgb);
	
	this = anchor;
	while (this) {	/* Free memory */
		old = this;
		this = old->next_label;
		GMT_free ((void *)old);
	}
}

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, *s;
	char txt[2][2];
	
	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 */
		
		ix = irint ((0.5 * (xmin + xmax) - h->x_min) / h->x_inc) + 1;
		iy = irint ((h->y_max - ymax) / h->y_inc) + 1;
		x0 = h->x_min + ix * h->x_inc;
		done = FALSE;
		while (!done && iy < h->ny) {
			y0 = h->y_max - iy * h->y_inc;
			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 dont 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 */
		
		GMT_geo_to_xy (save[i].x[0], save[i].y[0], &save[i].x[0], &save[i].y[0]);
		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++) {
			GMT_geo_to_xy (save[i].x[j], save[i].y[j], &save[i].x[j], &save[i].y[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;
		x0 = 0.5 * (save[i].x[0] + save[i].x[1]) + small * cos (a);
		y0 = 0.5 * (save[i].y[0] + save[i].y[1]) + small * sin (a);
		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.anot_font_size, gmtdefs.anot_font, 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;
				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 * cos (a), save[i].y[j] + tick_length * sin (a), 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 * cos (a), tick_length * sin (a), 2);
				}
				dist += inc;
			}
			x_back = x_front;
			y_back = y_front;
			j++;
		}
		GMT_free ((void *)s);
	}
}