/*--------------------------------------------------------------------
 *	$Id: grdmask.c,v 1.3.4.2 2002/02/27 17:41:10 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
 *--------------------------------------------------------------------*/
/*
 * grdmask defines a grid based on region and xinc/yinc values,
 * reads xy polygon files, and sets the grid nodes inside, on the
 * boundary, and outside of the polygons to the user-defined values
 * <in>, <on>, and <out>.  These may be any number, including NaN.
 *
 * Author:	Walter. H. F. Smith
 * Date:	23-May-1991
 * Version:	2.0
 * Modified:	PW: 12-MAY-1998, to GMT 3.1
 * 		PW: 18-FEB-2000, to GMT 3.3.4: Handle polarcap polygons
 *				and for -S check if points are inside -R.
 *				Also added -L for geographic data
 * 		PW: 13-JUL-2000, 3.3.5
 * Version:	3.4.1
 */
 
#include "gmt.h"

float *data;

main (int argc, char **argv)
{

	int	i, j, k, ij, nm, n_path, side, fno, n_files = 0, n_alloc = GMT_CHUNK, n_args;
	int	one_or_zero, di, dj, i0, j0, n_fields, n_read, n_expected_fields, polar;
	
	BOOLEAN	error = FALSE, done, nofile = TRUE, pixel = FALSE;
	BOOLEAN	do_coverage = FALSE, map_units = FALSE;
	
	double	*x, *y, xx, yy, xmin, xmax, ymin, ymax, offset, dx, dy, r, ca, sa, lon_sum, dlon;
	double xinc2, yinc2, distance, radius = 0.0, shrink, km_pr_deg, idx, idy, x0, y0, *in;

	float out_edge_in[3];
	
	char *maskfile = CNULL, line[BUFSIZ], *ptr;
	
	FILE *fp;
	
	struct GRD_HEADER header;
	
	argc = GMT_begin (argc, argv);
		
	GMT_grd_init (&header, argc, argv, FALSE);
	
	out_edge_in[2] = 1.0;			/* Default inside value */
	out_edge_in[0] = out_edge_in[1] = 0.0;	/* Default outside and edge value */

	/* Check command line arguments */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
			
				/* Common parameters */
			
				case 'H':
				case 'R':
				case 'V':
				case ':':
				case '\0':
					error += GMT_get_common_args (argv[i], &header.x_min, &header.x_max, &header.y_min, &header.y_max);
					break;
				
				/* Supplemental parameters */
				
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'N':
					strcpy (line, &argv[i][2]);
					ptr = strtok (line, "/");
					for (j = 0; j < 3; j++) {
						out_edge_in[j] = (ptr[0] == 'N' || ptr[0] == 'n') ? GMT_f_NaN : (float)atof (ptr);
						ptr = strtok (CNULL, "/");
					}
					break;
				case 'F':
					pixel = TRUE;
					break;
				case 'G':
					maskfile = &argv[i][2];
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
					break;
				case 'L':
					GMT_geographic_in = TRUE;
					break;
				case 'M':
					GMT_multisegment (&argv[i][2]);
					break;
				case 'S':
					GMT_getinc (&argv[i][2], &radius, &radius);
					do_coverage = TRUE;
					map_units = (argv[i][strlen(argv[i])-1] == 'k');
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdmask %s - Create mask grdfile from polygons or point coverage\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdmask [<xyfiles>] -G<mask_grd_file> -I<xinc[m|c]>[/<yinc>[m|c]]\n");
		fprintf (stderr, "\t-R<west/east/south/north> [-F] [-H[<nrec>]] [-L] [-N<out>/<edge>/<in>] [-M[<flag>]] [-S<radius>[k] [-V] [-:] [-bi[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\txyfiles is one or more polygon [or point] files\n");
		fprintf (stderr, "\t-G Specify file name for output mask grd file.\n");
		fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
		fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or c] to indicate minutes [or seconds].\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-F Force pixel registration for output grid [Default is gridline orientation]\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-L means that x is longitude, i.e. assumed to be periodic in 360\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N sets values to use if point is outside, on the path, or inside.\n");
		fprintf (stderr, "\t   NaN is a valid entry.  Default values are 0/0/1.\n");
		fprintf (stderr, "\t-S Input data is points and the mask nodes are set to <in> or <out> depending on\n");
		fprintf (stderr, "\t    whether they are inside a circle of specified radius [0] from the nearest data point.\n");
		fprintf (stderr, "\t    Append k to indicated km (implying -R -I are in degrees)\n");
		fprintf (stderr, "\t    [Default is to treat xyfiles as polygons and use inside/outside searching].\n");
		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 ('.');
		exit (EXIT_FAILURE);
	}

	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (header.x_inc <= 0.0 || header.y_inc <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", GMT_program);
		error = TRUE;
	}
	if (!maskfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", GMT_program);
		error = TRUE;
	}
	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 (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]]);
	}

	GMT_grd_RI_verify (&header, 1);

	if (pixel) {
		one_or_zero = 0;
		offset = 0.0;
		xinc2 = 0.5 * header.x_inc;
		yinc2 = 0.5 * header.y_inc;
	}
	else {
		one_or_zero = 1;
		offset = 0.5;
		xinc2 = yinc2 = 0.0;
	}
	header.nx = irint ((header.x_max - header.x_min) / header.x_inc) + one_or_zero;
	header.ny = irint ((header.y_max - header.y_min) / header.y_inc) + one_or_zero;
	header.node_offset = pixel;
	idx = 1.0 / header.x_inc;
	idy = 1.0 / header.y_inc;
	nm = header.nx * header.ny;
	data = (float *) GMT_memory (VNULL, (size_t)nm, sizeof(float), GMT_program);
	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);

	sprintf (line, "%s\n", gmtdefs.d_format);
	if (gmtdefs.verbose) {
		fprintf (stderr, "%s: Nodes completely outside the polygons will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[0])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[0]);
		fprintf (stderr, "%s: Nodes completely inside the polygons will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[2])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[2]);
		fprintf (stderr, "%s: Nodes on the polygons boundary will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[1])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[1]);
	}
	
	if (map_units) {
		km_pr_deg = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0;
		shrink = cos (0.5 * (header.y_min + header.y_max) * D2R);
		di = (int)ceil (radius / (km_pr_deg * header.x_inc * shrink));
		dj = (int)ceil (radius / (km_pr_deg * header.y_inc));
	}
	else {
		di = (int)ceil (radius * idx);
		dj = (int)ceil (radius * idy);
	}
	/* Initialize all nodes to the 'outside' value */
	
	for (ij = 0; ij < nm; ij++) data[ij] = out_edge_in[0];
	
	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
	
	done = FALSE;
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 2;
	
	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over all input files */
		if (!nofile && argv[fno][0] == '-') continue;
		if (nofile) {
			fp = GMT_stdin;
			done = TRUE;
#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 (line, BUFSIZ, fp);
		
		n_read = 0;
		n_fields = GMT_input (fp,  &n_expected_fields, &in);
		while (! (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 (fp, &n_expected_fields, &in);
			if ((GMT_io.status & GMT_IO_EOF)) continue;	/* At EOF */

			n_path = polar = 0;
			xmin = ymin = 1.0e+100;
			xmax = ymax = -1.0e+100;
			lon_sum = 0.0;
	
			while (!(GMT_io.status & (GMT_IO_SEGMENT_HEADER | GMT_IO_EOF))) {	/* Keep going until FALSE  */
				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);
				}
				x[n_path] = in[0];
				y[n_path] = in[1];
	
				if (x[n_path] > xmax) xmax = x[n_path];
				if (x[n_path] < xmin) xmin = x[n_path];
				if (y[n_path] > ymax) ymax = y[n_path];
				if (y[n_path] < ymin) ymin = y[n_path];
				if (n_path > 0) {	/* Do longitude-difference check sum */
					dlon = x[n_path] - x[n_path-1];
					if (fabs (dlon) > 180.0) dlon = copysign (360.0 - fabs (dlon), -dlon);
					lon_sum += dlon;
				}
				n_path++;
		
				if (n_path == (n_alloc-1)) {	/* n_alloc-1 since we may need 1 more to close polygon */
					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);
				}
				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}

			if (do_coverage) {	/* Assign 'inside' to nodes within given distance of data constrains */
				for (k = 0; k < n_path; k++) {

					if (y[k] < header.y_min || y[k] > header.y_max) continue;
					if (GMT_geographic_in) {	/* Make sure we are not off by 360 somewhere */
						x[k] -= 360.0;
						while (x[k] < header.x_min) x[k] += 360.0;
						if (x[k] > header.x_max) continue;
					}
					else if (x[k] < header.x_min || x[k] > header.x_max) continue;
					
					/* OK, this point is within bounds */
					
					i0 = (int)floor (((x[k] - header.x_min) * idx) + offset);
					j0 = (int)floor (((header.y_max - y[k]) * idy) + offset);
					data[j0 * header.nx + i0] = out_edge_in[2];	/* This is the nearest node */
					if (radius == 0.0) continue;

					for (j = j0 - dj; j <= (j0 + dj); j++) {
						if (j < 0 || j >= header.ny) continue;
						for (i = i0 - di; i <= (i0 + di); i++) {
							if (i < 0 || i >= header.nx) continue;
							ij = j * header.nx + i;
							x0 = header.x_min + i * header.x_inc + xinc2;
							y0 = header.y_max - j * header.y_inc - yinc2;
							dx = x[k] - x0;	dy = y[k] - y0;
							distance = (map_units) ? km_pr_deg * hypot (dx * shrink, dy) : hypot (dx, dy);
							if (distance > radius) continue;
							data[ij] = out_edge_in[2];	/* The in value */
						}
					}
				}
			}
			else {	/* assign 'inside' to nodes if they are inside given polygon */
				if (x[n_path-1] != x[0] || y[n_path-1] != y[0]) {
					x[n_path] = x[0];
					y[n_path] = y[0];
					dlon = x[n_path] - x[n_path-1];
					n_path++;
					if (fabs (dlon) > 180.0) dlon = copysign (360.0 - fabs (dlon), -dlon);
					lon_sum += dlon;
				}
				if (fabs (fabs (lon_sum) - 360.0) < GMT_CONV_LIMIT) {	/* Contains a pole, convert to polar x,y */
					polar = (y[0] < 0.0) ? -1 : +1;	/* S or N pole */
					for (j = 0; j < n_path; j++) {
						r = 90.0 - fabs (y[j]);
						sincos (x[j] * D2R, &sa, &ca);
						x[j] = r * ca;
						y[j] = r * sa;
					}
				}

				for (j = 0; j < header.ny; j++) {
		
					yy = header.y_max - j * header.y_inc;
					if (pixel) yy -= yinc2;
					if (!polar && (yy < ymin || yy > ymax)) continue;	/* Point outside, no need to assign value */
			
					for (i = 0; i < header.nx; i++) {
						xx = header.x_min + i * header.x_inc;
						if (pixel) xx += xinc2;
						if (GMT_geographic_in) {	/* Make sure we are not off by 360 somewhere */
							xx -= 360.0;
							while (xx < xmin) xx += 360.0;
						}
						if (!polar && (xx < xmin || xx > xmax)) continue;	/* Point outside, no need to assign value */
				
						if (polar == -1 && yy <= 0.0) {	/* Point on S hemisphere, check */
							r = 90.0 + yy;
							sincos (xx * D2R, &sa, &ca);
							side = GMT_non_zero_winding (r * ca, r * sa, x, y, n_path);
						}
						else if (polar == +1 && yy >= 0.0) {	/* Point on N hemisphere, check */
							r = 90.0 - yy;
							sincos (xx * D2R, &sa, &ca);
							side = GMT_non_zero_winding (r * ca, r * sa, x, y, n_path);
						}
						else if (polar)
							side = 0;
						else
							side = GMT_non_zero_winding (xx, yy, x, y, n_path);
							
						if (side == 0) continue;	/* Outside */
				
						/* Here, point is inside or on edge, we must assign value */
				
						data[j*header.nx+i] = out_edge_in[side];
					}
					if (gmtdefs.verbose) fprintf (stderr, "%s: scanning row %5.5d\r", GMT_program, j);
				}
			}
		}
		if (gmtdefs.verbose) fprintf (stderr, "\n");
		
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	if (GMT_write_grd (maskfile, &header, data, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error writing file %s\n", GMT_program, maskfile);
		exit (EXIT_FAILURE);
	}

	GMT_free ((void *)data);
	GMT_free ((void *)x);
	GMT_free ((void *)y);
	
	GMT_end (argc, argv);
}