/*--------------------------------------------------------------------
 *    $Id: blockmedian.c,v 1.5.4.3 2002/02/27 17:58:55 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
 *--------------------------------------------------------------------*/

/*
   blockmedian.c
   Takes lon, lat, data, [weight] on GMT_stdin or file and writes out one value
   per cell, where cellular region is bounded by West East South North and
   cell dimensions are delta_x, delta_y.
      
   Author: 	Walter H. F. Smith
   Date:	28 June, 1988
   Modified	26 April, 1991 for gmt v2.0 by whfs smith;
		added dynamic memory allocation.
   Modified:	3 Jan 1995 by PW for gmt 3.0
   Modified:	3 May 1998 by PW for gmt 3.1
   Modified:	29 Oct 1998 by WHFS to add -Q option
   Modified:	3.3.5: 10 Jul 2000 by PW to add -L
   Version:	3.4 01-MAR-2001 by PW, Use -F instead of -N, and add -C
   Version:	3.4.1
*/

#include "gmt.h"

struct DATA {
        int     i;
        float   w;
        float   x;
        float   y;
        float   z;
}       *data;


int	compare_x(const void *point_1, const void *point_2);
int	compare_y(const void *point_1, const void *point_2);
int	compare_index_z(const void *point_1, const void *point_2);
void median_output (int first_in_cell, int first_in_new_cell, double west, double south, double delta_x, double delta_y, double xy_off, int nx, double weight_sum, double *xx, double *yy, double *zz, int go_quickly);

main (int argc, char **argv)
{	
	
	BOOLEAN	error, weighted, offset, report, nofile = TRUE, done = FALSE, first = TRUE, median_xy = TRUE;
	/* Default value for go_quickly = FALSE for backward compatibility with 3.0  */
	
	FILE *fp = NULL;
	
	double	west, east, south, north, delta_x, delta_y, del_off, xy_off;
	double	*in, out[4], idx, idy;
	
	int	i, n_x, n_y, ix, iy, index, first_in_cell, first_in_new_cell, fno, n_files = 0, n_args, n_req;
	int	n_lost, n_read, n_pitched, n_cells_filled, n_alloc, n_expected_fields, n_fields, n_out, go_quickly = 0;
	
	char	modifier, buffer[BUFSIZ], format[BUFSIZ];

	argc = GMT_begin (argc, argv);
	
	west = east = south = north = delta_x = delta_y = 0.0;
	del_off = 0.5;
	error = weighted = offset = report = FALSE;

	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], &west, &east, &south, &north);
                                      break;
                              
				/* Supplemental parameters */
                              
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
					
				case 'C':
					median_xy = FALSE;
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &delta_x, &delta_y);
					break;
				case 'L':
					GMT_geographic_in = TRUE;
					break;
				case 'N':	/* Backward compatible with 3.3.6 */
				case 'F':
					offset = TRUE;
					break;
				case 'Q':
					go_quickly = 1;
					break;
				case 'W':
					if ( (modifier = argv[i][2]) == 'i' || modifier == 'I') {
						weighted = TRUE;
						report = FALSE;
					}
					else if (modifier == 'O' || modifier == 'o') {
						report = TRUE;
						weighted = FALSE;
					}
					else
						weighted = report = TRUE;
					break;
					
				default:
					error = TRUE;
                                        GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "blockmedian %s - Block averaging by L1 norm\n\n", GMT_VERSION);
		fprintf (stderr, "usage: blockmedian [infile(s)] -I<xinc[m]>[/<yinc>[m]] -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-C] [-F] [-H[<nrec>]] [-L] [-Q] [-V] [-W[i][o] ] [-:] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		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 s] to xinc or yinc to indicate minutes [or seconds],\n");
		fprintf (stderr, "\t   e.g., -I10m/5m grids longitude every 10 minutes, latitude every 5 minutes.\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C Output center of block and median z-value  [Default is median location (but see -Q)]\n");
		fprintf (stderr, "\t-F Offsets registration so block edges are on gridlines (pixel reg.).  [Default: grid reg.]\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-L means that x is longitude, i.e. assumed to be periodic in 360\n");
		fprintf (stderr, "\t-Q Quicker; get median z and x,y at that z.  [Default gets median x, median y, median z.]\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W sets Weight options.  -WI reads Weighted Input (4 cols: x,y,z,w) but writes only (x,y,z) Output.\n");
		fprintf (stderr, "\t   -WO reads unWeighted Input (3 cols: x,y,z) but reports sum (x,y,z,w) Output.\n");
		fprintf (stderr, "\t   -W with no modifier has both weighted Input and Output; Default is no weights used.\n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 3 (or 4 if -W is set).\n");
		GMT_explain_option ('o');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	if (!median_xy && go_quickly == 1) {
		fprintf (stderr, "%s: GMT WARNING:  -C overrides -Q\n", GMT_program);
		go_quickly = 0;
	}

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

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	idx = 1.0 / delta_x;
	idy = 1.0 / delta_y;
	n_x = irint ((east - west) * idx) + 1;
	n_y = irint ((north - south) * idy) + 1;
	
	if (!median_xy) go_quickly = 2;	/* Flag used in output calculation */
	if (offset) {
		n_x--;
		n_y--;
		del_off = 0.0;
	}
	xy_off = 0.5 - del_off;
	
	if (gmtdefs.verbose) {
		sprintf (format, "%%s: W: %s E: %s S: %s N: %s nx: %%d ny: %%d\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, west, east, south, north, n_x, n_y);
	}

	n_read = n_pitched = 0;
        n_alloc = GMT_CHUNK;
        data = (struct DATA *) GMT_memory (VNULL, (size_t)n_alloc, sizeof(struct DATA), GMT_program);
	
	/* Read the input data  */

	n_expected_fields = (GMT_io.binary[0]) ? GMT_io.ncol[0] : 3 + weighted;

	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
	
	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++) {
				GMT_fgets (buffer, BUFSIZ, fp);
				buffer[strlen(buffer)-1] = 0;
				if (first) (report && !(weighted)) ? fprintf (GMT_stdout, "%s weights\n", buffer) : fprintf (GMT_stdout, "%s\n", buffer);
			}
			first = FALSE;
		}
				
		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */
				
			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);
			}

			if (GMT_is_dnan (in[2])) continue;	/* Skip when z = NaN */

			n_read++;
		
			ix = (int)floor(((in[0] - west) * idx) + del_off);
			if ( ix < 0 || ix >= n_x ) continue;
			iy = (int)floor(((in[1] - south) * idy) + del_off);
			if ( iy < 0 || iy >= n_y ) continue;
		
			index = iy * n_x + ix;
		
			data[n_pitched].i = index;
			data[n_pitched].w = (float)((weighted) ? in[3] : 1.0);
			if (median_xy) {
				data[n_pitched].x = (float)in[0];
				data[n_pitched].y = (float)in[1];
			}
			data[n_pitched].z = (float)in[2];
		
			n_pitched++;
			if (n_pitched == n_alloc) {
				n_alloc += GMT_CHUNK;
				data = (struct DATA *) GMT_memory ((void *)data, (size_t)n_alloc, sizeof(struct DATA), GMT_program);
			}
		}
		if (fp != GMT_stdin) GMT_fclose(fp);
		
	}
	
	data = (struct DATA *) GMT_memory ((void *)data, (size_t)n_pitched, sizeof(struct DATA), GMT_program);
	n_lost = n_read - n_pitched;
	if (gmtdefs.verbose) fprintf(stderr,"%s: N read: %d N used: %d N outside_area: %d\n", GMT_program, 
		n_read,n_pitched,n_lost);


	/* Ready to go. */
	
	n_out = (report) ? 4 : 3;
	
	/* Sort on index and Z value */
	
	qsort((void *)data, (size_t)n_pitched, sizeof (struct DATA), compare_index_z);
	
	/* Find n_in_cell and write appropriate output  */

	first_in_cell = 0;
	n_cells_filled = 0;
	while (first_in_cell < n_pitched) {
		out[3] = data[first_in_cell].w;
		first_in_new_cell = first_in_cell + 1;
		while ( (first_in_new_cell < n_pitched) && (data[first_in_new_cell].i == data[first_in_cell].i) ) {
			out[3] += data[first_in_new_cell].w;
			first_in_new_cell++;
		}
		median_output(first_in_cell, first_in_new_cell, west, south, delta_x, delta_y, xy_off, n_x, out[3], &out[0], &out[1], &out[2], go_quickly);
		
		GMT_output (GMT_stdout, n_out, out);
		
		n_cells_filled++;
		first_in_cell = first_in_new_cell;
	}

	if (gmtdefs.verbose) {
		fprintf(stderr,"%s: N_cells_filled: %d\n", GMT_program, n_cells_filled);
	}

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

void median_output (int first_in_cell, int first_in_new_cell, double west, double south, double delta_x, double delta_y, double xy_off, int nx, double weight_sum, double *xx, double *yy, double *zz, int go_quickly)
{
	double	weight_half, weight_count;
	int	index, n_in_cell;
	
	weight_half  = 0.5 * weight_sum;
	n_in_cell = first_in_new_cell - first_in_cell;
	
	/* Data are already sorted on z  */

	index = first_in_cell;
	weight_count = data[first_in_cell].w;
	while (weight_count < weight_half) {
		index++;
		weight_count += data[index].w;
	}
	if ( weight_count == weight_half ) {
		*xx = 0.5 * (data[index].x + data[index + 1].x);
		*yy = 0.5 * (data[index].y + data[index + 1].y);
		*zz = 0.5 * (data[index].z + data[index + 1].z);
	}
	else {
		*xx = data[index].x;
		*yy = data[index].y;
		*zz = data[index].z;
	}
	
	/* Now get median x and median y if quick x and quick y not wanted:  */

	if (go_quickly == 1) return;
	
	if (go_quickly == 2) {	/* Get center of block */
		int i, j;
		j = data[index].i / nx;
		i = data[index].i % nx;
		*xx = west  + (i + xy_off) * delta_x;
		*yy = south + (j + xy_off) * delta_y;
		return;
	}
	
	if (n_in_cell > 2)
		qsort((void *)&data[first_in_cell], (size_t)n_in_cell, sizeof (struct DATA), compare_x);
	index = first_in_cell;
	weight_count = data[first_in_cell].w;
	while (weight_count < weight_half) {
		index++;
		weight_count += data[index].w;
	}
	if ( weight_count == weight_half )
		*xx = 0.5 * (data[index].x + data[index + 1].x);
	else
		*xx = data[index].x;

	if (n_in_cell > 2)
		qsort((void *)&data[first_in_cell], (size_t)n_in_cell, sizeof (struct DATA), compare_y);
	index = first_in_cell;
	weight_count = data[first_in_cell].w;
	while (weight_count < weight_half) {
		index++;
		weight_count += data[index].w;
	}
	if ( weight_count == weight_half )
		*yy = 0.5 * (data[index].y + data[index + 1].y);
	else
		*yy = data[index].y;
}

int	compare_index_z (const void *point_1, const void *point_2)
{
	int	index_1, index_2;
	float	data_1, data_2;
	struct DATA *p1, *p2;
	
	p1 = (struct DATA *)point_1;
	p2 = (struct DATA *)point_2;
	index_1 = p1->i;
	index_2 = p2->i;

	if (index_1 < index_2)
		return (-1);
	else if (index_1 > index_2)
		return (1);
	else {
		data_1 = p1->z;
		data_2 = p2->z;
		if (data_1 < data_2)
			return (-1);
		else if (data_1 > data_2)
			return (1);
		else
			return (0);
	}
}
int	compare_x(const void *point_1, const void *point_2)
{
	int	index_1, index_2;
	float	x_1, x_2;
	struct DATA *p1, *p2;
	
	p1 = (struct DATA *)point_1;
	p2 = (struct DATA *)point_2;
	index_1 = p1->i;
	index_2 = p2->i;

	if (index_1 < index_2)
		return (-1);
	else if (index_1 > index_2)
		return (1);
	else {
		x_1 = p1->x;
		x_2 = p2->x;
		if (x_1 < x_2)
			return (-1);
		else if (x_1 > x_2)
			return (1);
		else
			return (0);
	}
}

int	compare_y(const void *point_1, const void *point_2)
{
	int	index_1, index_2;
	float	y_1, y_2;
	struct DATA *p1, *p2;
	
	p1 = (struct DATA *)point_1;
	p2 = (struct DATA *)point_2;
	index_1 = p1->i;
	index_2 = p2->i;

	if (index_1 < index_2)
		return (-1);
	else if (index_1 > index_2)
		return (1);
	else {
		y_1 = p1->y;
		y_2 = p2->y;
		if (y_1 < y_2)
			return (-1);
		else if (y_1 > y_2)
			return (1);
		else
			return (0);
	}
}