/* imgbm_combine.c < binary input > binary output
 *
 * Read and write 8-byte binary records consisting of two 4-byte ints;
 * the first is the index to the img pixel, and the second is the data value.
 * This program re-sorts and re-medians.  It is designed to operate on the
 * output of a job consisting of multiple calls to gmtplus_imgtcbm
 *
 * W H F Smith,  16 July 1996.
 *
 *
 */

#include "../../include/gmt.h"

#define NX 10800
#define NY 6336

struct IMGPAIR {
	int	index;
	int	depth;
};

main(argc, argv)
int	argc;
char	argv;
{
	
	int	i, j, n_alloc, n_in, n_out, n_this_box, n_half, top, compare_imgpairs();
	struct IMGPAIR *imgpair;

	if (argc != 1) {
		fprintf(stderr,"usage:  imgbm_combine < <binary input> > <binary output>\n\n");
		exit(-1);
	}

	n_alloc = GMT_CHUNK;
	imgpair = (struct IMGPAIR *)GMT_memory(CNULL, (size_t)n_alloc, sizeof(struct IMGPAIR), "imgbm_combine");
	n_in = 0;

	while ( (fread((char *)&imgpair[n_in], sizeof(struct IMGPAIR), 1, stdin)) == 1) {
		n_in++;
		if (n_in == n_alloc) {
			n_alloc += GMT_CHUNK;
			imgpair = (struct IMGPAIR *)GMT_memory((char *)imgpair, (size_t)n_alloc, sizeof(struct IMGPAIR), "imgbm_combine");
		}
	}

	qsort((char *)imgpair, n_in, sizeof(struct IMGPAIR), compare_imgpairs);

	i = 0;
	n_out = 0;
	while (i < n_in) {
		j = i + 1;
		while ( (j < n_in) && (imgpair[j].index == imgpair[i].index) ) j++;
		n_this_box = j - i;
		n_half = n_this_box/2;
		if (n_this_box%2) {
			top = imgpair[i+n_half].depth;
		}
		else {
			top = (int)rint(0.5*(imgpair[i+n_half].depth + imgpair[i+n_half-1].depth));
		}
		imgpair[i].depth = top;
		if ( (fwrite((char *)&imgpair[i], sizeof(struct IMGPAIR), 1, stdout)) != 1) {
			fprintf(stderr,"imgbm_combine:  Write error.\n");
			exit(-1);
		}
		n_out++;
		i = j;
	}
	free((char *)imgpair);
	fprintf(stderr,"imgbm_combine read %d wrote %d\n", n_in, n_out);
	exit(0);
}

int	compare_imgpairs(p1, p2)
struct IMGPAIR *p1, *p2;
{
	if (p1->index < p2->index)
		return(-1);
	else if (p1->index > p2->index)
		return(1);
	else if (p1->depth < p2->depth)
		return(-1);
	else  if (p1->depth > p2->depth)
		return(1);
	else
		return(0);
}