/*--------------------------------------------------------------------
 *    The GMT-system:	@(#)nearneighbor.c	2.70  10/29/99
 *
 *	Copyright (c) 1991-1999 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: www.soest.hawaii.edu/gmt
 *--------------------------------------------------------------------*/
/*
 * Based on a specified grid size, nearneighbor reads an xyz file and
 * determines the nearest points to each node in sectors.  The default
 * looks for the nearest point for each quadrant.  The points must also
 * be within a maximum search-radius from the node.  For the nodes that
 * have a full set of nearest neighbors, a weighted average value is
 * computed.  New feature is full support for boundary conditions so
 * that geographic or periodic conditions are explicitly dealt with
 * in the sense that a data point may wrap around to serve as a
 * constraint on the other side of the periodic boundary.
 *
 * Author:	Paul Wessel
 * Date:	02-JUN-1999
 * Version:	3.3
 */
 
#include "gmt.h"
#include "gmt_boundcond.h"

float *grd;

struct NODE {	/* Structure with point id and distance pairs for all sectors */
	float *distance;	/* Distance of nearest datapoint to this node per sector */
	int *datum;		/* Point id of this data point */
} **grid_node;

struct POINT {	/* Structure with input data constraints */
	float x, y, z, w;
} *point;

struct NODE *add_new_node(int n);
void assign_node (struct NODE **node, int n_sector, int sector, double distance, int id);

main (int argc, char **argv)
{
	int i, j, k, ij, i0, j0, *di, dj, n_sectors = 4, sector, n, n_alloc = 5 * GMT_CHUNK, n_fields, nx_2;
	int ix, iy, n_set, n_almost, n_none, n_files = 0, n_args, fno, one_or_zero, n_expected_fields, pad[4], distance_flag = 0;
	int max_di, actual_max_di, ii, jj, n_req;

	BOOLEAN go, error = FALSE, done = FALSE, first = TRUE, nofile = TRUE;
	BOOLEAN set_empty = FALSE, weighted = FALSE, skip, wrap_180;

	double radius = 0.0, weight, weight_sum, *x0, *y0, dx, dy, delta, distance, factor;
	double *in, *shrink, km_pr_deg, x_left, x_right, y_top, y_bottom, offset, xinc2, yinc2, idx, idy;
	double half_y_width, y_width, half_x_width, x_width;

	float empty = 0.0;
	char *outfile = CNULL, line[BUFSIZ];

	FILE *fp = NULL;

	struct GRD_HEADER header;
	
	struct GMT_EDGEINFO edgeinfo;

	argc = GMT_begin (argc, argv);
	
	GMT_boundcond_init (&edgeinfo);

	GMT_grd_init (&header, argc, argv, FALSE);

	pad[0] = pad[1] = pad[2] = pad[3] = 0;
	
	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':	/* Input triplets [quadruplets] are binary, not ascii */
					error += GMT_io_selection (&argv[i][2]);
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
					break;
				case 'E':
					if (!argv[i][2]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Must specify value or NaN\n", GMT_program);
						error++;
					}
					else
						empty = (argv[i][2] == 'N' || argv[i][2] == 'n') ? GMT_f_NaN : (float)atof (&argv[i][2]);	
					set_empty = TRUE;
					break;
				case 'F':
					header.node_offset = TRUE;
					break;
				case 'G':
					outfile = &argv[i][2];
					break;
				case 'L':
					error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]);
					break;
				case 'N':
					n_sectors = atoi (&argv[i][2]);
					break;
				case 'S':
					GMT_getinc (&argv[i][2], &radius, &radius);
					if (argv[i][strlen(argv[i])-1] == 'k') distance_flag = 1;
					if (argv[i][strlen(argv[i])-1] == 'K') distance_flag = 2;
					break;
				case 'W':
					weighted = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else 
			n_files++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "nearneighbor %s - A \"Nearest neighbor\" gridding algorithm\n\n", GMT_VERSION);
		fprintf(stderr, "usage: nearneighbor [xyzfile(s)] -G<out_grdfile> -I<dx>[m|c][/<dy>[m|c]]\n");
		fprintf(stderr, "	-N<sectors> -R<west/east/south/north> -S<radius>[m|c|k|K] [-E<empty>] [-F]\n");
		fprintf(stderr, "	[-H ] [-L<flags>] [-V ] [-W] [-:] [-bi[s][<n>]]\n\n");
		if (GMT_quick) exit (EXIT_FAILURE);
		fprintf(stderr, "	-G name of output grid.\n");
		fprintf(stderr, "	-I sets the grid spacing for the grid.  Append m for minutes, c for seconds.\n");
		fprintf(stderr, "	-N sets number of sectors. Default is quadrant search [4].\n");
		GMT_explain_option ('R');
		fprintf(stderr, "	-S sets search radius in -R, -I units; append m or c for minutes or seconds.\n");
		fprintf(stderr, "	   Append k for km (implies -R,-I in degrees), use flat Earth approximation.\n");
		fprintf(stderr, "	   Append K for km (implies -R,-I in degrees), use great circle distances.\n");
		fprintf(stderr, "	   \n");
		fprintf(stderr, "\n\tOPTIONS:\n");
		fprintf(stderr, "	-E value to use for empty nodes [Default is NaN].\n");
		fprintf(stderr, "	-F Force pixel registration [Default is gridline registration].\n");
		GMT_explain_option ('H');
		fprintf(stderr, "	-L sets boundary conditions.  <flags> can be either\n");
		fprintf(stderr, "	   g for geographic boundary conditions, or one or both of\n");
		fprintf(stderr, "	   x for periodic boundary conditions on x\n");
		fprintf(stderr, "	   y for periodic boundary conditions on y\n");
		GMT_explain_option ('V');
		fprintf(stderr, "	-W input file has observation weights in 4th column.\n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf(stderr, "	   Default is 3 (or 4 if -W is set) 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 (!outfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR option -G:  Must specify output file\n", GMT_program);
		error++;
	}
	if (n_sectors <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Must specify a positive number of sectors\n", GMT_program);
		error++;
	}
	if (radius <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Must specify a positive search radius\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++;
	}	
	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++;
	}
	n_req = (weighted) ? 4 : 3;
	if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = n_req;
	if (GMT_io.binary[0] && n_req > GMT_io.ncol[0]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  binary input data must have at least %d columns\n", GMT_program, n_req);
		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);

	n_expected_fields = (GMT_io.binary[0]) ? GMT_io.ncol[0] : BUFSIZ;
	
        ix = (gmtdefs.xy_toggle);       iy = 1 - ix;
        
        if (n_files > 0)
        	nofile = FALSE;
        else
        	n_files = 1;
        n_args = (argc > 1) ? argc : 2;
        
	if (header.node_offset) {
		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;
	}
	
	idx = 1.0 / header.x_inc;
	idy = 1.0 / header.y_inc;

	header.nx = irint ( (header.x_max - header.x_min) * idx) + one_or_zero;
	header.ny = irint ( (header.y_max - header.y_min) * idy) + one_or_zero;
	
	GMT_boundcond_param_prep (&header, &edgeinfo);

	if (gmtdefs.verbose) fprintf (stderr, "%s:  Grid dimensions are nx = %d, ny = %d\n", GMT_program,
		header.nx, header.ny);

	grid_node = (struct NODE **) GMT_memory (VNULL, (size_t)(header.nx * header.ny), sizeof (struct NODE *), GMT_program);
	point = (struct POINT *) GMT_memory (VNULL, (size_t)n_alloc, sizeof (struct POINT), GMT_program);
	
	di = (int *) GMT_memory (VNULL, (size_t)header.ny, sizeof (int), GMT_program);
	shrink = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof (double), GMT_program);

	x0 = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof (double), GMT_program);
	y0 = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof (double), GMT_program);
	for (i = 0; i < header.nx; i++) x0[i] = header.x_min + i * header.x_inc + xinc2;
	for (j = 0; j < header.ny; j++) y0[j] = header.y_max - j * header.y_inc - yinc2;
	if (distance_flag) {	/* Input data is geographical */
		km_pr_deg = 0.001 * 2.0 * M_PI * gmtdefs.ellipse[gmtdefs.ellipsoid].eq_radius / 360.0;
		max_di = (int) (ceil (header.nx / 2.0) + 0.1);
		actual_max_di = 0;
		for (j = 0; j < header.ny; j++) {
			shrink[j] = cosd (y0[j]);
			di[j] = (fabs (y0[j]) == 90.0) ? max_di : (int)(ceil (radius / (km_pr_deg * header.x_inc * shrink[j])) + 0.1);
			if (di[j] > max_di) di[j] = max_di;
			if (di[j] > actual_max_di) actual_max_di = di[j];
		}
		dj = (int) (ceil (radius / (km_pr_deg * header.y_inc)) + 0.1);
	}
	else {	/* Plain Cartesian data */
		max_di = (int) (ceil (radius * idx) + 0.1);
		for (j = 0; j < header.ny; j++) di[j] = max_di;
		dj = (int) (ceil (radius * idy) + 0.1);
		actual_max_di = max_di;
	}
	factor = n_sectors / (2.0 * M_PI);

	x_left = header.x_min - actual_max_di * header.x_inc;	x_right = header.x_max + actual_max_di * header.x_inc;
	y_top = header.y_max + dj * header.y_inc;	y_bottom = header.y_min - dj * header.y_inc;
	x_width = header.x_max - header.x_min;		y_width = header.y_max - header.y_min;
	half_x_width = 0.5 * x_width;			half_y_width = 0.5 * y_width;
	nx_2 = edgeinfo.nxp / 2;
	n = 0;
	
	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;
#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++) {
				fgets (line, BUFSIZ, fp);
				if (first) printf ("%s", line);
			}
			first = FALSE;
		}

		n_fields = GMT_input (fp, &n_expected_fields, &in);
			
		while (! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */
	
			skip = 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);
				exit (EXIT_FAILURE);
			}

			if (in[ix] < x_left || in[ix] > x_right) skip = TRUE;
			if (in[iy] < y_bottom || in[iy] > y_top) skip = TRUE;
		
			if (!skip) {
				point[n].x = (float)in[ix];
				point[n].y = (float)in[iy];
				point[n].z = (float)in[2];
				if (weighted) point[n].w = (float)in[3];
		
				/* Find indeces of the node closest to this data point */

				i0 = (int)floor (((in[ix] - header.x_min) * idx) + offset);
				j0 = (int)floor (((header.y_max - in[iy]) * idy) + offset);
		
				/* Loop over all nodes within radius of this node */

				for (j = j0 - dj; j <= (j0 + dj); j++) {
					
					wrap_180 = FALSE;

					if (j < 0) {	/* Depending on BC's we wrap around or skip */
						if (edgeinfo.gn) {	/* N Polar condition */
							jj = abs (j) - header.node_offset;
							wrap_180 = TRUE;
						}
						else if (edgeinfo.nyp) {	/* Periodic in y */
							jj = j + edgeinfo.nyp;
						}
						else
							continue;
					}
					else if (j >= header.ny) {	/* Depending on BC's we wrap around or skip */
						if (edgeinfo.gs) {	/* S Polar condition */
							jj = j - 2 + header.node_offset;
							wrap_180 = TRUE;
						}
						else if (edgeinfo.nyp) {	/* Periodic in y */
							jj = j - edgeinfo.nyp;
						}
						else
							continue;
					}
					else
						jj = j;

					for (i = i0 - di[jj]; i <= (i0 + di[jj]); i++) {

						if (i < 0) {	/* Depending on BC's we wrap around or skip */
							if (edgeinfo.nxp) {	/* Periodic in x */
								ii = abs (i) - header.node_offset;
							}
							else
								continue;
						}
						else if (i >= header.nx) {	/* Depending on BC's we wrap around or skip */
							if (edgeinfo.nxp) {	/* Periodic in x */
								ii = i - edgeinfo.nxp;
							}
							else
								continue;
						}
						else
							ii = i;

						if (wrap_180) ii = (ii + nx_2) % edgeinfo.nxp;
						k = jj * header.nx + ii;
						dx = in[ix] - x0[ii];	dy = in[iy] - y0[jj];

						/* Check for wrap-around in x or y */
						if (edgeinfo.nxp && fabs (dx) > half_x_width) dx -= copysign (x_width, dx);
						if (edgeinfo.nyp && fabs (dy) > half_y_width) dy -= copysign (y_width, dy);

						switch (distance_flag) {	/* Take different action depending on how we want distances calculated */

							case 0:		/* Cartesian distance */

								distance = hypot (dx, dy);
								break;
							case 1:		/* Flat Earth Approximation */

								distance = km_pr_deg * hypot (dx * shrink[jj], dy);
								break;
							case 2:		/* Full spherical calculation */

								distance = km_pr_deg * GMT_great_circle_dist (x0[ii], y0[jj], in[ix], in[iy]);
								break;
							default:
								break;
						}

						if (distance > radius) continue;

						sector = ((int)((d_atan2 (dy, dx) + M_PI) * factor)) % n_sectors;
						assign_node (&grid_node[k], n_sectors, sector, distance, n);

						if (edgeinfo.nxp && ii == 0) {	/* Must replicate to redundant column */
							assign_node (&grid_node[k+header.nx-1], n_sectors, sector, distance, n);
							if (edgeinfo.nyp && jj == 0)	/* Must replicate to redundant row */
								assign_node (&grid_node[header.nx-1], n_sectors, sector, distance, n);
						}
						if (edgeinfo.nyp && jj == 0)	/* Must replicate to redundant row */
							assign_node (&grid_node[ii], n_sectors, sector, distance, n);
					}
				}
				n++;
				if (n == n_alloc) {
					n_alloc += GMT_CHUNK;
					point = (struct POINT *) GMT_memory ((void *)point, (size_t)n_alloc, sizeof (struct POINT), GMT_program);
				}
			}
			n_fields = GMT_input (fp, &n_expected_fields, &in);
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	point = (struct POINT *) GMT_memory ((void *)point, (size_t)n, sizeof (struct POINT), GMT_program);
	grd = (float *) GMT_memory (VNULL, (size_t)(header.nx * header.ny), sizeof (float), GMT_program);

	/* Compute weighted averages based on the nearest neighbors */
	
	n_set = n_almost = n_none = 0;

	if (!set_empty) empty = GMT_f_NaN;
	
	for (j = ij = 0; j < header.ny; j++) {
		for (i = 0; i < header.nx; i++, ij++) {
			grd[ij] = empty;

			if (!grid_node[ij]) {	/* No nearest neighbors */
				n_none++;
				continue;
			}

			for (k = 0, go = TRUE; go && k < n_sectors; k++) if (grid_node[ij]->datum[k] < 0) go = FALSE;
			if (!go) { 	/* Not full set of neighbors */
				n_almost++;
				continue;
			}

			n_set++;
			weight_sum = grd[ij] = 0.0;	/* Replace the empty so that we may compute a sum */
			for (k = 0; k < n_sectors; k++) {
				delta = 3.0 * grid_node[ij]->distance[k] / radius;
				weight = 1.0 / (1.0 + delta * delta);	/* This is distance weight */
				if (weighted) weight *= point[grid_node[ij]->datum[k]].w;	/* This is observation weight */
				grd[ij] += (float)(weight * point[grid_node[ij]->datum[k]].z);
				weight_sum += weight;
			}
			grd[ij] /= (float)weight_sum;
		}
	}
	
	if (GMT_write_grd (outfile, &header, grd, 0.0, 0.0, 0.0, 0.0, pad, FALSE)) {
		fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile);
		exit (EXIT_FAILURE);
	}
	
	if (gmtdefs.verbose) {
		sprintf (line, "%s)\n\0", gmtdefs.d_format);
		fprintf (stderr, "%s: %d nodes were assigned an average value\n", GMT_program, n_set);
		fprintf (stderr, "%s: %d nodes failed sector criteria and %d nodes had no neighbor points (all set to ", GMT_program, n_almost, n_none);
		(GMT_is_dnan (empty)) ? fprintf (stderr, "NaN)\n") : fprintf (stderr, line, empty);
	}

	GMT_free ((void *)grd);
	GMT_free ((void *)point);
	GMT_free ((void *)grid_node);
	GMT_free ((void *)shrink);
	GMT_free ((void *)di);
	GMT_free ((void *)x0);
	GMT_free ((void *)y0);
	
	GMT_end (argc, argv);
}

struct NODE *add_new_node(int n)
{
	struct NODE *new;
	
	new = (struct NODE *) GMT_memory (VNULL, (size_t)1, sizeof (struct NODE), GMT_program);
	new->distance = (float *) GMT_memory (VNULL, (size_t)n, sizeof (float), GMT_program);
	new->datum = (int *) GMT_memory (VNULL, (size_t)n, sizeof (int), GMT_program);
	while (n > 0) new->datum[--n] = -1;
	
	return (new);
}

void assign_node (struct NODE **node, int n_sector, int sector, double distance, int id)
{
	/* Allocates node space if not already used and updates the value if closer to node */

	if (!(*node)) *node = add_new_node (n_sector);
	if ((*node)->datum[sector] == -1 || (*node)->distance[sector] > distance) {
		(*node)->distance[sector] = (float)distance;
		(*node)->datum[sector] = id;
	}
}