/*
# This file is part of the Astrometry.net suite.
# Licensed under a 3-clause BSD style license - see LICENSE
*/


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

#include "healpix.h"
#include "mathutil.h"
#include "pnpoly.h"

/* Size of the line buffer */
#define BUFSIZE 1000

/* Need to extend this structure to handle more complicated polygons */
typedef struct {
	double *corners;
} rect_field;

/* Dead simple linked list */
struct ll_node {
	int data;
	struct ll_node *next;
};

/* TODO: factor this out into mathutil or something similar. Eventually
 * this file won't need it anymore anyway if the healpix code is rewritten
 * to abandon the current behaviour of using hierarchical numbering for 
 * powers of 4. */
static Inline anbool ispowerof4(unsigned int x) {
	if (x >= 0x40000)
		return (					x == 0x40000   ||
				 x == 0x100000   || x == 0x400000  ||
				 x == 0x1000000  || x == 0x4000000 ||
				 x == 0x10000000 || x == 0x40000000);
	else
		return (x == 0x1		|| x == 0x4	   ||
				x == 0x10	   || x == 0x40	  ||
				x == 0x100	  || x == 0x400	 ||
				x == 0x1000	 || x == 0x4000	||
				x == 0x10000);
}

/* Compute the min of n doubles */
static double nmin(double *p, int n)
{
	int i;
	double cmin;
	assert (n > 0);
	cmin = p[0];
	for (i = 1; i < n; i++)
	{   
		if (p[i] < cmin)
		cmin = p[i];
	}
	return cmin;
}

/* Compute the min of n doubles and return the index rather than value */
static int nminind(double *p, int n)
{
	int i;
	int cminind;
	
	assert (n > 0);
	cminind = 0;
	for (i = 1; i < n; i++)
	{   
		if (p[i] < p[cminind])
		{   
			cminind = i;
		}
	}
	return cminind;
}

/* Compute the max of n doubles */
static double nmax(double *p, int n)
{
	int i;
	double cmax;

	assert (n > 0);
	cmax = p[0];
	for (i = 1; i < n; i++)
	{   
		if (p[i] > cmax)
			cmax = p[i];
	}
	return cmax;															  
}

/* Compute the range from max to min of n points stored in a double vector. */
static double range(double *p, int n)
{
	double max = nmax (p, n);
	double min = nmin (p, n);
	return fabs (max - min);
}


/* Wraps around point_in_poly to choose the axis along which the corners 
 * have minimum variance in order to do the chop-off projection.
 *
 * TODO: generalize this to n-point convex polygons. Notes have been made
 * as to where changes need to occur and what.
 */
int is_inside_field(rect_field *f, double *p)
{
	// DIMQUADS
	/* FIXME assumes 4 points, will need malloc'ing */
	double coordtrans[12];
	
	double *coord1;
	double *coord2;

	/* FIXME assumes 4 points */
	double *x = coordtrans;
	double *y = coordtrans + 4;
	double *z = coordtrans + 8;
	double t_coord1, t_coord2;
	double ranges[3];
	int i;
	int xgot = 0, ygot = 0, zgot = 0;
	/* FIXME Assumes 4 points */
	for (i = 0; i < 12; i++)
	{
		int j = i % 3;
		switch (j)
		{
			case 0:
				coordtrans[xgot++] = f->corners[i];
				break;
			case 1:
				/* FIXME 4+ should be n+ */
				coordtrans[4 + ygot++] = f->corners[i];
				break;
			case 2:
				/* FIXME 4+ should be n+ */
				coordtrans[8 + zgot++] = f->corners[i];
				break;
		}
	}
	/* FIXME all these calls should use n */
	ranges[0] = range(x, 4);
	ranges[1] = range(y, 4);
	ranges[2] = range(z, 4);
	
	switch (nminind(ranges, 3))
	{
		case 0:
			coord1 = y;
			coord2 = z;
			t_coord1 = p[1];
			t_coord2 = p[2];
			break;
		case 1:
			coord1 = x;
			coord2 = z;
			t_coord1 = p[0];
			t_coord2 = p[2];
			break;
		case 2:
			coord1 = x;
			coord2 = y;
			t_coord1 = p[0];
			t_coord2 = p[1];
			break;
		default:
			return -1;
	}
	/* FIXME should call with n rather than 4 */
	return point_in_poly(coord1, coord2, 4, t_coord1, t_coord2);

}
/* An ugly recursive implementation written while still trying to get 
 * the algorithm right, left here for... some reason */
#if 0
void fill_maps_recursive(char *minmap, char *maxmap, uint hpx, uint Nside,
		rect_field *curfield, char *visited)
{
	double thishpx_coords[3];
	uint neighbours[8];
	uint nn;
	if (visited[hpx / 8] & (1 << (hpx % 8)))
		return;
	visited[hpx / 8] |= 1 << (hpx % 8);
	healpix_to_xyzarr_lex(0.5, 0.5, hpx, Nside, thishpx_coords);
	//printf("Examining healpix %d, centered at (%f, %f, %f)\n", hpx, 
	//		thishpx_coords[0], thishpx_coords[1], thishpx_coords[2]);
	if (is_inside_field(curfield, thishpx_coords))
	{
		int j;
		maxmap[hpx / 8] |= (1 << (hpx % 8));
		nn = healpix_get_neighbours_nside(hpx, neighbours, Nside);
		for (j = 0; j < nn; j++)
		{
				
			double ncoords[3];
			healpix_to_xyzarr_lex(0.5, 0.5, neighbours[j], Nside, ncoords);

			//printf("- Examining neighbour healpix %d, centered at (%f, %f, %f)\n", neighbours[j],
			//	ncoords[0], ncoords[1], ncoords[2]);
			
			
			if (!is_inside_field(curfield, ncoords)) {
				//printf("-- Not in field, breaking off neighbour search\n");
				break;
			}
		}
		if (j == nn)
			minmap[(hpx / 8)] |= (1 << (hpx % 8));
		for (j = 0; j < nn; j++)
		{
			//printf("Recursing on neighbour of %d, %d\n", hpx, neighbours[j]);
			fill_maps_recursive(minmap, maxmap, neighbours[j], Nside,
					curfield, visited);
		}
	}
}
void fill_maps(char *minmap, char *maxmap, uint hpx, uint Nside,
		rect_field *curfield)
{
	char *visited = malloc(2 * Nside * Nside * sizeof(char));
	uint i;
	uint visitedcnt = 0;
	for (i = 0; i < 2 * Nside * Nside; i++)
		visited[i] = 0;

	fill_maps_recursive(minmap, maxmap, hpx, Nside, curfield, visited);
}
#endif


/* This basically takes two appropriately sized char[] arrays to be used
 * as bitmaps for the min/max healpix mapping of the sky, as well as a 
 * field we'd like to process and add to the bitmaps, and a healpix number
 * to start adding from, and an Nside factor.
 *
 * This basically starts at the healpix closest to the center of a field 
 * (actually it starts at hpx, see todo below) and percolates outward by
 * checking neighbouring healpixes. A healpix is included in the upper bound
 * map if it's center is inside the field boundaries, and included in the
 * lower bound map only if all of its neighbours' centers are inside the
 * field boundaries as well.
 * 
 * TODO: bring the center-computing code from main() up here
 * TODO: factor that bitmap access ugliness out into some damned macros
 */
void fill_maps(char *minmap, char *maxmap, uint hpx, uint Nside,
				rect_field *curfield)

{
	/* Gotta love "are we done" switch variables */
	anbool done = FALSE;
	
	/* Bitmap we'll use to keep from revisiting the same healpixes */ 
	char *visited = malloc(2 * Nside * Nside * sizeof(char));
	
	/* Store the head of the queue (actually a LIFO) of healpixes to examine */
	struct ll_node *queue = NULL;
	double thishpx_coords[3];

	int i;

	/* Initialize the visited map */	
	for (i = 0; i < 2 * Nside * Nside; i++)
		visited[i] = 0;

	do {
		/* nn = "number of neighbours */
		uint nn;
		anbool found_neighbour_outside = FALSE;
		
		/* always need room for at least 8 neighbours; actual number
		 * gets stored in nn */
		uint neighbours[8];
		if (visited[hpx / 8] & (1 << (hpx % 8)))
		{
			/* should never happen */
			assert(1 == 0);
		}
		/* set that we've visited hpx */
		visited[hpx / 8] |= (1 << (hpx % 8));
		
		/* compute the xyz location of the center of hpx */
		healpix_to_xyzarr(hpx, Nside, 0.5, 0.5, thishpx_coords);
		
		/* skip the body of this loop if we can */
		if (!is_inside_field(curfield, thishpx_coords))
			goto getnext;
		
		/* always include in maxmap */
		maxmap[hpx / 8] |= 1 << (hpx % 8);

		nn = healpix_get_neighbours(hpx, neighbours, Nside);

		/* check inclusion for each neighbour and enqueue it if unvisited */
		for (i = 0; i < nn; i++)
		{
			double ncoords[3];
			healpix_to_xyzarr(neighbours[i], Nside, 0.5, 0.5, ncoords);
			
			if (!is_inside_field(curfield, ncoords))
				found_neighbour_outside = TRUE;
			
			if (!(visited[neighbours[i] / 8] & (1 << (neighbours[i] % 8))))
			{
				struct ll_node *newnode = malloc(sizeof(struct ll_node));
				newnode->next = queue;
				newnode->data = neighbours[i];
				queue = newnode;
			}
		}

		/* If no neighbours lie outside then include in the minmap */
		if (!found_neighbour_outside)
			minmap[hpx / 8] |= 1 << (hpx % 8);
		
		/* dequeue the next one to be checked and set hpx */
		getnext:
		if (queue != NULL)
		{
			struct ll_node *newhead = queue->next;
			hpx = queue->data;
			free(queue);
			queue = newhead;
		}
		else {
			/* If the queue is empty and we've gone through the last iteration
			 * then it's okay to terminate. note that right after dequeueing
			 * you can have queue == NULL but still not be done, which is 
			 * why we can't just use "queue != NULL" as our loop condition */
			done = TRUE;
		}
	} while (!done);
	
	free(visited);
}

static void print_help(FILE *f, char *name)
{
	fprintf(f, "This program computes statistics about a set of (for the moment) rectangular\nfields on the sky.\n\n");
	fprintf(f, "\tUsage: %s -N <Nside> [-I <inputfile>]\n\n", name);
	fprintf(f, "In the absence of a -I switch, reads coordinates from standard input.\n");
	fprintf(f, "Input should be 4 XYZ coordinates per line, with components and coordinates\nseparated by tabs (i.e. 12 tab-delimited doubles)\n");
}

int main(int argc, char **argv)
{
	double max;
	rect_field curfield;
	int filled_min = 0, filled_max = 0;
	char *hpmap_min, *hpmap_max;
	char *buf = malloc(BUFSIZE * sizeof(char));
	int ich, i;
	int Nside = -1;
	uint fields;
	char *infilename = NULL;
	FILE *input;
	if (argc == 0)
	{
		print_help(stderr, argv[0]);
		exit(1);
	}
	while ((ich = getopt(argc, argv, "N:I:")) != EOF)
	{
		switch (ich)
		{
			case 'N':
				Nside = atoi(optarg);
				break;
			case 'I':
				if (optarg == NULL)
				{
					print_help(stderr, argv[0]);
					fprintf(stderr, "Error: -I requires argument");
					exit(1);
				}
				infilename = strdup(optarg);
				break;
		}
	}
	if (Nside < 1)
	{
		print_help(stderr, argv[0]);
		fprintf(stderr, "\nError: specify a positive Nside value with -N\n");
		exit(1);
	}
	else if (ispowerof4(Nside))
	{
		print_help(stderr, argv[0]);
		fprintf(stderr, "Error: Nside values that are powers of 4 \
				are bad news at\n the moment, choose a different one\n");
		exit(1);
	}
	
	if (infilename) {
		input = fopen(infilename, "r");
		if (input == NULL)
		{
			perror(argv[0]);
			exit(1);
		}
	}
	else {
		input = stdin;
	}
	/* We could get away with allocating ceil(2/3 * Nside * Nside) */
	hpmap_min = malloc(2 * Nside * Nside * sizeof(char));
	hpmap_max = malloc(2 * Nside * Nside * sizeof(char));
	
	if (hpmap_min == NULL || hpmap_max == NULL)
	{
		fprintf(stderr, "malloc failed!\n");
		exit(1);
	}	
	for (i = 0; i < 2 * Nside * Nside; i++) {
		hpmap_min[i] = 0;
		hpmap_max[i] = 0;
	}
	
	curfield.corners = malloc(3 * 4 * sizeof(double));
	fields = 0;
	while (fgets(buf, BUFSIZE, input) != NULL)
	{
		uint centerhp;
		int i, j;
		double center[3];
		center[0] = center[1] = center[2] = 0;
		fields++;
		//printf("Doing field %d\n",fields);
		curfield.corners[0] = atof(strtok(buf, "\t"));
		
		/* 12 = 3 coords x 4 pts, got 1 */
		for (j = 1; j < 12; j++)
		{
			char *tok = strtok(NULL, "\t");
			if (tok == NULL)
			{
				fprintf(stderr, "Premature end of line!\n");
				exit(1);
			}
			curfield.corners[j] = atof(tok);
		}
		
		for (i = 0; i < 4; i++)
		{
			for (j = 0; j < 3; j++)
			{
				center[j] += curfield.corners[3*i + j];
			}
		}
		for (i = 0; i < 3; i++)
			center[i] /= 4;
		normalize_3(center);
		
		centerhp = xyzarrtohealpix(center, (uint)Nside);
		fill_maps(hpmap_min, hpmap_max, centerhp, (uint)Nside, &curfield);
	}
	for (i = 0; i < 12 * Nside * Nside; i++)
	{
		if (hpmap_min[i / 8] & (1 << (i % 8)))
			filled_min++;
		if (hpmap_max[i / 8] & (1 << (i % 8)))
			filled_max++;
	}
	max = 12 * Nside * Nside;

	printf("Min: %f, Max: %f\n",
			((double)filled_min) / max, 
			((double)filled_max) / max);

	fclose(input);
	return 0;
}