/*--------------------------------------------------------------------
 *	$Id: grdsample.c,v 1.35 2006/04/13 08:34:41 pwessel Exp $
 *
 *	Copyright (c) 1991-2006 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
 *--------------------------------------------------------------------*/
/*
 * grdsample reads a grdfile and evaluates the grid at new grid positions
 * specified by new dx/dy values using a 2-D Taylor expansion of order 3.
 * In order to evaluate derivatives along the edges of the surface, I assume 
 * natural bi-cubic spline conditions, i.e. both the second and third normal 
 * derivatives are zero, and that the dxdy derivative in the corners are zero, too.
 *
 * Author:	Paul Wessel
 * Date:	19-JUL-1989
 * Revised:	6-JAN-1990	PW: Updated to v.2.0
 * Revised:	16-JUN-1998	PW: Updated to v.3.1
 * Version:	4
 */

#include "gmt.h"

struct GRDSAMPLE_CTRL {
	struct F {	/* -F */
		BOOLEAN active;
	} F;
	struct G {	/* -G<grdfile> */
		BOOLEAN active;
		char *file;
	} G;
	struct I {	/* -Idx[/dy] */
		BOOLEAN active;
		double xinc, yinc;
	} I;
	struct L {	/* -L<flag> */
		BOOLEAN active;
		char mode[4];
	} L;
	struct Q {	/* -Q[<value>] */
		BOOLEAN active;
		double value;
	} Q;
	struct T {	/* -T */
		BOOLEAN active;
	} T;
};

int main (int argc, char **argv)
{
	BOOLEAN error = FALSE, greenwich = FALSE, area_set = FALSE, n_set = FALSE, node = FALSE;

	char *infile = CNULL, format[BUFSIZ];

	int i, j, ii, jj, mx = 0;

	size_t ij, k, nm;
	
	float *a, *b;

	double *lon, lat;

	struct GRD_HEADER grd_a, grd_b;
	struct GMT_EDGEINFO edgeinfo;
	struct GMT_BCR bcr;
	struct GRDSAMPLE_CTRL *Ctrl;

	void *New_Grdsample_Ctrl (), Free_Grdsample_Ctrl (struct GRDSAMPLE_CTRL *C);

	argc = GMT_begin (argc, argv);

	Ctrl = (struct GRDSAMPLE_CTRL *)New_Grdsample_Ctrl ();	/* Allocate and initialize a new control structure */
	
	GMT_grd_init (&grd_b, argc, argv, FALSE);

	GMT_boundcond_init (&edgeinfo);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {

				/* Common parameters */

				case 'R':
				case 'V':
				case 'f':
				case '\0':
					error += GMT_parse_common_options (argv[i], &grd_b.x_min, &grd_b.x_max, &grd_b.y_min, &grd_b.y_max);
					break;

				/* Supplemental parameters */

				case 'F':
					Ctrl->F.active = TRUE;
					break;
				case 'G':
					Ctrl->G.active = TRUE;
					Ctrl->G.file = strdup (&argv[i][2]);
					break;
				case 'N':	/* Backwards compatible.  nx/ny can now be set with -I */
					sscanf (&argv[i][2], "%d/%d", &grd_b.nx, &grd_b.ny);
					if (grd_b.ny == 0) grd_b.ny = grd_b.nx;
					n_set = TRUE;
					break;
				case 'I':
					Ctrl->I.active = TRUE;
					if (GMT_getinc (&argv[i][2], &Ctrl->I.xinc, &Ctrl->I.yinc)) {
						GMT_inc_syntax ('I', 1);
						error = TRUE;
					}
					break;
				case 'L':
					Ctrl->L.active = TRUE;
					strncpy (Ctrl->L.mode, &argv[i][2], 4);
					break;
				case 'Q':
					Ctrl->Q.active = TRUE;
					Ctrl->Q.value = (argv[i][2]) ? atof (&argv[i][2]) : 1.0;
					break;
				case 'T':	/* Convert from pixel file <-> gridfile */
					Ctrl->T.active = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else 
			infile = argv[i];
	}

	if (argc == 1 || GMT_give_synopsis_and_exit) {
		fprintf (stderr, "grdsample %s - Resample a gridded file onto a new grid\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdsample <old_grdfile> -G<new_grdfile> [-F] [%s] [-L<flag>]\n", GMT_I_OPT);
		fprintf (stderr, "\t[-Q[<value>]] [%s] [-T] [-V] [%s]\n", GMT_Rgeo_OPT, GMT_f_OPT);

		if (GMT_give_synopsis_and_exit) exit (EXIT_FAILURE);

		fprintf (stderr, "\t<old_grdfile> is data set to be resampled\n");
		fprintf (stderr, "\t-G sets the name of the interpolated output grdfile\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-F force pixel node registration  [Default is centered]\n");
		GMT_inc_syntax ('I', 0);
		fprintf (stderr, "\t-L sets boundary conditions.  <flag> can be either\n");
		fprintf (stderr, "\t   g for geographic boundary conditions\n");
		fprintf (stderr, "\t   or one or both of\n");
		fprintf (stderr, "\t   x for periodic boundary conditions on x\n");
		fprintf (stderr, "\t   y for periodic boundary conditions on y\n");
		fprintf (stderr, "\t-Q Quick mode, use bilinear rather than bicubic interpolation.\n");
		fprintf (stderr, "\t   Optionally, append <value> in the 0 < value <= 1 range.\n");
		fprintf (stderr, "\t   [Default = 1 requires all 4 nodes to be non-NaN.], <value> = 0.5\n");
		fprintf (stderr, "\t   will interpolate about 1/2 way from a non-NaN to a NaN node, while\n");
		fprintf (stderr, "\t   0.1 will go about 90%% of the way, etc. -Q0 will return the value of the \n");
		fprintf (stderr, "\t   nearest node instead of interpolating.\n");
		fprintf (stderr, "\t-R specifies a subregion [Default is old region]\n");
		fprintf (stderr, "\t-T Toggles between grid registration and pixel registration\n");
		fprintf (stderr, "\t   In this case, -I is ignored.\n");
		GMT_explain_option ('V');
		GMT_explain_option ('f');
		fprintf (stderr, "\t   One only of -N -I must be specified\n");

		exit (EXIT_FAILURE);
	}

	if (Ctrl->Q.active && (Ctrl->Q.value < 0.0 || Ctrl->Q.value > 1.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Q:  threshold must be in [0,1] range\n", GMT_program);
		error++;
	}
	if (!infile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (!Ctrl->G.file) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", GMT_program);
		error++;
	}
	if (Ctrl->L.active && GMT_boundcond_parse (&edgeinfo, Ctrl->L.mode)) error++;

	if (grd_b.x_min != grd_b.x_max && grd_b.y_min != grd_b.y_max) area_set = TRUE;
	if (Ctrl->I.active) {
		grd_b.x_inc = Ctrl->I.xinc;
		grd_b.y_inc = Ctrl->I.yinc;
	}
	
	if (GMT_read_grd_info (infile, &grd_a)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, infile);
		exit (EXIT_FAILURE);
	}

	if (!area_set) {
		grd_b.x_min = grd_a.x_min;
		grd_b.x_max = grd_a.x_max;
		grd_b.y_min = grd_a.y_min;
		grd_b.y_max = grd_a.y_max;
	}
	if (!Ctrl->T.active) {
		if (Ctrl->I.active && n_set) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Only one of -I, -N may be specified\n", GMT_program);
			error++;
		}
		if (n_set && (grd_b.nx <= 0 || grd_b.ny <= 0)) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -N:  Must specify positive integers\n", GMT_program);
			error++;
		}
		if (Ctrl->I.active) GMT_RI_prepare (&grd_b);	/* Ensure -R -I consistency and set nx, ny */
		if (Ctrl->I.active && (grd_b.x_inc <= 0.0 || grd_b.y_inc <= 0.0)) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -I:  Must specify positive increments\n", GMT_program);
			error++;
		}
		if (!(Ctrl->I.active || n_set)) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  One of -I, -N must be specified\n", GMT_program);
			error++;
		}
	}

	if (error) exit (EXIT_FAILURE);

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

	GMT_boundcond_param_prep (&grd_a, &edgeinfo);

	if (Ctrl->T.active) {
		Ctrl->F.active = !grd_a.node_offset;	/* Change to the opposite of what it is */
		grd_b.nx = (Ctrl->F.active) ? grd_a.nx - 1 : grd_a.nx + 1;
		grd_b.ny = (Ctrl->F.active) ? grd_a.ny - 1 : grd_a.ny + 1;
		area_set = Ctrl->I.active = FALSE;
	}

	a = (float *) GMT_memory (VNULL, (size_t)((grd_a.nx + 4) * (grd_a.ny + 4)), sizeof(float), GMT_program);

	if (area_set) {
		if (grd_b.y_min < grd_a.y_min || grd_b.y_max > grd_a.y_max) {
			fprintf (stderr, "%s:  Selected region exceeds the boundaries of the grdfile!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		else if (!edgeinfo.nxp && (grd_b.x_min < grd_a.x_min || grd_b.x_max > grd_a.x_max)) {
			fprintf (stderr, "%s:  Selected region exceeds the boundaries of the grdfile!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
	}

	if (!Ctrl->F.active && !Ctrl->T.active) Ctrl->F.active = grd_a.node_offset;
	grd_b.node_offset = Ctrl->F.active;
	grd_b.xy_off = 0.5 * grd_b.node_offset;

	if (edgeinfo.nxp && grd_b.x_min < 0.0 && grd_b.x_max > 0.0)
		greenwich = TRUE;
	else if (edgeinfo.nxp && grd_b.x_max > 360.0) {
		greenwich = TRUE;
		grd_b.x_min -= 360.0;
		grd_b.x_max -= 360.0;
	}
	if (!Ctrl->I.active) {
		grd_b.x_inc = GMT_get_inc (grd_b.x_min, grd_b.x_max, grd_b.nx, grd_b.node_offset);
		grd_b.y_inc = GMT_get_inc (grd_b.y_min, grd_b.y_max, grd_b.ny, grd_b.node_offset);
	}

	GMT_grd_RI_verify (&grd_b, 1);

	nm = ((size_t)grd_b.nx) * ((size_t)grd_b.ny);
	b = (float *) GMT_memory (VNULL, nm, sizeof(float), GMT_program);

	sprintf (format, "%%s: New grid (%s/%s/%s/%s) nx = %%d ny = %%d dx = %s dy = %s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
	if (gmtdefs.verbose) fprintf (stderr, format, GMT_program, grd_b.x_min, grd_b.x_max, grd_b.y_min, grd_b.y_max, grd_b.nx, grd_b.ny, grd_b.x_inc, grd_b.y_inc);

	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2;	/* Leave room for 2 empty boundary rows/cols */

	if (GMT_read_grd (infile, &grd_a, a, grd_a.x_min, grd_a.x_max, grd_a.y_min, grd_a.y_max, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, infile);
		exit (EXIT_FAILURE);
	}

	/* Initialize bcr structure:  */

	GMT_bcr_init (&grd_a, GMT_pad, Ctrl->Q.active, Ctrl->Q.value, &bcr);

	/* Set boundary conditions  */

	GMT_boundcond_set (&grd_a, &edgeinfo, GMT_pad, a);

	/* Precalculate longitudes */

	lon = (double *) GMT_memory (VNULL, (size_t)grd_b.nx, sizeof (double), GMT_program);
	for (i = 0; i < grd_b.nx; i++) {
		lon[i] = GMT_i_to_x (i, grd_b.x_min, grd_b.x_max, grd_b.x_inc, grd_b.xy_off, grd_b.nx);
		if (edgeinfo.nxp && greenwich && lon[i] > 180.0) lon[i] -= 360.0;
	}

	if (Ctrl->Q.value == 0.0) {
		node = TRUE;
		mx = grd_a.nx + 4;
	}
	
	for (j = ij = 0; j < grd_b.ny; j++) {
		lat = GMT_j_to_y (j, grd_b.y_min, grd_b.y_max, grd_b.y_inc, grd_b.xy_off, grd_b.ny);
		for (i = 0; i < grd_b.nx; i++, ij++) {
			if (node) {
				ii = GMT_x_to_i (lon[i], grd_a.x_min, grd_a.x_inc, grd_a.xy_off, grd_a.nx);
				jj = GMT_y_to_j (lat, grd_a.y_min, grd_a.y_inc, grd_a.xy_off, grd_a.ny);
				k = ((size_t)(j+GMT_pad[3])) * ((size_t)mx) + ((size_t)i) + ((size_t)GMT_pad[0]);

				b[ij] = a[k];
			}
			else
				b[ij] = (float)GMT_get_bcr_z (&grd_a, lon[i], lat, a, &edgeinfo, &bcr);
		}
	}

	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0;	/* No boundary rows/cols on output */
	if (GMT_write_grd (Ctrl->G.file, &grd_b, b, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error writing file %s\n", GMT_program, Ctrl->G.file);
		exit (EXIT_FAILURE);
	}

	GMT_free ((void *)a);
	GMT_free ((void *)b);
	GMT_free ((void *)lon);

	Free_Grdsample_Ctrl (Ctrl);	/* Deallocate control structure */

	GMT_end (argc, argv);

	exit (EXIT_SUCCESS);
}

void *New_Grdsample_Ctrl () {	/* Allocate and initialize a new control structure */
	struct GRDSAMPLE_CTRL *C;
	
	C = (struct GRDSAMPLE_CTRL *) GMT_memory (VNULL, 1, sizeof (struct GRDSAMPLE_CTRL), "New_Grdsample_Ctrl");
	
	/* Initialize values whose defaults are not 0/FALSE/NULL */
	C->Q.value = 1.0;
	return ((void *)C);
}

void Free_Grdsample_Ctrl (struct GRDSAMPLE_CTRL *C) {	/* Deallocate control structure */
	if (C->G.file) GMT_free ((void *)C->G.file);	
	GMT_free ((void *)C);	
}
