/*--------------------------------------------------------------------
 *	$Id: gmtselect.c,v 1.5.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
 *--------------------------------------------------------------------*/
/*
 * gmtselect is a general-purpose spatial filter.  Data pass or fail based
 * on one or more conditions.  5 conditions may be set:
 *
 *	1. Only data inside a rectangular area may pass
 *	2. Only data within a certain distance from given points may pass
 *	3. Only data within a certain distance from given lines may pass
 *	4. Only data within given polygons may pass
 *	5. Only data within the coastline may pass
 *
 * Distances are calculated in the users units using Euclidian geometry
 * unless a map projection and region (-R -J) are used.  Then, distances
 * are calculated using spherical geometry and converted to km, and any
 * distances given in options or via headers are assumed to be in km.
 *
 * Any one of these conditions may be negated for the opposite result
 * Both binary and ascii data files are accomodated
 *
 * Author:	Paul Wessel
 * Date:	25-AUG-1998
 * Version:	3.1
 *		3.2, 15-MAR-1999
 *		3.3, 13-APR-1999.	Added quick check for inside a square before checking
 *					for inside a circle (suggested by Joaquim Luis)
 *					Also added option -Cf for faster, approximate distances
 *		3.3.2 8-SEPT-1999.	Sorting -C points in x to speed up search as
 *				 	suggested by Joaquim Luis.
 *		3.3.4 17-FEB-2000.	Polygons containing either S or N pole will work correctly
 *		3.3.5 10-JUL-2000	Added plain -L for periodicity
 * Version:	3.4.1
 */
 
#include "gmt.h"

#define DEG_TO_KM (6371.0087714 * D2R)

char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"};

struct GMT_SHORE c;

struct LINES {
	double *lon, *lat;
	double dist;
	int np;
	int polar;
};

struct DATA {
	int      i;
	double   x;
	double   y;
	double   d;
} *data;

double dist_to_point (double lon, double lat, double *xp, double *yp, int np, int *id);
double GMT_cartesian_dist (double x0, double y0, double x1, double y1);
double GMT_flatearth_dist (double x0, double y0, double x1, double y1);
double GMT_km_dist (double x0, double y0, double x1, double y1);
int points_init (char *file, double **xp, double **yp, double **dp, double dist, BOOLEAN greenwich);
int lines_init (char *file, struct LINES **p, double dist, BOOLEAN greenwich);
int near_a_point_cart (double x, double y, double *xp, double *yp, double *dp, int np);
int near_a_point (double x, double y, double *xp, double *yp, double *dp, int np);
int near_a_line (double lon, double lat, struct LINES *p, int np);
int compare_x(const void *point_1, const void *point_2);

PFD distance_func;

main (int argc, char **argv)
{
	int i, j, k, fno, n_read = 0, n_pass = 0, n_files = 0, n_args, n_points, n_lines, n_pols;
	int n_fields, n_expected_fields, out_edge_in[5], ind, bin, last_bin = -1;
	int min_level = 0, max_level = MAX_LEVEL, np[2], base = 3, wd[2], id, this_node;
	int side, is_inside = 1, row, col;

	
	BOOLEAN error = FALSE, limit_line_distance = FALSE, limit_polygon = FALSE, limit_point_distance = FALSE, limit_point_distance_cart = FALSE;
	BOOLEAN want_inside[5] = {TRUE,TRUE,TRUE,TRUE,TRUE}, nofile = TRUE, done = FALSE, first = TRUE, inside, greenwich = FALSE;
	BOOLEAN flat_earth = FALSE, do_project = FALSE, binary = FALSE, limit_coast = FALSE, dry_wet_only;
	
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, xx, yy, min_area = 0.0;
	double *xp, *yp, *dp, pt_dist = 0.0, ln_dist = 0.0, *in, west_border, east_border;
	double xmin, xmax, ymin, ymax, lon, r, sa, ca;
	
	char buffer[BUFSIZ], *pt_file, *ln_file, *pol_file, *ptr, res = 'l';
	
	FILE *fp = NULL;
	
	struct LINES *pol, *line;
	struct POL *p[2];

	pt_file = ln_file = pol_file = CNULL;
	memset ((void *)out_edge_in, 0, (size_t)(5 * sizeof (int)));	/* Default for "wet" areas = 0 (outside) */
	out_edge_in[1] = out_edge_in[3] = 1;			/* Default for "dry" areas = 1 (inside) */

	argc = GMT_begin (argc, argv);
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'H':
				case 'J':
				case 'R':
				case 'V':
				case ':':
				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':
					j = sscanf (&argv[i][2], "%lf/%d/%d", &min_area, &min_level, &max_level);
					if (j == 1) min_level = 0, max_level = MAX_LEVEL;
					break;
				case 'C':
					limit_point_distance = TRUE;
					k = 2;
					if (argv[i][2] == 'f') flat_earth = TRUE, k = 3;
					for (j = k; argv[i][j] && argv[i][j] != '/'; j++);
					if (!argv[i][j]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -C:  Expects -C[f]dist/file\n", GMT_program);
						error++;
					}
					else {
						pt_file = &argv[i][j+1];
						pt_dist = atof (&argv[i][k]);
					}
					break;
				case 'D':
					res = argv[i][2];
					base = GMT_set_resolution (&res, 'D');
					break;
				case 'L':
					if (argv[i][2]) {	/* Set line options */
						limit_line_distance = TRUE;
						for (j = 2; argv[i][j] && argv[i][j] != '/'; j++);
							if (!argv[i][j]) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -L:  Expects -Ldist/file\n", GMT_program);
							error++;
						}
						else {
							ln_file = &argv[i][j+1];
							ln_dist = atof (&argv[i][2]);
						}
					}
					else	/* Flag to process geographic data */
						GMT_geographic_in = GMT_geographic_out = TRUE;
					break;
				case 'F':
					limit_polygon = TRUE;
					pol_file = &argv[i][2];
					break;
				case 'I':
					for (j = 2; argv[i][j]; j++) {
						switch (argv[i][j]) {
							case 'r':
								want_inside[0] = FALSE;
								break;
							case 'c':
								want_inside[1] = FALSE;
								break;
							case 'l':
								want_inside[2] = FALSE;
								break;
							case 'f':
								want_inside[3] = FALSE;
								break;
							case 's':
								want_inside[4] = FALSE;
								break;
							default:
								fprintf (stderr, "%s: GMT SYNTAX ERROR -I:  Expects -Icflrs\n", GMT_program);
								error++;
								break;
						}
					}
					break;
				case 'M':               /* Multiple line segments input */
					GMT_multisegment (&argv[i][2]);
					break;
				case 'N':
					limit_coast = TRUE;
					strcpy (buffer, &argv[i][2]);
					if (buffer[strlen(buffer)-1] == 'o') { /* Edge is considered outside */
						is_inside = 2;
						buffer[strlen(buffer)-1] = 0;
					}
					ptr = strtok (buffer, "/");
					j = 0;
					while (j < 5 && ptr) {
						switch (ptr[0]) {
							case 's':	/* Skip points in this level */
								out_edge_in[j] = 0;
								break;
							case 'k':
								out_edge_in[j] = 1;
								break;
							default:
								fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Bad modifier (use s or k)\n", GMT_program);
								error++;
						}
						ptr = strtok (CNULL, "/");
						j++;
					}
					if (!(j == 2 || j == 5)) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Specify 2 or 5 arguments\n", GMT_program);
						exit (EXIT_FAILURE);
					}
					dry_wet_only = (j == 2);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "gmtselect %s - Select data subsets based on multiple spatial criteria\n\n", GMT_VERSION);
		fprintf (stderr, "usage: gmtselect <infiles> [-A<min_area>[/<min_level>/<max_level>]] [-C[f]<dist>/<ptfile>] [-D<resolution>]\n");
		fprintf (stderr, "\t[-F<polygon>] [-J<parameters>] [-L[<dist>/<lfile>]] [-H[<nrec>]] [-I[cflrs]\n");
		fprintf (stderr, "\t[-M[<flag>]]  [-N<maskvalues>[o]] [-R<west/east/south/north>] [-V] [-:] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tinfiles (in ASCII or binary) have 2 or more columns with (x,y) or (y,x) in first columns.\n");
		fprintf (stderr, "\t  If no file(s) is given, standard input is read.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A coastline features smaller than <min_area> (in km^2) or of levels (0-4) outside the min-max levels\n");
		fprintf (stderr, "\t   will be skipped [0/4] (see pscoast for details)].  Ignored unless -N is set\n");
		fprintf (stderr, "\t-C pass locations that are within <dist> of any point in ASCII <ptfile>\n");
		fprintf (stderr, "\t   Give 0 distance if 3rd column of ptfile has individual distances.\n");
		fprintf (stderr, "\t   Distances are in users units [or km if -R -J are used].\n");
		fprintf (stderr, "\t   Prepend f to use approximate (flat earth) rather than exact great circle distances.\n");
		fprintf (stderr, "\t-D Choose one of the following resolutions: (Ignored unless -N is set)\n");
		fprintf (stderr, "\t   f - full resolution (may be very slow for large regions)\n");
		fprintf (stderr, "\t   h - high resolution (may be slow for large regions)\n");
		fprintf (stderr, "\t   i - intermediate resolution\n");
		fprintf (stderr, "\t   l - low resolution [Default]\n");
		fprintf (stderr, "\t   c - crude resolution, for tasks that need crude continent outlines only\n");
		fprintf (stderr, "\t-L with no arguments indicates that x is longitude and periodic in 360.\n");
		fprintf (stderr, "\t   With arguments, pass locations that are within <dist> of any line in ASCII <linefile>\n");
		fprintf (stderr, "\t   Give 0 distance if 2nd column of segment headers have individual distances.\n");
		fprintf (stderr, "\t   Distances are in users units [or km if -R -J are used].\n");
		fprintf (stderr, "\t-F pass locations that are inside the polygons in the ASCII <polygon> file\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I Used to reverse the tests, i.e. pass locations outside the region\n");
		fprintf (stderr, "\t   Supply a combination of cflr where each flag means:\n");
		fprintf (stderr, "\t   c will pass locations beyond the minimum distance to the points in -C\n");
		fprintf (stderr, "\t   f will pass locations outside the polygons in -F\n");
		fprintf (stderr, "\t   l will pass locations beyond the minimum distance to the lines in -L\n");
		fprintf (stderr, "\t   r will pass locations outside the region given in -R [and -J]\n");
		fprintf (stderr, "\t   s will pass locations that otherwise would be skipped in -N\n");
		GMT_explain_option ('J');
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N set if a point outside or inside a geographic feature should be s(kipped) or k(ept) (requires -R).\n");
		fprintf (stderr, "\t   Append o to let feature boundary be considered outside [Default is inside].\n");
		fprintf (stderr, "\t   Specify this information with s or k using 1 of 2 formats:\n");
		fprintf (stderr, "\t   -N<wet>/<dry>.\n");
		fprintf (stderr, "\t   -N<ocean>/<land>/<lake>/<island>/<pond>.\n");
		fprintf (stderr, "\t   k means keep and s means skip [Default is s/k/s/k/s (i.e., s/k)]\n");
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		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 (limit_coast && !project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option when -N is set\n", GMT_program);
		error++;
	}
	if (GMT_geographic_in && !project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option when plain -L is set\n", GMT_program);
		error++;
	}
	if (limit_coast && !MAPPING) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify a map projection when -N is set\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 (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 (error) exit (EXIT_FAILURE);

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

	if (!MAPPING && limit_point_distance) limit_point_distance_cart = TRUE;

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

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	if (project_info.region_supplied) {
		if (project_info.projection < 0) {	/* Supply dummy linear proj */
			project_info.projection = project_info.xyz_projection[0] = project_info.xyz_projection[1] = LINEAR;
			project_info.pars[0] = project_info.pars[1] = 1.0;
		}
		do_project = TRUE;
		if (MAPPING) {
			if (west < 0.0 && east < 0.0) {
				west += 360.0;
				east += 360.0;
			}
			greenwich = (west < 0.0 && east > 0.0);
		}
		GMT_map_setup (west, east, south, north);
	}

	if (limit_coast) {
		if (dry_wet_only) {
			out_edge_in[3] = out_edge_in[1];
			out_edge_in[2] = out_edge_in[4] = out_edge_in[0];
		}
		if (GMT_init_shore (res, &c, west, east, south, north))  {
			fprintf (stderr, "%s: %s resolution shoreline data base not installed\n", GMT_program, shore_resolution[base]);
			exit (EXIT_FAILURE);
		}
		west_border = floor (project_info.w / c.bsize) * c.bsize;
		east_border = ceil (project_info.e / c.bsize) * c.bsize;
		wd[0] = 1;	wd[1] = -1;
		np[0] = np[1] = 0;
	}

	binary = (GMT_io.binary[0] || GMT_io.binary[1]);

	distance_func = (PFD) ((MAPPING) ? ((flat_earth) ? GMT_flatearth_dist : GMT_km_dist) : GMT_cartesian_dist);

	if (limit_point_distance) { 
		n_points = points_init (pt_file, &xp, &yp, &dp, pt_dist, greenwich);
		
		/* Copy xp into struct data, sort, and copy back */
		
		data = (struct DATA *) GMT_memory (VNULL, (size_t)n_points, sizeof(struct DATA), GMT_program);

		for (k = 0; k < n_points; k++) {
			data[k].i = k;
			data[k].x = xp[k];
			data[k].y = yp[k];
			data[k].d = dp[k];
		}
		/* Sort on x */
		qsort ((void *)data, (size_t)n_points, sizeof (struct DATA), compare_x);
		for (k = 0; k < n_points; k++) {	/* Put back the new order */
			xp[k] = data[k].x;
			yp[k] = data[k].y;
			dp[k] = data[k].d;
		}
		GMT_free ((void *)data);
	}

	if (limit_line_distance)  n_lines  = lines_init (ln_file, &line, ln_dist, greenwich);
	if (limit_polygon) {
		n_pols = lines_init (pol_file, &pol, -9999.0, greenwich);
		if (do_project) {	/* Convert all the polygons */
			for (i = 0; i < n_pols; i++) {
				if (pol[i].polar) {	/* Polygon is a polar cap, use simple polar (r-theta) transformation */
					for (j = 0; j < pol[i].np; j++) {
						r = 90.0 - fabs (pol[i].lat[j]);
						sincos (pol[i].lon[j] * D2R, &sa, &ca);
						pol[i].lon[j] = r * ca;
						pol[i].lat[j] = r * sa;
					}
				}
				else {
					for (j = 0; j < pol[i].np; j++) {
						GMT_geo_to_xy (pol[i].lon[j], pol[i].lat[j], &xx, &yy);
						pol[i].lon[j] = xx;
						pol[i].lat[j] = yy;
					}
				}
			}
		}
	}

	/* Now we are ready to take on some input values */
	
	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
		
	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over input files, if any */
		if (!nofile && argv[fno][0] == '-') continue;
		
		if (nofile) {	/* Just read standard input */
			fp = GMT_stdin;
			done = TRUE;
			if (gmtdefs.verbose) fprintf (stderr, "%s: Reading from standard input\n", GMT_program);
#ifdef SET_IO_MODE
			GMT_setmode (0);
#endif
		}
		else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) {
			fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]);
			continue;
		}

		if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]);
		
		if (gmtdefs.io_header) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (buffer, BUFSIZ, fp);
				if (first) fprintf (GMT_stdout, "%s", buffer);
			}
			first = FALSE;
		}
		n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 2;
	
		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && ! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

			while (GMT_io.status & GMT_IO_SEGMENT_HEADER && !(GMT_io.status & GMT_IO_EOF)) {
				GMT_write_segmentheader (GMT_stdout, n_expected_fields);
				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++;

			lon = in[0];
			if (project_info.region_supplied) {
				inside = !GMT_map_outside (lon, in[1]);
				if (inside != want_inside[0]) continue;
			}

			if (limit_point_distance_cart) {
				inside = near_a_point_cart (lon, in[1], xp, yp, dp, n_points); 
				if (inside != want_inside[1]) continue;
			}
			else if (limit_point_distance) {
				inside = near_a_point (lon, in[1], xp, yp, dp, n_points);
				if (inside != want_inside[1]) continue;
			}

			if (limit_line_distance) {
				inside = near_a_line (lon, in[1], line, n_lines);
				if (inside != want_inside[2]) continue;
			}

			if (limit_polygon) {
				if (do_project)
					GMT_geo_to_xy (lon, in[1], &xx, &yy);
				else {
					xx = lon;
					yy = in[1];
				}
				i = 0;
				inside = FALSE;
				while (i < n_pols && !inside) {
					if (pol[i].polar == -1 && in[1] <= 0.0) {	/* Point on S hemisphere, check */
						r = 90.0 + in[1];
						sincos (lon * D2R, &sa, &ca);
						inside = (GMT_non_zero_winding (r * ca, r * sa, pol[i].lon, pol[i].lat, pol[i].np) >= is_inside);
					}
					else if (pol[i].polar == +1 && in[1] >= 0.0) {	/* Point on N hemisphere, check */
						r = 90.0 - in[1];
						sincos (lon * D2R, &sa, &ca);
						inside = (GMT_non_zero_winding (r * ca, r * sa, pol[i].lon, pol[i].lat, pol[i].np) >= is_inside);
					}
					else
						inside = (GMT_non_zero_winding (xx, yy, pol[i].lon, pol[i].lat, pol[i].np) >= is_inside);
					i++;
				}
				if (inside != want_inside[3]) continue;
			}

			if (limit_coast) {
				xx = lon;
				while (xx < 0.0) xx += 360.0;
				row = ((int)floor ((90.0 - in[1]) / c.bsize));
				if (row >= c.bin_ny) row = c.bin_ny - 1;	/* Presumably only kicks in for south pole */
				col = (int)floor (xx / c.bsize);
				bin = row * c.bin_nx + col;
				if (bin != last_bin) {	/* Do this upon entering new bin */
					ind = 0;
					while (ind < c.nb && c.bins[ind] != bin) ind++;	/* Set ind to right bin */
					if (ind == c.nb) continue;			/* Bin not among the chosen ones */
					last_bin = bin;
					GMT_get_shore_bin (ind, &c, min_area, min_level, max_level);

					/* Must use polygons.  Go in both directions to cover both land and sea */
					for (id = 0; id < 2; id++) {
						GMT_free_polygons (p[id], np[id]);
						if (np[id]) GMT_free ((void *)p[id]);
						np[id] = GMT_assemble_shore (&c, wd[id],  min_level, TRUE, greenwich, west_border, east_border, &p[id]);
						np[id] = GMT_prep_polygons (&p[id], np[id], greenwich, FALSE, 0.0, -1);
					}
				}

				if (c.ns == 0) {	/* No lines go through, check node level */
					this_node = MIN (MIN (c.node_level[0], c.node_level[1]) , MIN (c.node_level[2], c.node_level[3]));
				}
				else {
					this_node = 0;
					GMT_geo_to_xy (lon, in[1], &xx, &yy);
					for (id = 0; id < 2; id++) {
		
						for (k = 0; k < np[id]; k++) {
		
							if (p[id][k].n == 0) continue;

							/* Find min/max of polygon */
			
							xmin = xmax = p[id][k].lon[0];
							ymin = ymax = p[id][k].lat[0];

							for (i = 1; i < p[id][k].n; i++) {
								if (p[id][k].lon[i] < xmin) xmin = p[id][k].lon[i];
								if (p[id][k].lon[i] > xmax) xmax = p[id][k].lon[i];
								if (p[id][k].lat[i] < ymin) ymin = p[id][k].lat[i];
								if (p[id][k].lat[i] > ymax) ymax = p[id][k].lat[i];
							}
							
							if (yy < ymin || yy > ymax) continue;
							if (xx < xmin || xx > xmax) continue;

							/* Must compare with polygon */

							if ((side = GMT_non_zero_winding (xx, yy, p[id][k].lon, p[id][k].lat, p[id][k].n)) < is_inside) continue;	/* Outside polygon */
				
							/* Here, point is inside, we must assign value */
				
							if (p[id][k].level > this_node) this_node = p[id][k].level;
						}
					}
				}
				inside = out_edge_in[this_node];
				if (inside != want_inside[4]) continue;
			}

			/* Here, we have passed all test and the point is output */

			if (binary)
				GMT_output (GMT_stdout, n_expected_fields, in);
			else
				fprintf (GMT_stdout, "%s", GMT_io.current_record);
			n_pass++;
			
			if (gmtdefs.verbose && n_read%5000 == 0) fprintf (stderr, "%s: Read %d records\r", GMT_program, n_read);
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Read %d records, passed %d records\n", GMT_program, n_read, n_pass);

	if (limit_point_distance) {
		GMT_free ((void *)xp);
		GMT_free ((void *)yp);
		GMT_free ((void *)dp);
	}
	
	GMT_end (argc, argv);
}

int points_init (char *file, double **xp, double **yp, double **dp, double dist, BOOLEAN greenwich)
{
	FILE *fp;
	double *x, *y, *d, *in;
	int i = 0, n_alloc = GMT_CHUNK, n_fields, n_expected_fields = BUFSIZ;

	x = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	y = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);
	d = (double *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (double), GMT_program);

	if ((fp = GMT_fopen (file, "r")) == NULL) {
		fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	while ((n_fields = GMT_input_ascii (fp, &n_expected_fields, &in)) >= 0 && ! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

		while (GMT_io.status & GMT_IO_SEGMENT_HEADER && !(GMT_io.status & GMT_IO_EOF)) {
			n_fields = GMT_input_ascii (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 in file %s\n", GMT_program, n_fields,  n_expected_fields, i, file);
			exit (EXIT_FAILURE);
		}

		if (n_expected_fields < 2) {
			fprintf (stderr, "%s: Failure to read file %s near line %d\n", GMT_program, file, i);
			exit (EXIT_FAILURE);
		}
		
		x[i] = in[0];	y[i] = in[1];
		d[i] = (n_expected_fields >= 3 && dist == 0.0) ? in[2] : dist;
		if (MAPPING) {
			if (greenwich  && x[i] > 180.0) x[i] -= 360.0;
			if (!greenwich && x[i] < 0.0)   x[i] += 360.0;
		}

		i++;
		if (i == n_alloc) {
			n_alloc += GMT_CHUNK;
			x = (double *) GMT_memory ((void *)x, (size_t)n_alloc, sizeof (double), GMT_program);
			y = (double *) GMT_memory ((void *)y, (size_t)n_alloc, sizeof (double), GMT_program);
			d = (double *) GMT_memory ((void *)d, (size_t)n_alloc, sizeof (double), GMT_program);
		}

	}
	GMT_fclose (fp);
	
	if (i == 0) {
		fprintf (stderr, "%s: File %s is empty!\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}
	
	x = (double *) GMT_memory ((void *)x, (size_t)i, sizeof (double), GMT_program);
	y = (double *) GMT_memory ((void *)y, (size_t)i, sizeof (double), GMT_program);
	d = (double *) GMT_memory ((void *)d, (size_t)i, sizeof (double), GMT_program);
	*xp = x;
	*yp = y;
	*dp = d;

	return (i);
}

int lines_init (char *file, struct LINES **p, double dist, BOOLEAN greenwich)
{
	FILE *fp;
	struct LINES *e;
	int i = -1, j = 0, n, i_alloc = GMT_CHUNK, n_read = 0, j_alloc = GMT_CHUNK;
	int n_fields, n_expected_fields = BUFSIZ;
	BOOLEAN poly = FALSE, check_cap, save;
	double d, dlon, lon_sum, *in;

	if (fabs (dist + 9999.0) < GMT_CONV_LIMIT) poly = TRUE;	/* A polygon */
	check_cap = (poly && MAPPING);
	
	e = (struct LINES *) GMT_memory (VNULL, (size_t)i_alloc, sizeof (struct LINES), GMT_program);

	if ((fp = GMT_fopen (file, "r")) == NULL) {
		fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}

	save = GMT_io.multi_segments;	/* Must set this to TRUE temporarily */
	GMT_io.multi_segments = TRUE;
	
	n_fields = GMT_input_ascii (fp, &n_expected_fields, &in);
	if (GMT_io.status & GMT_IO_EOF) {
		fprintf (stderr, "%s: File %s is empty!\n", GMT_program, file);
		exit (EXIT_FAILURE);
	}
	if (!(GMT_io.status & GMT_IO_SEGMENT_HEADER)) {	/* Quick check up front that it is a -M kind of file */
		fprintf (stderr, "%s: Files for -F or -L must be in multisegment format!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	while (n_fields >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */
		while (GMT_io.status & GMT_IO_SEGMENT_HEADER && !(GMT_io.status & GMT_IO_EOF)) {
			/* To use different line-distances for each segment, place the distance in the segment header */
			i++;
			n_read++;
			n = sscanf (&GMT_io.segment_header[1], "%lg", &d);
			e[i].dist = (n == 1 && dist == 0.0) ? d : dist;
			j_alloc = GMT_CHUNK;
			j = 0;
			lon_sum = 0.0;
			n_fields = GMT_input_ascii (fp, &n_expected_fields, &in);
		}
		if ((GMT_io.status & GMT_IO_EOF)) continue;	/* At EOF */
		
		e[i].lon = (double *) GMT_memory (VNULL, (size_t)j_alloc, sizeof (double), GMT_program);
		e[i].lat = (double *) GMT_memory (VNULL, (size_t)j_alloc, sizeof (double), GMT_program);

		while (! (GMT_io.status & (GMT_IO_SEGMENT_HEADER | GMT_IO_EOF))) {	/* Keep going until FALSE or = 2 segment header */
			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, i);
				exit (EXIT_FAILURE);
			}
			
			n_read++;
			if (n_expected_fields < 2) {
				fprintf (stderr, "%s: Failure to read file %s near line %d\n", GMT_program, file, n_read);
				exit (EXIT_FAILURE);
			}
			e[i].lon[j] = in[0];
			e[i].lat[j] = in[1];
			if (MAPPING) {
				if (greenwich && e[i].lon[j] > 180.0) e[i].lon[j] -= 360.0;
				if (!greenwich && e[i].lon[j] < 0.0) e[i].lon[j] += 360.0;
			}

			if (check_cap && j > 0) {	/* Keep track of sum (dlon) */
				dlon = e[i].lon[j] - e[i].lon[j-1];
				if (fabs (dlon) > 180.0) dlon = copysign (360.0 - fabs (dlon), -dlon);
				lon_sum += dlon;
			}
			j++;
			if (j == (j_alloc-1)) {
				j_alloc += GMT_CHUNK;
				e[i].lon = (double *) GMT_memory ((void *)e[i].lon, (size_t)j_alloc, sizeof (double), GMT_program);
				e[i].lat = (double *) GMT_memory ((void *)e[i].lat, (size_t)j_alloc, sizeof (double), GMT_program);
			}
			n_fields = GMT_input_ascii (fp, &n_expected_fields, &in);
		}
		e[i].np = j;

		/* If dist = -9999.0 then file is a polygon and we must close it if needed */

		if (poly && !(e[i].lon[0] == e[i].lon[j-1] && e[i].lat[0] == e[i].lat[j-1])) {
			e[i].lon[j] = e[i].lon[0];
			e[i].lat[j] = e[i].lat[0];
			e[i].np++;
			if (check_cap) {	/* Keep track of sum (dlon) */
				dlon = e[i].lon[j] - e[i].lon[j-1];
				if (fabs (dlon) > 180.0) dlon = copysign (360.0 - fabs (dlon), -dlon);
				lon_sum += dlon;
			}
		}
		if (check_cap && fabs (fabs (lon_sum) - 360.0) < GMT_CONV_LIMIT) {	/* Contains a pole */
			e[i].polar = (e[i].lat[0] < 0.0) ? -1 : +1;	/* S or N pole */
		}
		
		if (i == (i_alloc-1)) {
			i_alloc += GMT_CHUNK;
			e = (struct LINES *) GMT_memory ((void *)e, (size_t)i_alloc, sizeof (struct LINES), GMT_program);
		}
	}
	GMT_fclose (fp);
	GMT_io.multi_segments = save;

	i++;
	e = (struct LINES *) GMT_memory ((void *)e, (size_t)i, sizeof (struct LINES), GMT_program);
	*p = e;

	return (i);
}

int near_a_line (double lon, double lat, struct LINES *p, int np)
{
	int i = 0, j0, j1, inside = FALSE;
	double d_min, dx, xc, yc, a, a_inv, d;

	while (i < np && !inside) {

		/* Find nearest point on this line */

		d_min = dist_to_point (lon, lat, p[i].lon, p[i].lat, p[i].np, &j0);
		if (j0 == 0)
			j1 = j0 + 1;
		else if (j0 == p[i].np-1)
			j1 = j0 - 1;
		else
			j1 = ((*distance_func) (lon, lat, p[i].lon[j0-1], p[i].lat[j0-1]) < (*distance_func) (lon, lat, p[i].lon[j0+1], p[i].lat[j0+1])) ? j0 - 1 : j0 + 1;

		dx = p[i].lon[j1] - p[i].lon[j0];
		if (dx == 0.0) {
			xc = p[i].lon[j0];
			yc = lat;
		}
		else {
			a = (p[i].lat[j1] - p[i].lat[j0]) / dx;
			a_inv = -1.0 / a;
			xc = (lat - p[i].lat[j0] + a_inv * lon - a * p[i].lon[j0]) / (a_inv - a);
			yc = lat - a_inv * (xc - lon);
		}

		if (xc >= MIN(p[i].lon[j0], p[i].lon[j1]) && xc <= MAX(p[i].lon[j0], p[i].lon[j1]) &&
			yc >= MIN(p[i].lat[j0], p[i].lat[j1]) && yc <= MAX(p[i].lat[j0], p[i].lat[j1]))
			d = (*distance_func) (lon, lat, xc, yc);
		else
			d = d_min;
		inside = (d <= p[i].dist);
		i++;
	}
	return (inside);
}

double dist_to_point (double lon, double lat, double *xp, double *yp, int np, int *id)
{
	int i;
	double d, d_min;

	d_min = 1.0e100;
	for (i = 0; i < np; i++) {
		d = (*distance_func) (lon, lat, xp[i], yp[i]);
		if (d < d_min) {
			d_min = d;
			*id = i;
		}
	}
	return (d_min);
}

int near_a_point (double x, double y, double *xp, double *yp, double *dp, int np)
{
	int i = 0, inside = FALSE;
	double d;

	while (i < np && !inside) {
		d = (*distance_func) (x, y, xp[i], yp[i]);
		inside = (d <= dp[i]);
		i++;
	}
	return (inside);
}

int near_a_point_cart (double x, double y, double *xp, double *yp, double *dp, int np)
{
	int i = 0, inside = FALSE;
	double d;

	if ((x < (xp[0] - dp[0])) || (x > (xp[np-1]) + dp[np-1])) return (inside);
	while (i < np && !inside) {
		if (fabs (x - xp[i]) <= dp[i])
			if (fabs (y - yp[i]) <= dp[i]) {
				d = (*distance_func) (x, y, xp[i], yp[i]);
				inside = (d <= dp[i]);
			}
		i++;
	}
	return (inside);
}

double GMT_cartesian_dist (double x0, double y0, double x1, double y1)
{
	/* Calculates the good-old straight line distance in users units */

	return (hypot ( (x1 - x0), (y1 - y0)));
}

double GMT_flatearth_dist (double x0, double y0, double x1, double y1)
{
	/* Calculates the approximate flat earth distance in km */

	return (hypot ( (x1 - x0) * cosd (0.5 * (y1 + y0)), (y1 - y0)) * DEG_TO_KM);
}

double GMT_km_dist (double x0, double y0, double x1, double y1)
{
	/* Calculates the grdat circle distance in km */

	return (GMT_great_circle_dist (x0, y0, x1, y1) * DEG_TO_KM);
}

int compare_x (const void *point_1, const void *point_2)
{
	struct DATA *p1, *p2;
	
	p1 = (struct DATA *)point_1;
	p2 = (struct DATA *)point_2;

	if (p1->x < p2->x)
		return (-1);
	else if (p1->x > p2->x)
		return (1);
	else
		return (0);
}