/* $Id: img2mercgrd.c,v 1.2.4.1 2002/02/27 17:41:10 pwessel Exp $
 *
 * img2mercgrd.c
 *
 * img2mercgrd reads an "img" format file (used by Sandwell and Smith
 * in their estimates of gravity and depth from satellite altimetry),
 * extracts a Region [with optional N by N averaging], and writes out
 * the data [or Track control] as a pixel-registered GMT "grd" file, 
 * preserving the Mercator projection (gmtdefaults ELLIPSOID = Sphere) 
 * inherent in the img file.
 *
 * img file is read from dir GMT_IMGDIR if this is defined as an
 * environment variable; else img file name is opened directly.
 *
 * The coordinates in the img file are determined by the latitude and
 * longitude span of the img file and the width of an img pixel in
 * minutes of longitude.  These may be changed from their default values
 * using the lower-case [-x<maxlon>] [-y<minlat>/<maxlat>] [-m<minutes>]
 * but must be set to match the coverage of the input img file.
 *
 * The user-supplied -R<w>/<e>/<s>/<n> will be adjusted by rounding
 * outward so that the actual range spanned by the file falls exactly on
 * the edges of the input pixels.  If the user uses the option to 
 * average the input pixels into N by N square block averages, then the
 * -R will be adjusted so that the range spans the block averages.
 *
 * The coordinates in the output file are in spherical mercator projected
 * units ranging from 0,0 in the lower left corner of the output to
 * xmax, ymax in the upper right corner, where xmax,ymax are the width
 * and height of the (spherical) Mercator map with -Rww/ee/ss/nn and -Jm1.
 * Here ww/ee/ss/nn are the outward-rounded range values.
 *
 * Further details on the coordinate system can be obtained from the
 * comments in gmt_imgsubs.c and gmt_imgsubs.h
 *
 * This is a complete rebuild (for v3.1) of the old program by this name.
 * New functionality added here is the averaging option [-N] and the
 * options to define -m, -x, -y, which permit handling very early versions
 * of these files and anticipate higher-resolutions to come down the road.
 * Also added is the option to look for the img file in GMT_IMGDIR
 *
 * Author:  Walter H F Smith
 * Date:    8 October, 1998
 * 
 **WHFS 27 Nov 1998 fixed bug so that -T0 gives correct results
 **  PW 18 Oct 1999 Use WORDS_BIGENDIAN macro set by configure.
 *
 */


#include "gmt_imgsubs.h"

main (int argc, char **argv)
{
	int	i, ij, navgsq, error, irarg = -1, pad[4], tempint;
	int	navg = 1;	/* navg by navg pixels are averaged if navg > 1; else if navg == 1 do nothing */
	int	output_type = -1;	/* force user to choose something between 0 and 3  */
	int	iout, jout, iinstart, iinstop, jinstart, jinstop, k, kk, ion, jin, jj, iin, ii, kstart;
	double	scale_factor = 1.0;	/* use 0.1 for grav,curv,north,east files.  */
	double	west, east, south, north, toplat, botlat, dx, rnavgsq, csum, dsum;
	struct GRD_HEADER h;
	struct GMT_IMG_COORD imgcoord;
	float	*a, empty_val;
	short int *row;
	int	*ix;
	char	*infilename, *outfile, *dirstem, infile[320];
	FILE	*fp;

	struct GMT_IMG_RANGE imgrange = { GMT_IMG_MAXLON, GMT_IMG_MINLAT, GMT_IMG_MAXLAT, GMT_IMG_MPIXEL };

	argc = GMT_begin (argc, argv);
	
	GMT_grd_init (&h, argc, argv, FALSE);
	
	empty_val = GMT_f_NaN;

	error = FALSE;
	west = east = south = north = 0.0;
	fp = NULL;
	outfile = CNULL;
	infilename = CNULL;

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

				/* Common parameters:  */

				case 'R':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					if (argv[i][1] == 'R') irarg = i;
					break;

				/* Supplemental parameters:  */

				case 'G':
					outfile = &argv[i][2];
					break;
				case 'N':
					if ((sscanf(&argv[i][2], "%d", &navg)) != 1 || navg < 1) {
						error++;
						fprintf(stderr,"%s:  GMT SYNTAX ERROR -N requires an integer > 1.\n", GMT_program);
					}
					break;
				case 'S':
					if ((sscanf(&argv[i][2], "%lf", &scale_factor)) != 1) {
						error++;
						fprintf(stderr,"%s:  GMT SYNTAX ERROR -S requires a scale factor.\n", GMT_program);
					}
					break;				
				case 'T':
					if ((sscanf(&argv[i][2], "%d", &output_type)) != 1 || output_type < 0 || output_type > 3) {
						error++;
						fprintf(stderr,"%s:  GMT SYNTAX ERROR -T requires an output type 0,1,2 or 3.\n", GMT_program);
					}
					break;
				case 'm':
					if ((sscanf(&argv[i][2], "%lf", &imgrange.mpixel)) != 1 || imgrange.mpixel <= 0.0) {
						error++;
						fprintf(stderr,"%s:  GMT SYNTAX ERROR -m requires a positive value.\n", GMT_program);
					}
					break;
				case 'x':
					if ((sscanf(&argv[i][2], "%lf", &imgrange.maxlon)) != 1 || imgrange.maxlon < 360.0) {
						error++;
						fprintf(stderr,"%s:  GMT SYNTAX ERROR -x requires a number >= 360.0.\n", GMT_program);
					}
					break;
				case 'y':
					if ((sscanf(&argv[i][2], "%lf/%lf", &imgrange.minlat, &imgrange.maxlat)) != 2) {
						error++;
						fprintf(stderr,"%s:  GMT SYNTAX ERROR -y requires <minlat>/<maxlat>.\n", GMT_program);
					}
					break;
				default:
					error++;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else {
			infilename = argv[i];
		}
	}

	if (GMT_quick || argc == 1) {
		fprintf(stderr,"img2mercgrd %s - Extract a grd file from an img file, preserving Mercator projection.\n\n", GMT_VERSION);
		fprintf(stderr,"usage:  img2mercgrd <world_image_filename> -R<w>/<e>/<s>/<n> -G<grdfile> -T<type> \n");
		fprintf(stderr,"\t\t[-N<navg>] [-S<scale>] [-V] [-m<minutes>] [-x<maxlon>] [-y<minlat>/<maxlat>]\n\n");

		if (GMT_quick) exit (-1);
		
		fprintf (stderr, "\tREQUIRED ARGUMENTS:\n");
		fprintf(stderr,"\t<world_image_filename> gives location of img file.\n");
		fprintf(stderr,"\t-G sets filename for the output .grd format file.\n");
		fprintf(stderr,"\t-R specifies the region in degrees or degrees:minutes.\n");
		fprintf(stderr,"\t-T selects output type:\n");
		fprintf(stderr,"\t\t-T0 for obsolete img files w/ no constraint code, gets data.\n");
		fprintf(stderr,"\t\t-T1 for new img file w/ constraints coded, gets data at all points.\n");
		fprintf(stderr,"\t\t-T2 for new img file w/ constraints coded, gets data only at constrained points, NaN elsewhere.\n");
		fprintf(stderr,"\t\t-T3 for new img file w/ constraints coded, gets 1 at constraints, 0 elsewhere.\n\n");
		
		fprintf (stderr, "\tOPTIONAL ARGUMENTS:\n");
		fprintf(stderr,"\t-N<navg> will ouput averages of input in navg by navg squares.  [no averaging]\n");
		fprintf(stderr,"\t-S<scale> will multiply img integer values by scale for output [%lg]\n", scale_factor);
		GMT_explain_option ('V');
		
		fprintf (stderr, "\n\tOPTIONAL ARGUMENTS WHICH DEPEND ON THE img FILE VERSION:\n");
		fprintf (stderr, "\t-m<minutes> input img pixels are <minutes> minutes of longitude wide. [2.0]\n");
		fprintf (stderr, "\t-x<maxlon> input img file runs from 0 to <maxlon> longitude. [360.0]\n");
		fprintf (stderr, "\t-y<minlon>/<maxlon> input img file bottom and top latitudes. [-72.006/72.006]\n");
		exit(-1);
	}

	if (irarg < 0 || west >= east || south >= north) {
		fprintf(stderr, "%s: GMT SYNTAX ERROR:  Must specify -R with west < east and south < north.\n", GMT_program);
		error++;
	}
	if (outfile == NULL) {
		fprintf(stderr, "%s: GMT SYNTAX ERROR:  Must specify output grdfile name with -G.\n", GMT_program);
		error++;
	}
	if (infilename == NULL) {
		fprintf(stderr, "%s: GMT SYNTAX ERROR:  Must specify input imgfile name.\n", GMT_program);
		error++;
	}
	if (output_type < 0 || output_type > 3) {
		fprintf(stderr, "%s: GMT SYNTAX ERROR:  Must specify output type in the range 0-3 with -T.\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	if (GMT_img_setup_coord ( &imgrange, &imgcoord) ) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Error in img coordinate specification [-m -x or -y]\n", GMT_program);
		error++;
	}
	else if (navg && (imgcoord.nx360%navg != 0 || imgcoord.nyrow%navg != 0) ) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Bad choice of navg in -N.  Must divide %d and %d\n", GMT_program,
			imgcoord.nx360, imgcoord.nyrow);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

	if ( (dirstem = getenv ("GMT_IMGDIR")) == CNULL) {
		sprintf (infile, "%s", infilename);
	}
	else {
		sprintf (infile, "%s%c%s", dirstem, DIR_DELIM, infilename);
	}
	if ((fp = fopen(infile, "rb")) == NULL) {
		fprintf(stderr, "%s: GMT SYNTAX ERROR:  Cannot open %s for binary read.\n", GMT_program, infile);
		error++;
	}
	if (error) exit (EXIT_FAILURE);
	
	
	/* Expected edges of input image based on coordinate initialization (might not exactly match user spec):  */
	toplat = GMT_img_ypix_to_lat (0.0, &imgcoord);
	botlat = GMT_img_ypix_to_lat ( (double)imgcoord.nyrow, &imgcoord);
	dx = 1.0 / ( (double)imgcoord.nx360 / 360.0);
	if (gmtdefs.verbose)
		fprintf (stderr, "%s expects %s to be %d by %d pixels spanning 0/%5.1lf/%.8lg/%.8lg\n",
			GMT_program, infile, imgcoord.nxcol, imgcoord.nyrow, dx*imgcoord.nxcol, botlat, toplat);

	if (toplat < north)
		fprintf(stderr, "%s:  WARNING:  Your top latitude (%.12lg) lies outside top latitude of input (%.12lg).\n",
			GMT_program, north, toplat);
	if (botlat > south)
		fprintf(stderr, "%s:  WARNING:  Your bottom latitude (%.12lg) lies outside bottom latitude of input (%.12lg).\n",
			GMT_program, south, botlat);
	
	/* Re-adjust user's -R so that it falls on pixel coordinate boundaries:  */
	
	jinstart = navg * (int)floor (GMT_img_lat_to_ypix (north, &imgcoord) / navg);
	jinstop  = navg * (int)ceil  (GMT_img_lat_to_ypix (south, &imgcoord) / navg);
	/* jinstart <= jinputrow < jinstop  */
	h.ny = (jinstop - jinstart) / navg;
	north = GMT_img_ypix_to_lat ((double)jinstart, &imgcoord);
	south = GMT_img_ypix_to_lat ((double)jinstop,  &imgcoord);

	iinstart = navg * (int)floor (west/(dx*navg));
	iinstop  = navg * (int)ceil  (east/(dx*navg));
	/* iinstart <= ipixelcol < iinstop, but modulo all with imgcoord.nx360  */
	/* Reset left and right edges of user area:  */
	west = iinstart * dx;
	east = iinstop  * dx;
	h.nx = (iinstop - iinstart) / navg;

	if (gmtdefs.verbose) {
		fprintf(stderr, "%s:  To fit [averaged] input, your %s is adjusted to -R%.12lg/%.12lg/%.12lg/%.12lg\n",
			GMT_program, argv[irarg], west, east, south, north);
		fprintf (stderr, "%s:  The output will be %d by %d pixels.\n", GMT_program, h.nx, h.ny);
	}

	/* Set iinstart so that it is non-negative, for use to index pixels.  */
	while (iinstart < 0) iinstart += imgcoord.nx360;
	
	/* Set navgsq, rnavgsq, for the averaging:  */
	navgsq = navg * navg;
	rnavgsq = 1.0 / navgsq;

	/* Set up header with Mercatorized dimensions assuming -Jm1  */
	h.x_min = 0.0;
	h.x_max = h.nx * navg * dx;
	h.y_min = 0.0;
	h.y_max = h.ny * navg * dx;
	h.x_inc = navg * dx;
	h.y_inc = navg * dx;
	h.node_offset = 1;
	h.z_scale_factor = 1.0;
	h.z_add_offset = 0.0;
	strcpy (h.x_units, "Spherical Mercator projected Longitude, -Jm1, length from West Edge.");
	strcpy (h.y_units, "Spherical Mercator projected Latitude, -Jm1, length from South Edge.");
	if (output_type < 3)
		strcpy (h.z_units, "meters, mGal, Eotvos, or micro-radians, depending on img file and -S.");
	else
		strcpy (h.z_units, "T/F, one or more constraints fell in this pixel.");
	strcpy (h.title, "Data from Altimetry");
	sprintf (h.remark, "Spherical Mercator Projected with -Jm1 -R%.12lg/%.12lg/%.12lg/%.12lg", west, east, south, north);
	h.z_min = DBL_MAX;
	h.z_max = -DBL_MAX;

	GMT_grd_RI_verify (&h, 1);

	/* Now malloc some space for float grd array, integer pixel index, and short integer data buffer.  */

	row = (short int *) GMT_memory (VNULL, (size_t)(navg * imgcoord.nxcol), sizeof (short int), GMT_program);
	ix = (int *) GMT_memory (VNULL, (size_t)(navgsq * h.nx), sizeof (int), GMT_program);
	a = (float *) GMT_memory (VNULL, (size_t)(h.nx * h.ny), sizeof (float), GMT_program);

	/* Load ix with the index to the correct column, for each output desired.  This helps for Greenwich, 
		also faster averaging of the file, etc.  Note for averaging each n by n block is looped in turn. */
	
	if (navg > 1) {
		k = 0;
		for (iout = 0; iout < h.nx; iout++) {
			ion = iout * navg;
			for (jin = 0; jin < navg; jin++) {
				jj = jin * imgcoord.nxcol;
				for (iin = 0; iin < navg; iin++) {
					ii = (iin + iinstart + ion) % imgcoord.nx360;
					ix[k] = ii + jj;
					k++;
				}
			}
		}
	}
	else {
		for (iout = 0; iout < h.nx; iout++) {
			ix[iout] = (iout + iinstart) % imgcoord.nx360;
		}
	}


	/* Now before beginning data loop, fseek if needed.  */
	if (jinstart > 0 && jinstart < imgcoord.nyrow) {
		fseek (fp, 2 * imgcoord.nxcol * jinstart, SEEK_SET);
	}
	
	/* Now loop over output points, reading and handling data as needed */

	for (ij = 0, jout = 0; jout < h.ny; jout++) {
		jin = jinstart + navg * jout;
		if (jin < 0 || jin >= imgcoord.nyrow) {
			for (iout = 0; iout < h.nx; iout++, ij++) {
				a[ij] = empty_val;
			}
		}
		else {
			if ( (fread((void *)row, sizeof (short int), (size_t)(navg * imgcoord.nxcol), fp) ) != (size_t)(navg * imgcoord.nxcol) ) {
				fprintf(stderr,"%s:  ERROR:  Read failure at jin = %d.\n", GMT_program, jin);
				exit (EXIT_FAILURE);
			}

#if defined(_WIN32) || WORDS_BIGENDIAN == 0
			for (iout = 0; iout < navg * imgcoord.nxcol; iout++) row[iout] = GMT_swab2 (row[iout]);
#endif

			for (iout = 0, kstart = 0; iout < h.nx; iout++, ij++, kstart += navgsq) {
				if (navg) {
					csum = 0.0;
					dsum = 0.0;
					for (k = 0, kk = kstart; k < navgsq; k++, kk++) {
						tempint = (int)row[ix[kk]];
						if (output_type) {
							if ( ( (abs(tempint))%2) != 0) {
								csum += 1.0;
								tempint--;
							}
						}
						dsum += (double) tempint;
					}
					csum *= rnavgsq;
					dsum *= rnavgsq;
				}
				else {
					tempint = (int)row[ix[iout]];
					if (output_type) {
						if ( ( (abs(tempint))%2) != 0) {
							csum = 1.0;
							tempint--;
						}
						else {
							csum = 0.0;
						}
					}
					dsum = (double) tempint;
				}
				
				dsum *= scale_factor;
				
				switch (output_type) {
					case 0:
					case 1:
						a[ij] = (float) dsum;
						break;
					case 2:
						a[ij] = (float)((csum >= 0.5) ? dsum : GMT_f_NaN);
						break;
					case 3:
						a[ij] = (float)csum;
						break;
				}
				

				if (output_type != 2 || csum >= 0.5) {
					if (h.z_min > a[ij]) h.z_min = a[ij];
					if (h.z_max < a[ij]) h.z_max = a[ij];
				}
			}
		}
	}
	
	fclose (fp);
	if (gmtdefs.verbose)
		fprintf(stderr,"Created %d by %d Mercatorized .grd file.  Min, Max values are %.8lg  %.8lg\n", h.nx, h.ny, h.z_min, h.z_max);

	if (GMT_write_grd (outfile, &h, a, 0.0, 0.0, 0.0, 0.0, pad, FALSE)) {
		fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile);
		exit (EXIT_FAILURE);
	}

	GMT_free((void *)a);
	GMT_free((void *)ix);
	GMT_free((void *)row);
	
        GMT_end (argc, argv);
}