/*
# This file is part of the Astrometry.net suite.
# Licensed under a 3-clause BSD style license - see LICENSE
*/
#include <string.h>
#include <math.h>
#include <assert.h>

#include "plotindex.h"
#include "cairoutils.h"
#include "ioutils.h"
#include "log.h"
#include "errors.h"
#include "starutil.h"
#include "index.h"
#include "qidxfile.h"
#include "permutedsort.h"

DEFINE_PLOTTER(index);

plotindex_t* plot_index_get(plot_args_t* pargs) {
	return plotstuff_get_config(pargs, "index");
}

void* plot_index_init(plot_args_t* plotargs) {
	plotindex_t* args = calloc(1, sizeof(plotindex_t));
	args->indexes = pl_new(16);
	args->qidxes = pl_new(16);
	args->stars = TRUE;
	args->quads = TRUE;
	args->fill = FALSE;
	return args;
}

static void pad_qidxes(plotindex_t* args) {
	while ((pl_size(args->qidxes)) < pl_size(args->indexes))
		pl_append(args->qidxes, NULL);
}

void plot_quad_xy(cairo_t* cairo, double* quadxy, int dimquads) {
	int k;
	double cx, cy;
	double theta[DQMAX];
	int* perm;

	cx = cy = 0.0;
	for (k=0; k<dimquads; k++) {
		cx += quadxy[2*k+0];
		cy += quadxy[2*k+1];
	}
	cx /= dimquads;
	cy /= dimquads;

	for (k=0; k<dimquads; k++)
		theta[k] = atan2(quadxy[2*k+1] - cy, quadxy[2*k+0] - cx);
	perm = permuted_sort(theta, sizeof(double), compare_doubles_asc, NULL, dimquads);
	for (k=0; k<dimquads; k++) {
		double px,py;
		px = quadxy[2 * perm[k] + 0];
		py = quadxy[2 * perm[k] + 1];
		if (k == 0) {
			cairo_move_to(cairo, px, py);
		} else {
			cairo_line_to(cairo, px, py);
		}
	}
	free(perm);
	cairo_close_path(cairo);
}

static void plotquad(cairo_t* cairo, plot_args_t* pargs, plotindex_t* args, index_t* index, int quadnum, int DQ) {
	int k;
	unsigned int stars[DQMAX];
	double ra, dec;
	double px, py;
	double xy[DQMAX*2];
	int N;

	quadfile_get_stars(index->quads, quadnum, stars);
	N = 0;
	for (k=0; k<DQ; k++) {
		if (startree_get_radec(index->starkd, stars[k], &ra, &dec)) {
			ERROR("Failed to get RA,Dec for star %i\n", stars[k]);
			continue;
		}
		if (!plotstuff_radec2xy(pargs, ra, dec, &px, &py)) {
			ERROR("Failed to convert RA,Dec %g,%g to pixels for quad %i\n", ra, dec, quadnum);
			continue;
		}
		xy[2*k + 0] = px;
		xy[2*k + 1] = py;
		N++;
	}
	if (N < 3)
		return;
	plot_quad_xy(cairo, xy, N);
	if (args->fill)
		cairo_fill(cairo);
	else
		cairo_stroke(cairo);
}

void plot_index_plotquad(cairo_t* cairo, plot_args_t* pargs, plotindex_t* args, index_t* index, int quadnum, int DQ) {
	plotquad(cairo, pargs, args, index, quadnum, DQ);
}

int plot_index_plot(const char* command,
					cairo_t* cairo, plot_args_t* pargs, void* baton) {
	plotindex_t* args = (plotindex_t*)baton;
	int i;
	double ra, dec, radius;
	double xyz[3];
	double r2;

	pad_qidxes(args);

	plotstuff_builtin_apply(cairo, pargs);

	if (plotstuff_get_radec_center_and_radius(pargs, &ra, &dec, &radius)) {
		ERROR("Failed to get RA,Dec center and radius");
		return -1;
	}
	radecdeg2xyzarr(ra, dec, xyz);
	r2 = deg2distsq(radius);
	logmsg("Field RA,Dec,radius = (%g,%g), %g deg\n", ra, dec, radius);
	logmsg("distsq: %g\n", r2);

	for (i=0; i<pl_size(args->indexes); i++) {
		index_t* index = pl_get(args->indexes, i);
		int j, N;
		int DQ;
		double px,py;

		if (args->stars) {
			// plot stars
			double* radecs = NULL;
			startree_search_for(index->starkd, xyz, r2, NULL, &radecs, NULL, &N);
			if (N) {
				assert(radecs);
			}
			logmsg("Found %i stars in range in index %s\n", N, index->indexname);
			for (j=0; j<N; j++) {
				logverb("  RA,Dec (%g,%g) -> x,y (%g,%g)\n", radecs[2*j], radecs[2*j+1], px, py);
				if (!plotstuff_radec2xy(pargs, radecs[j*2], radecs[j*2+1], &px, &py)) {
					ERROR("Failed to convert RA,Dec %g,%g to pixels\n", radecs[j*2], radecs[j*2+1]);
					continue;
				}
				cairoutils_draw_marker(cairo, pargs->marker, px, py, pargs->markersize);
				cairo_stroke(cairo);
			}
			free(radecs);
		}
		if (args->quads) {
			DQ = index_get_quad_dim(index);
			qidxfile* qidx = pl_get(args->qidxes, i);
			if (qidx) {
				int* stars;
				int Nstars;
				il* quadlist = il_new(256);

				// find stars in range.
				startree_search_for(index->starkd, xyz, r2, NULL, NULL, &stars, &Nstars);
				logmsg("Found %i stars in range of index %s\n", N, index->indexname);
				logmsg("Using qidx file.\n");
				// find quads that each star is a member of.
				for (j=0; j<Nstars; j++) {
					uint32_t* quads;
					int Nquads;
					int k;
					if (qidxfile_get_quads(qidx, stars[j], &quads, &Nquads)) {
						ERROR("Failed to get quads for star %i\n", stars[j]);
						return -1;
					}
					for (k=0; k<Nquads; k++)
						il_insert_unique_ascending(quadlist, quads[k]);
				}
				for (j=0; j<il_size(quadlist); j++) {
					plotquad(cairo, pargs, args, index, il_get(quadlist, j), DQ);
				}

			} else {
				// plot quads
				N = index_nquads(index);
				for (j=0; j<N; j++) {
					plotquad(cairo, pargs, args, index, j, DQ);
				}
			}
		}
	}
	return 0;
}

int plot_index_add_qidx_file(plotindex_t* args, const char* fn) {
	int i;
	qidxfile* qidx = qidxfile_open(fn);
	if (!qidx) {
		ERROR("Failed to open quad index file \"%s\"", fn);
		return -1;
	}
	pad_qidxes(args);
	i = pl_size(args->indexes) - 1;
	pl_set(args->qidxes, i, qidx);
	return 0;
}

int plot_index_add_file(plotindex_t* args, const char* fn) {
	index_t* index = index_load(fn, 0, NULL);
	if (!index) {
		ERROR("Failed to open index \"%s\"", fn);
		return -1;
	}
	pl_append(args->indexes, index);
	return 0;
}

int plot_index_command(const char* cmd, const char* cmdargs,
					   plot_args_t* pargs, void* baton) {
	plotindex_t* args = (plotindex_t*)baton;
	if (streq(cmd, "index_file")) {
		const char* fn = cmdargs;
		return plot_index_add_file(args, fn);
	} else if (streq(cmd, "index_qidxfile")) {
		const char* fn = cmdargs;
		return plot_index_add_qidx_file(args, fn);
	} else if (streq(cmd, "index_draw_stars")) {
		args->stars = atoi(cmdargs);
	} else if (streq(cmd, "index_draw_quads")) {
		args->quads = atoi(cmdargs);
	} else if (streq(cmd, "index_fill")) {
		args->fill = atoi(cmdargs);
	} else {
		ERROR("Did not understand command \"%s\"", cmd);
		return -1;
	}
	return 0;
}

void plot_index_free(plot_args_t* plotargs, void* baton) {
	plotindex_t* args = (plotindex_t*)baton;
	int i;
	for (i=0; i<pl_size(args->indexes); i++) {
		index_t* index = pl_get(args->indexes, i);
		index_free(index);
	}
	pl_free(args->indexes);
	for (i=0; i<pl_size(args->qidxes); i++) {
		qidxfile* qidx = pl_get(args->qidxes, i);
		qidxfile_close(qidx);
	}
	pl_free(args->qidxes);
	free(args);
}