/*--------------------------------------------------------------------
 *	$Id: psmask.c,v 1.4.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
 *--------------------------------------------------------------------*/
/*
 * psmask tries to achieve masking using one of two different approachs:
 * The default way of operation is: Instead of painting tiles where there's no
 * data [or where there is data] psmask uses contouring to find the polygons
 * that contain the data [or contain the regions with no data].  For many types
 * of data coverage, this results in a managable path instead of thousands of
 * tiles.  So, instead of painting polygons, psmask sets up actual clip paths.
 * As an option, the user may specify a rgb combination to fill the clipped areas.
 * To avoid having to deal with the problems that arise if the data distribution
 * is such that part of the map boundary should be part of the clip paths, we
 * internally enlarge the grid by one gridsize unit so no nodes along the edges
 * have data.  Before using the clippaths we move points outside the real region
 * onto the boundary.
 * Should the paths be too long for PostScript to handle, the user may override
 * the default with the -T switch.  Then, masking is achieved using tiling
 *
 * Author:	Paul Wessel
 * Date:	01-JUL-2000
 * Version:	3.4.1
 *
 */

#include "gmt.h"

#define IJ(i,j) ((i) + (j) * nx)

/* Global variables for contour_sub routines */

int p[5], i_off[5], j_off[5], k_off[5], offset;
unsigned int bit[32];

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

void draw_clip_contours (double *xx, double *yy, int nn, int rgb[], int id, int flag);
void dump_clip_countours (double *xx, double *yy, int nn, int id, char *file);
void shrink_clip_contours (double *x, double *y, int n, double w, double e);
int trace_clip_contours (char *grd, int nx, int *edge, double west, double north, double dx, double dy, double xinc2, double yinc2, double **xx, double **yy, int i, int j, int kk, int *max);
int clip_contours (char *grd, int nx, int ny, double west, double north, double dx, double dy, double xinc2, double yinc2, int *edge, int first, double **x, double **y, int *max);

main (int argc, char **argv)
{
	int i, j, ij, n, nm, n_edges, di, dj, ii, jj, n_expected_fields;
	int nx, ny, section, n_alloc, n_read, n_fields, one_or_zero = 1;
	
	BOOLEAN error = FALSE, first = TRUE, dump = FALSE, invert = FALSE;
	BOOLEAN end_of_clip = FALSE, find_paths = TRUE, use_poly;
	BOOLEAN map_units = FALSE, pixel = FALSE, node_only;
	
	char dfile[BUFSIZ], line[BUFSIZ];
	
	double west, east, north, south, dx, dy, radius = 0.0, *in, distance;
	double x0, y0, x1, y1, xx[4], yy[4], dummy, shrink, km_pr_deg, offset = 0.5;
	double xinc2, yinc2, idx, idy;
	
	FILE *fp = NULL;
	
	struct GMT_FILL fill;
	
	west = east = north = south = xinc2 = yinc2 = 0.0;

	argc = GMT_begin (argc, argv);

	GMT_init_fill (&fill, -1, -1, -1);
	strcpy (dfile, "mask");
	dx = dy = 0.0;
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'B':
				case 'J':
				case 'K':
				case 'O':
				case 'P':
				case 'R':
				case 'U':
				case 'V':
				case 'X':
				case 'x':
				case 'Y':
				case 'y':
				case 'c':
				case ':':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
			
				case 'b':	/* Input pairs are binary, not ascii */
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'C':
					end_of_clip = TRUE;
					break;
				case 'D':	/* Dump the polygons to files */
					strcpy (dfile, &argv[i][2]);
					dump = TRUE;
					break;
				case 'E':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'F':
					pixel = TRUE;
					break;
				case 'G':
					if (GMT_getfill (&argv[i][2], &fill)) {
						GMT_fill_syntax ('G');
						error++;
					}
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &dx, &dy);
					break;
				case 'N':
					invert = TRUE;
					break;
				case 'S':	/* Radius of influence */
					GMT_getinc(&argv[i][2], &radius, &radius);
					map_units = (argv[i][strlen(argv[i])-1] == 'k');
					break;
				case 'T':
					find_paths = FALSE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else {
			if ((fp = GMT_fopen(argv[i], GMT_io.r_mode)) == NULL) {
				fprintf (stderr, "%s: Could not open file %s\n", GMT_program, argv[i]);
				exit (EXIT_FAILURE);
			}
		}
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"psmask %s - Clipping of 2-D data sets\n\n", GMT_VERSION);
		fprintf (stderr, "usage: psmask <xyz-file> -I<xinc>[m|c][/<yinc>[m|c]] -J<params> -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-B<tickinfo>] [-C] -D<file> [-Eaz/el] [-F] [-G<fill>] [-H[<nrec>]] [-K] [-M[<flag>]]\n");
		fprintf (stderr, "\t[-N] [-O] [-P] [-S<radius>[k]] [-T] [-U[<label>]] [-V] [-X<x_shift>] [-Y<y_shift>]\n");
		fprintf (stderr, "\t[-c<ncopies>][ -:] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<xyz-file> is the datafile.  If not given, read standard input\n");
		fprintf (stderr, "\t-I sets grid increments; enter xinc, optionally xinc/yinc.\n");
		fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or s] to xinc or yinc to indicate minutes [or seconds]\n");
		GMT_explain_option ('j');
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		GMT_explain_option ('b');
		fprintf (stderr, "\t-C means stop existing clip-path.  No other options required\n");
		fprintf (stderr, "\t-D dumps the clip-paths to files using the prefix <file>_ [mask_]\n");
		fprintf (stderr, "\t   Ignored if -T is specified\n");
		fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
		fprintf (stderr, "\t-F Force pixel registration [Default is gridline registration]\n");
		fprintf (stderr, "\t-G sets fill color r/g/b [Default is no fill]\n");
		GMT_explain_option ('H');
		GMT_explain_option ('K');
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N will invert the sense of the clipping [or tiling]\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr, "\t-S Search radius: Nodes inside circles of <radius> centered on\n");
		fprintf (stderr, "\t   the input data poins are considered to be reliable estimates of the surface\n");
		fprintf (stderr, "\t   Default is -S0, i.e., only the nearest node is considered reliable\n");
		fprintf (stderr, "\t   Append k to indicate km (hence -R -I are in degrees)\n");
		fprintf (stderr, "\t-T will paint tiles.  [Default will trace data outline]\n");
		fprintf (stderr, "\t   If set you must also specify a color/fill with -G\n");
		fprintf (stderr, "\t   !May only be used with linear, Mercator, and basic cylindrical projections!\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		GMT_explain_option ('X');
		GMT_explain_option ('c');
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		GMT_explain_option ('o');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!end_of_clip) {
		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 (!project_info.region_supplied) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
			error++;
		}
		if (!find_paths && !RECT_GRATICULE) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Only available with Linear, Mercator, or basic cylindrical projections\n", GMT_program);
			error++;
		}
		if (!find_paths && fill.rgb[0] == -1) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -T option:  Must also specify a tile color with -G\n", GMT_program);
			error++;
		}
		if (dx <= 0.0 || dy <= 0.0) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify positive increments\n", GMT_program);
			error++;
		}
		if (radius < 0.0) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Radius must be positive\n", GMT_program);
			error++;
		}
		if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 2;
		if (GMT_io.binary[0] && GMT_io.ncol[0] < 2) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 2 columns\n", GMT_program);
			error++;
		}
		if (GMT_io.binary[0] && gmtdefs.io_header) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
			error++;
		}
	}

	if (error) exit (EXIT_FAILURE);

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

	if (GMT_io.binary[0] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]);
	}

	if (pixel) {
		xinc2 = 0.5 * dx;
		yinc2 = 0.5 * dy;
		one_or_zero = 0;
		offset = 0.0;
	}
	idx = 1.0 / dx;
	idy = 1.0 / dy;

	if (!end_of_clip) GMT_map_setup (west, east, south, north);

	if (!project_info.x_off_supplied && gmtdefs.overlay) gmtdefs.x_origin = 0.0;
	if (!project_info.y_off_supplied && gmtdefs.overlay) gmtdefs.y_origin = 0.0;

	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 (!end_of_clip) {	/* Start new clip_path */
	
		if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);

		if (gmtdefs.verbose) fprintf (stderr, "%s: Allocate memory, read and process data file\n", GMT_program);
	
		/* Enlarge region by 1 row/column */

		west -= dx;	east += dx;	south -= dy;	north += dy;

		nx = irint ((east - west)   * idx) + one_or_zero;
		ny = irint ((north - south) * idy) + one_or_zero;
	
		nm = nx * ny;
		grd = (char *) GMT_memory (VNULL, (size_t)nm, sizeof (char), GMT_program);
	
		if (map_units) {
			km_pr_deg = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0;
			shrink = cos (0.5 * (south + north) * D2R);
			di = (int)ceil (radius / (km_pr_deg * dx * shrink));
			dj = (int)ceil (radius / (km_pr_deg * dy));
		}
		else {
			di = irint (0.5 * radius / dx);
			dj = irint (0.5 * radius / dy);
		}
		node_only = (di == 0 && dj == 0);

		if (fp == NULL) {
			fp = GMT_stdin;
#ifdef SET_IO_MODE
			GMT_setmode (0);
#endif
		}
		
		if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp);

		n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 2;
		n_read = 0;
		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {

			while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) { /* Skip headers */
				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
			if ((GMT_io.status & GMT_IO_EOF)) continue;	/* At EOF */

			if (GMT_io.status & GMT_IO_MISMATCH) {
				fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n_read);
				exit (EXIT_FAILURE);
			}
			n_read++;

			/* Determine the node closest to the data point */

			i = (int)floor (((in[0] - west) * idx) + offset);
			if (i < 0 || i >= nx) continue;
			j = (int)floor (((north - in[1]) * idy) + offset);
			if (j < 0 || j >= ny) continue;

			if (node_only) {
				grd[j*nx+i] = 1;
			}
			else {

				/* Set coordinate of this node */

				x0 = west + i * dx + xinc2;
				y0 = north - j * dy - yinc2;
		
				/* Set this and all nodes within radius distance to 1 */

				for (ii = i - di; ii <= i + di; ii++) {
					if (ii < 0 || ii >= nx) continue;
					x1 = west + ii * dx + xinc2;
					for (jj = j - dj; jj <= j + dj; jj++) {
						if (jj < 0 || jj >= ny) continue;
						y1 = north - jj * dy - yinc2;
						distance = (map_units) ? km_pr_deg * hypot ((x1 - x0) * shrink, y1 - y0) : hypot (x1 - x0, y1 - y0);
						if (distance > radius) continue;
						grd[jj*nx+ii] = 1;
					}
				}
			}
		}

		if (fp != GMT_stdin) GMT_fclose (fp);
		if (gmtdefs.verbose) fprintf (stderr, "%s: Read %d data points\n", GMT_program, n_read);

		if (invert) for (i = 0; i < nm; i++) grd[i] = 1 - grd[i];	/* Reverse sense of test */

		/* Force perimeter nodes to be FALSE */
		
		for (i = 0, ij = (ny-1) * nx; i < nx; i++) grd[i] = grd[i+ij] = FALSE;
		for (j = 0; j < ny; j++) grd[j*nx] = grd[(j+1)*nx-1] = FALSE;
		
		if (find_paths) {	/* Must trace the outline of ON/OFF values in grd */
			x = (double *) GMT_memory (VNULL, (size_t)GMT_CHUNK, sizeof (double), GMT_program);
			y = (double *) GMT_memory (VNULL, (size_t)GMT_CHUNK, sizeof (double), GMT_program);
			n_alloc = GMT_CHUNK;
	
			n_edges = ny * (int )ceil (nx / 16.0);
			edge = (int *) GMT_memory (VNULL, (size_t)n_edges, sizeof (int), GMT_program);
	
			if (frame_info.plot) {
				ps_setpaint (gmtdefs.basemap_frame_rgb);
				GMT_map_basemap ();
				ps_setpaint (gmtdefs.background_rgb);
			}

			if (gmtdefs.verbose) fprintf (stderr, "%s: Tracing the clip path\n", GMT_program);
	
			section = 0;
			first = TRUE;
			while ((n = clip_contours (grd, nx, ny, west, north, dx, dy, xinc2, yinc2, edge, first, &x, &y, &n_alloc)) > 0) {
				shrink_clip_contours (x, y, n, west, east);
				if (dump) dump_clip_countours (x, y, n, section, dfile);
				draw_clip_contours (x, y, n, fill.rgb, section, first);
				first = FALSE;
				section++;
			}
	
			draw_clip_contours (x, y, 0, fill.rgb, section, 2);	/* Activate clip-path */
		
			GMT_free ((void *)edge);
			GMT_free ((void *)x);
			GMT_free ((void *)y);
		}
		else {	/* Just paint tiles */
			if (gmtdefs.verbose) fprintf (stderr, "%s: Tiling...\n", GMT_program);
			x_grid = (double *) GMT_memory (VNULL, (size_t)nx, sizeof (double), GMT_program);
			y_grid = (double *) GMT_memory (VNULL, (size_t)ny, sizeof (double), GMT_program);
			for (i = 2; i < nx - 1; i++) GMT_geo_to_xy (west + (i - 0.5) * dx + xinc2, south + dy + yinc2, &x_grid[i], &dummy);
			for (j = 2; j < ny - 1; j++) GMT_geo_to_xy (west + dx + xinc2, north - (j - 0.5) * dy - yinc2, &dummy, &y_grid[j]);
			GMT_geo_to_xy (west + dx + xinc2, north - dy - yinc2, &x_grid[1], &y_grid[1]);
			GMT_geo_to_xy (east - dx - xinc2, south + dy + yinc2, &x_grid[nx-1], &y_grid[ny-1]);
			
			use_poly = (fill.use_pattern || project_info.three_D);

			for (j = 1; j < ny; j++) {
				ij = j * nx + 1;
				for (i = 1; i < nx; i++, ij++) {
					if (((int)grd[ij]) == 0) continue;
					
					if (use_poly) {
						xx[0] = xx[3] = x_grid[i];	xx[1] = xx[2] = x_grid[i+1];
						yy[0] = yy[1] = y_grid[j+1];	yy[2] = yy[3] = y_grid[j];
						if (project_info.three_D) GMT_2D_to_3D (xx, yy, 4);
						GMT_fill (xx, yy, 4, &fill, FALSE);
					}
					else
						ps_rect (x_grid[i], y_grid[j+1], x_grid[i+1], y_grid[j], fill.rgb, FALSE);
				}
			}
			
			GMT_free ((void *)x_grid);
			GMT_free ((void *)y_grid);
			
			if (frame_info.plot) {
				ps_setpaint (gmtdefs.basemap_frame_rgb);
				GMT_map_basemap ();
			}
		}
		
		if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);

		GMT_free ((void *)grd);
		if (find_paths && gmtdefs.verbose) fprintf (stderr, "%s: clipping on!\n", GMT_program);
	}
	else {	/* Just undo previous clip-path */
		ps_clipoff ();
		if (frame_info.plot) {
			ps_setpaint (gmtdefs.basemap_frame_rgb);
			GMT_map_basemap ();
		}
		if (gmtdefs.verbose) fprintf (stderr, "%s: clipping off!\n", GMT_program);
	}
	
	ps_setpaint (gmtdefs.background_rgb);
	
	ps_plotend (gmtdefs.last_page);
	
	
	GMT_end (argc, argv);
}

void draw_clip_contours (double *xx, double *yy, int nn, int rgb[], int id, int flag)
{
	int i;
	double x, y;
	char comment[64];
	
	if (nn < 2 && flag < 2) return;
	
	for (i = 0; i < nn; i++) {
		x = xx[i];	y = yy[i];
		GMT_geo_to_xy (x, y, &xx[i], &yy[i]);
	}
	nn = GMT_compact_line (xx, yy, nn, FALSE, 0);
	
	if (project_info.three_D) GMT_2D_to_3D (xx, yy, nn);
	
	if (nn > 0) {
		sprintf (comment, "Start of clip path sub-segment %d", id);
		ps_comment (comment);
	}
	ps_clipon (xx, yy, nn, rgb, flag);
	if (nn > 0) {
		sprintf (comment, "End of clip path sub-segment %d", id);
		ps_comment (comment);
	}		
}

void dump_clip_countours (double *xx, double *yy, int nn, int id, char *file)
{
	int i;
	double out[2];
	char fname[BUFSIZ];
	FILE *fp;
	
	if (nn < 2) return;
	
	if (GMT_io.binary[1])
		sprintf (fname, "%s_%d.b", file, id);
	else
		sprintf (fname, "%s_%d.xy", file, id);

	if ((fp = GMT_fopen (fname, GMT_io.w_mode)) == NULL) {
		fprintf (stderr, "%s: Unable to create file %s - exiting\n", GMT_program, fname);
		exit (EXIT_FAILURE);
	}
	for (i = 0; i < nn; i++) {
		out[0] = xx[i];	out[1] = yy[i];
		GMT_output (fp, 2, out);
	}
	GMT_fclose (fp);
}

int clip_contours (char *grd, int nx, int ny, double west, double north, double dx, double dy, double xinc2, double yinc2, int *edge, int first, double **x, double **y, int *max)
{
	/* The routine finds the zero-contour in the grd dataset.  it assumes that
	 * no node has a value exactly == 0.0.  If more than max points are found
	 * trace_clip_contours will try to allocate more memory in blocks of GMT_CHUNK points
	 */
	 
	 static int i0, j0;
	 int i, j, ij, n = 0, n_edges, edge_word, edge_bit;
	 BOOLEAN go_on = TRUE;
	 
	 
	 n_edges = ny * (int) ceil (nx / 16.0);
	 offset = n_edges / 2;
	 
	 /* Reset edge-flags to zero, if necessary */
	 if (first) {
	 	i0 = 0;	/* Begin with upper left bin which is i = 0 and j = 1 */
	 	j0 = 1;
		p[0] = p[4] = 0;	p[1] = 1;	p[2] = 1 - nx;	p[3] = -nx;
		i_off[0] = i_off[2] = i_off[3] = i_off[4] = 0;	i_off[1] =  1;
		j_off[0] = j_off[1] = j_off[3] = j_off[4] = 0;	j_off[2] = -1;
		k_off[0] = k_off[2] = k_off[4] = 0;	k_off[1] = k_off[3] = 1;
		for (i = 1, bit[0] = 1; i < 32; i++) bit[i] = bit[i-1] << 1;
	 }

	/* Loop over interior boxes */
		
	for (j = j0; go_on && j < ny; j++) {
		ij = i0 + j * nx;
		for (i = i0; go_on && i < nx-1; i++, ij++) {	/* nx-1 since the last bin starts at nx-2 and ends at nx-1 */
			edge_word = ij / 32 + offset;
			edge_bit = ij % 32;
			if (!(edge[edge_word] & bit[edge_bit]) && ((grd[ij]+grd[ij-nx]) == 1)) { /* Start tracing countour */
				*x[0] = west + i*dx + xinc2;
				*y[0] = north - (j-0.5) * dy - yinc2;
				edge[edge_word] |= bit[edge_bit];
				n = trace_clip_contours (grd, nx, edge, west, north, dx, dy, xinc2, yinc2, x, y, i, j, 3, max);
				go_on = FALSE;
				i0 = i + 1;
				j0 = j;	/* Return to finish this row later */
			}
		}
		if (go_on) i0 = 0;	/* Go to start of next row unless we found something */
	}
	
	return (n);
}

int trace_clip_contours (char *grd, int nx, int *edge, double west, double north, double dx, double dy, double xinc2, double yinc2, double **xx, double **yy, int i, int j, int kk, int *max)
{
	int n = 1, k, k0, ij, ij0, n_cuts, kk_opposite, first_k, more;
	int edge_word, edge_bit, m;
	double xk[4], yk[4], x0, y0;
	
	m = *max - 2;
	
	more = TRUE;
	do {
		ij = i + j * nx;
		x0 = west + i * dx + xinc2;
		y0 = north - j * dy - yinc2;
		n_cuts = 0;
		k0 = kk;

		for (k = 0; k < 4; k++) {	/* Loop over box sides */

			/* Skip where we already have a cut (k == k0) */
			
			if (k == k0) continue;
			
			/* Skip edge already has been used */
			
			ij0 = IJ (i + i_off[k], j + j_off[k]);
			edge_word = ij0 / 32 + k_off[k] * offset;
			edge_bit = ij0 % 32;
			if (edge[edge_word] & bit[edge_bit]) continue;
			
			/* Skip if no zero-crossing on this edge */
			
			if ((grd[ij+p[k+1]] + grd[ij+p[k]]) != 1) continue;
			
			/* Here we have a crossing */
			
			if (k%2) {	/* Cutting a S-N line */
				if (k == 1) {
					xk[1] = x0 + dx;
					yk[1] = y0 + 0.5*dy;
				}
				else {
					xk[3] = x0;
					yk[3] = y0 + 0.5*dy;
				}
			}
			else {	/* Cutting a E-W line */
				if (k == 0) {
					xk[0] = x0 + 0.5*dx;
					yk[0] = y0;
				}
				else {
					xk[2] = x0 + 0.5*dx;
					yk[2] = y0 + dy;
				}
			}
			kk = k;
			n_cuts++;
		}
		
		if (n > m) {	/* Must try to allocate more memory */
			*max += GMT_CHUNK;
			m += GMT_CHUNK;
			*xx = (double *) GMT_memory ((void *)*xx, (size_t)(*max), sizeof (double), "trace_clip_contours");
			*yy = (double *) GMT_memory ((void *)*yy, (size_t)(*max), sizeof (double), "trace_clip_contours");
		}
		if (n_cuts == 0) {	/* Close interior contour and return */
			if (fmod ((*xx[0] - xinc2), dx) == 0.0)	/* On side 1 or 3 */
				first_k = ((*xx[0] - x0) == 0.0) ? 3 : 1;
			else 	/* On side 0 or 2 */
				first_k = ((*yy[0] - y0) == 0.0) ? 0 : 2;
			kk_opposite = (first_k + 2) % 4;
			if (k0 != kk_opposite) {
				(*xx)[n] = x0 + 0.5*dx;
				(*yy)[n] = y0 + 0.5*dy;
				n++;
			}
			(*xx)[n] = (*xx)[0];
			(*yy)[n] = (*yy)[0];
			n++;
			more = FALSE;
		}
		else if (n_cuts == 1) {	/* Draw a line to this point and keep tracing */
			/* Add center of box if this and previous cut NOT on opposite edges */
			kk_opposite = (k0 + 2) % 4;
			if (kk != kk_opposite) {
				(*xx)[n] = x0 + 0.5*dx;
				(*yy)[n] = y0 + 0.5*dy;
				n++;
			}
			(*xx)[n] = xk[kk];
			(*yy)[n] = yk[kk];
			n++;
		}
		else {	/* Saddle point, we decide to connect to the point nearest previous point */
			kk = (k0 + 1)%4;	/* Pick next edge since it is arbitrarely where we go */
			/* First add center of box */
			(*xx)[n] = x0 + 0.5*dx;
			(*yy)[n] = y0 + 0.5*dy;
			n++;
			(*xx)[n] = xk[kk];
			(*yy)[n] = yk[kk];
			n++;
		}
		if (more) {	/* Mark this edge as used */
			ij0 = IJ (i + i_off[kk], j + j_off[kk]);
			edge_word = ij0 / 32 + k_off[kk] * offset;
			edge_bit = ij0 % 32;
			edge[edge_word] |= bit[edge_bit];	
		}
		
		/* Get next box (i,j,kk) */
		
		i -= (kk-2)%2;
		j -= (kk-1)%2;
		kk = (kk+2)%4;
		
	} while (more);
	return (n);
}

void shrink_clip_contours (double *x, double *y, int n, double w, double e)
{
	/* Moves outside points to boundary */
	int i;

	for (i = 0; i < n; i++) {
		if (x[i] < w)
			x[i] = w;	
		if (x[i] > e)
			x[i] = e;
		if (y[i] < project_info.s)
			y[i] = project_info.s;
		if (y[i] > project_info.n)
			y[i] = project_info.n;
	}
}