#define NX 10800
#define NY 6336
#define NPAIR 10000000

#include <stdio.h>
#include <math.h>

struct IMGPAIR {
        int     index;
        int     depth;
} imgpair[NPAIR];

main(argc, argv)
int     argc;
char    **argv;
{
        
        int     i, j, k, n_in, n_out, n_this_box, n_half, top, compare_imgpairs();
        double  radlat, radius, ypix, x,y,z;
        char    buffer[128];

        radius = (NX/2)/M_PI;


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

        k = 0;
	while(gets(buffer)) {
		if ( (sscanf(buffer, "%lf %lf %lf", &x, &y, &z)) != 3) {
			fprintf(stderr,"xyz2imgbm:  error reading x,y,z from line %d\n", k+1);
			exit(-1);
		}
		if (fabs(y) >= 90.0) {
			fprintf(stderr,"xyz2imgbm:  y range error at line %d\n", k+1);
			exit(-1);
		}
		if (z < -12000.0 || z > 9000.) {
			fprintf(stderr,"xyz2imgbm:  z range error at line %d\n", k+1);
			exit(-1);
		}
		if (k == NPAIR) {
			fprintf(stderr,"xyz2imgbm:  NPAIR is too small.  Recompile.\n");
			exit(-1);
		}
		while (x < 0.0) x += 360.0;
		while (x >= 360.0) x -= 360.0;
		i = floor(30.0 * x);
		if (i == NX) i = 0;
		radlat = M_PI * y / 180.0;
                ypix = NY/2 - radius*log(tan(M_PI_4 + 0.5*radlat));
                if (ypix < 0.0 || ypix >= (double)NY) continue;
                j = floor(ypix);
                imgpair[k].index = j * NX + i;
                imgpair[k].depth = (int)rint(z);
                k++;
        }

        n_in = k;
        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,"xyz2imgbm:  Write error.\n");
                        exit(-1);
                }
                n_out++;
                i = j;
        }
        free((char *)imgpair);
        fprintf(stderr,"xyz2imgbm 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);
}