/* img_rescale_prediction.c unscaled.img -N<negative_scale.grd> -P<positive_scale.grd>  > scaled.img
 *
 * 23 Feb 2009 rewritten for new img grid sizes -- WHFS
 *
 * Given global grids -R0/360/-y/y of scale factors for
 * positive and negative predicted values, with y sufficient
 * to cover latitude range of img pixels, read an unscaled
 * img file (highpass topo prediction from g_to_h.com) and
 * rescale it by the values in the grids, which are obtained
 * from the nettleton.com  Bilinear interpolation of the grid
 * is used, and grids are expected to be grid-registered.
 *
 * Requires link with GMT for the grd handling routines.
* On duck, this compiles with
gcc -Wall -O2 img_rescale_prediction.c img_subs.c -I$GMTHOME/src -I/sw/include -L$GMTHOME/lib -L/sw/lib -lgmt -lnetcdf -lm -o img_rescale_prediction
 *
 * Compiles on raptor with:
cc -I$GMTHOME/include -I$NETCDFHOME/include -Aa -g -L$GMTHOME/lib -L$NETCDFHOME/lib img_rescale_prediction.c -lgmt -lnetcdf -lm -o img_rescale_prediction
 *
 * W H F Smith, 03 November 2000
 */
 
#include "img_predict.h"
#include "gmt.h"

int main (int argc, char **argv) {
	
	double	lat, radius, dy, dy1, dx1, a, b, x;
	short int *hbuf;
	float	*pz, *nz;
	double	*dx;
	int	*igrd;
	int	error = 0, i, ii, jimg, jgrd0, nximg, nyimg, nyimghalf;
	size_t	nsize;
	char	*pfile, *nfile, *hfile;
	struct GRD_HEADER pgrd, ngrd;
	FILE	*fp;
	
	pfile = CNULL;
	nfile = CNULL;
	hfile = CNULL;
	

	argc = GMT_begin (argc, argv);
	GMT_grd_init (&ngrd, argc, argv, FALSE);
	GMT_grd_init (&pgrd, argc, argv, FALSE);
	
	/* Parse arg list:  */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				case 'N':
					nfile = &argv[i][2];
					break;
				case 'P':
					pfile = &argv[i][2];
					break;
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], VNULL, VNULL, VNULL, VNULL);
					break;
				default:
					fprintf (stderr, "img_rescale_prediction:  Unsupported option:  %s\n", argv[i]);
					error++;
					break;
			}
		}
		else {
			hfile = argv[i];
		}
	}
	
	if (argc == 1 || error) {
		fprintf (stderr, "img_rescale_prediction %s - Read img on stdin, write rescaled img on stdout\n\n", GMT_VERSION);
		fprintf (stderr, "usage:  img_rescale_prediction input.img -N<negative_scale.grd> -P<positive_scale.grd> [-V] > output.img\n");
		fprintf (stderr, "\t-N and -P must be grid-registered grd files covering -R0/360/-y,y sufficient to cover img area.\n");
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}
	
	if (check_imgfile_size (hfile, &nximg, &nyimg) ) {
		fprintf (stderr, "img_rescale_prediction:  FATAL ERROR.  Cannot understand size of %s\n", hfile);
		exit (EXIT_FAILURE);
	}
	
	if ( (fp = fopen (hfile, "r")) == NULL) {
		fprintf (stderr, "img_rescale_prediction:  FATAL ERROR.  Cannot open %s\n", hfile);
		exit (EXIT_FAILURE);
	}
	
	radius = nximg/(2.0 * M_PI);
	nyimghalf = nyimg/2;
	lat = R2D * (2.0 * atan(exp( (nyimghalf - 0.5) / radius) ) - M_PI_2);     
	nsize = nximg;

 	if ( (dx = (double *) GMT_memory (VNULL, nsize, sizeof (double), GMT_program) ) == (double *)NULL) {
		fprintf (stderr, "%s:  Error allocating array for dx.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
 	if ( (igrd = (int *) GMT_memory (VNULL, nsize, sizeof (int), GMT_program) ) == (int *)NULL) {
		fprintf (stderr, "%s:  Error allocating array for igrd.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
 	if ( (hbuf = (short int *) GMT_memory (VNULL, nsize, sizeof (short int), GMT_program) ) == (short int *)NULL) {
		fprintf (stderr, "%s:  Error allocating array for one row of the img files.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	
	
	if (nfile == (char *)NULL || nfile[0] == '\0') {
		fprintf (stderr, "img_rescale_prediction:  You must supply -N<negative_scale.grd> grd file name.\n");
		error ++;
	}
	else if (GMT_read_grd_info (nfile, &ngrd)) {
		fprintf (stderr, "%s: Error reading header of file %s\n", GMT_program, nfile);
		error++;
	}
	else if (ngrd.x_min != 0.0 || ngrd.x_max != 360.0 || ngrd.y_min > lat || ngrd.y_max < lat) {
		fprintf (stderr, "%s: Error in -R coverage of file %s\n", GMT_program, nfile);
		error++;
	}
	else if (ngrd.node_offset != 0) {
		fprintf (stderr, "%s: Error.  File %s must be grid-registered.\n", GMT_program, nfile);
		error++;
	}
	
	if (pfile == (char *)NULL || pfile[0] == '\0') {
		fprintf (stderr, "img_rescale_prediction:  You must supply -P<positive_scale.grd> grd file name.\n");
		error ++;
	}
	else if (GMT_read_grd_info (pfile, &pgrd)) {
		fprintf (stderr, "%s: Error reading header of file %s\n", GMT_program, pfile);
		error++;
	}
	else if (pgrd.x_min != 0.0 || pgrd.x_max != 360.0 || pgrd.y_min > lat || pgrd.y_max < lat) {
		fprintf (stderr, "%s: Error in -R coverage of file %s\n", GMT_program, pfile);
		error++;
	}
	else if (pgrd.node_offset != 0) {
		fprintf (stderr, "%s:  Error.  File %s must be grid-registered.\n", GMT_program, pfile);
		error++;
	}
	else if (pgrd.x_inc != ngrd.x_inc || pgrd.y_inc != ngrd.y_inc || pgrd.y_min != ngrd.y_min || pgrd.y_max != ngrd.y_max) {
		fprintf (stderr, "%s:  Error.  Grid files must have matching grid meshes.\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);
	
	nsize = ngrd.nx * ngrd.ny;
	
	if ( (nz = (float *) GMT_memory (VNULL, nsize, sizeof (float), GMT_program) ) == (float *)NULL) {
		fprintf (stderr, "%s:  Error allocating array for negative scale grid.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if ( (pz = (float *) GMT_memory (VNULL, nsize, sizeof (float), GMT_program) ) == (float *)NULL) {
		fprintf (stderr, "%s:  Error allocating array for positive scale grid.\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (GMT_read_grd (nfile, &ngrd, nz, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, nfile);
		exit (EXIT_FAILURE);
	}
	if (GMT_read_grd (pfile, &pgrd, pz, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, pfile);
		exit (EXIT_FAILURE);
	}


	if (gmtdefs.verbose) {
		fprintf (stderr, "%s:  Grids are %d by %d and img files are %d by %d\n", 
			GMT_program, pgrd.nx, pgrd.ny, nximg, nyimg);
	}
        
	a = 360.0 / nximg;		/* Pixel width in degrees  */
	b = a / ngrd.x_inc;	/* Pixel width in steps of grd file  */
	for (i = 0; i < nximg; i++) {
		x = (i + 0.5) * b;	/* Pixel center distance from Greenwich in grid steps  */
		igrd[i] = (int) floor(x);
		dx[i] = x - igrd[i];
	}
        
	ii = ngrd.nx * ngrd.ny;
        
	for (i = 0; i < ii; i++) {
		if (GMT_is_fnan(nz[i]) ) {
			fprintf (stderr, "%s:  Error.  Grid file %s contains NaNs.\n", GMT_program, nfile);
			exit (EXIT_FAILURE);
		}
		if (GMT_is_fnan(pz[i]) ) {
			fprintf (stderr, "%s:  Error.  Grid file %s contains NaNs.\n", GMT_program, pfile);
			exit (EXIT_FAILURE);
		}
	}        

	nsize = nximg;
	for (jimg = 0; jimg < nyimg; jimg++) {
        
		if ( (fread ( (void *)hbuf, sizeof (short int), nsize, fp)) != nsize) {
			fprintf (stderr, "%s:  Error reading img row from %s at jimg = %d\n", GMT_program, hfile, jimg);
			exit (EXIT_FAILURE);
		}
        	
		lat = R2D * (2.0 * atan(exp( (nyimghalf - jimg - 0.5) / radius) ) - M_PI_2);
		jgrd0 = (int) floor ( (ngrd.y_max - lat) / ngrd.y_inc);
        	
		/* Special cases for the top and bottom row, which may fall very slightly outside 72:  */
		if (jgrd0 < 0) {
			jgrd0 = 0;
			dy = 0.0;
		}
		else if (jgrd0 > ngrd.ny-2) {
			jgrd0 =  ngrd.ny-2;
			dy = 1.0;
		}
		else {
			dy = (ngrd.y_max - lat) / ngrd.y_inc - jgrd0;
		}
		jgrd0 *= ngrd.nx;
		dy1 = 1.0 - dy;
        	
		for (i = 0; i < nximg; i++) {
			if (hbuf[i] == 0) continue;
			ii = jgrd0 + igrd[i];
			dx1 = (1.0 - dx[i]);
			if (hbuf[i] > 0) {
				a = pz[ii] * dx1 + pz[ii+1] * dx[i];
				ii += ngrd.nx;
				b = pz[ii] * dx1 + pz[ii+1] * dx[i];
			}
			else {
				a = nz[ii] * dx1 + nz[ii+1] * dx[i];
				ii += ngrd.nx;
				b = nz[ii] * dx1 + nz[ii+1] * dx[i];
			}
			x = (int) rint ( (b * dy + a * dy1) * hbuf[i]);
			if (fabs(x) > 32767.0) {
				fprintf (stderr, "%s:  Short Int Overflow at jimg = %d  i = %d\n", GMT_program, jimg, i);
				exit (EXIT_FAILURE);
			}
			hbuf[i] = (short int)x;
		}
        	
		if ( (fwrite ( (void *)hbuf, sizeof (short int), nsize, stdout)) != nsize) {
			fprintf (stderr, "%s:  Error writing img row to stdout at jimg = %d\n", GMT_program, jimg);
			exit (EXIT_FAILURE);
		}
	}
	
	fclose (fp);
	free ( (void *)pz);
	free ( (void *)nz);
	free ( (void *)hbuf);
	free ( (void *)igrd);
	free ( (void *)dx);
	GMT_end (argc, argv);
	exit (EXIT_SUCCESS);
}