#
#include <stdio.h>
#include <stdint.h>
#include "defs.h"
#include <math.h>
struct Site *nextone();

int triangulate, sorted, plot, debug;
double xmin, xmax, ymin, ymax, deltax, deltay;
int vert_count, edge_count, line_count, dual_count, site_count;
double *site_list, *vert_list, *line_list;
unsigned int *edge_list, *dual_list;

/* sort sites on y, then x, coord */
int scomp(s1, s2) struct Point *s1, *s2;
{
	if (s1->y < s2->y)
		return (-1);
	if (s1->y > s2->y)
		return (1);
	if (s1->x < s2->x)
		return (-1);
	if (s1->x > s2->x)
		return (1);
	return (0);
}

/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
void readsites(double *lattice) {

	sites = (struct Site *)myalloc(nsites * sizeof *sites);
	for (unsigned int j = 0; j < nsites; j++) {
		sites[j].coord.x = lattice[2 * j];
		sites[j].coord.y = lattice[2 * j + 1];
		sites[j].sitenbr = j;
		sites[j].refcnt = 0;
	}
	qsort(sites, nsites, sizeof(struct Site), scomp);
	xmin = sites[0].coord.x;
	xmax = sites[0].coord.x;
	for (unsigned int i = 1; i < nsites; i += 1) {
		if (sites[i].coord.x < xmin)
			xmin = sites[i].coord.x;
		if (sites[i].coord.x > xmax)
			xmax = sites[i].coord.x;
	};
	ymin = sites[0].coord.y;
	ymax = sites[nsites - 1].coord.y;
}

unsigned int delete_duplicates(unsigned int ne, unsigned int *etv) {
	unsigned int ndup = 0;
	bool *duplicates = (bool *)calloc(ne, sizeof(bool));
	for (unsigned int i = 0; i < ne; i++) {
		unsigned int v1 = etv[2 * i];
		unsigned int v2 = etv[2 * i + 1];
		for (unsigned int j = 0; j < i; j++) {
			unsigned int w1 = etv[2 * j];
			unsigned int w2 = etv[2 * j + 1];
			if (w1 == v1 && w2 == v2) {
				duplicates[j] = true;
				ndup++;
				break;
			}
		}
	}
	unsigned int newsize = (int)ne - (int)ndup;
	unsigned int count = 0;
	for (unsigned int i = 0; i < ne; i++) {
		if (!duplicates[i]) {
			etv[2 * count] = etv[2 * i];
			etv[2 * count + 1] = etv[2 * i + 1];
			count++;
		}
	}
	free(duplicates);
	return newsize;
}

intptr_t *run_voronoi(unsigned int num, double *lattice, bool periodic, double xmin, double xmax, double ymin, double ymax)
{
	struct Site *(*next)();

	sorted = 0;
	triangulate = 0;
	plot = 0;
	debug = 0;

	alloclist = (void **)malloc(9 * num * sizeof(void *));
	allocnum = 0;

	freeinit(&sfl, sizeof *sites);

	unsigned int eff_num = num;
	double *eff_lattice = lattice;

	if (periodic) {
		eff_num = 9 * num;
		eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));

		for (unsigned int i = 0; i < num; i++) {
			// original sites - our baby boys
			eff_lattice[2*i] = lattice[2*i];
			eff_lattice[2*i+1] = lattice[2*i+1];

			// sites to the right
			eff_lattice[2*(num+i)+1] = lattice[2*i+1] - xmin + xmax;
			eff_lattice[2*(num+i)] = lattice[2*i];

			// sites to the left
			eff_lattice[2*(2*num+i)+1] = lattice[2*i+1] - xmax + xmin;
			eff_lattice[2*(2*num+i)] = lattice[2*i];

			// sites to the top
			eff_lattice[2*(3*num+i)+1] = lattice[2*i+1];
			eff_lattice[2*(3*num+i)] = lattice[2*i] - ymin + ymax;

			// sites to the bottom
			eff_lattice[2*(4*num+i)+1] = lattice[2*i+1];
			eff_lattice[2*(4*num+i)] = lattice[2*i] - ymax + ymin;

			// sites to the upper right
			eff_lattice[2*(5*num+i)+1] = lattice[2*i+1] - xmin + xmax;
			eff_lattice[2*(5*num+i)] = lattice[2*i] - ymin + ymax;

			// sites to the upper left
			eff_lattice[2*(6*num+i)+1] = lattice[2*i+1] - xmax + xmin;
			eff_lattice[2*(6*num+i)] = lattice[2*i] - ymin + ymax;

			// sites to the lower left
			eff_lattice[2*(7*num+i)+1] = lattice[2*i+1] - xmax + xmin;
			eff_lattice[2*(7*num+i)] = lattice[2*i] + ymin - ymax;

			// sites to the lower right
			eff_lattice[2*(8*num+i)+1] = lattice[2*i+1] + xmax - xmin;
			eff_lattice[2*(8*num+i)] = lattice[2*i] + ymin - ymax;
		}
	}

	nsites = eff_num;
	readsites(eff_lattice);
	if (periodic) free(eff_lattice);
	next = nextone;

	siteidx = 0;
	geominit();

	site_list = malloc(2 * nsites * sizeof(double));
	vert_list = NULL;
	line_list = NULL;
	edge_list = NULL;
	dual_list = NULL;
	site_count = 0;
	vert_count = 0;
	edge_count = 0;
	line_count = 0;
	dual_count = 0;

	voronoi(triangulate, next);

	unsigned int real_vert_count = vert_count;
	unsigned int real_edge_count = edge_count;
	double *real_vert_list = vert_list;
	unsigned int *real_edge_list = edge_list;
	unsigned int *real_dual_list = dual_list;

	if (periodic) {
		real_vert_count = 0;
		real_edge_count = 0;
		real_vert_list = (double *)malloc(2 * vert_count * sizeof(double));
		real_edge_list = (unsigned int *)malloc(2 * edge_count * sizeof(unsigned int));
		real_dual_list = (unsigned int *)malloc(3 * dual_count * sizeof(unsigned int));
		unsigned int *triangle_analyzed = (unsigned int *)malloc(4 * dual_count * sizeof(unsigned int));
		unsigned int q = 0;
		int *new_vert_ind = (int *)malloc(vert_count * sizeof(int));

		for (unsigned int i = 0; i < vert_count; i++) {
			unsigned int t1, t2, t3;
			t1 = dual_list[3*i]; t2 = dual_list[3*i+1]; t3 = dual_list[3*i+2];
			if (t1 >= num && t2 >= num && t3 >= num) {
				new_vert_ind[i] = -1;
			}
			else if (t1 < num && t2 < num && t3 < num) {
				real_vert_list[2*real_vert_count] = vert_list[2*i];
				real_vert_list[2*real_vert_count+1] = vert_list[2*i+1];
				real_dual_list[3*real_vert_count] = t1;
				real_dual_list[3*real_vert_count+1] = t2;
				real_dual_list[3*real_vert_count+2] = t3;
				new_vert_ind[i] = real_vert_count;
				real_vert_count++;
			}
			else {
				unsigned int tt1, tt2, tt3;
				tt1 = t1 % num; tt2 = t2 % num; tt3 = t3 % num;
				unsigned int s1, s2, s3;
				s1 = tt1 < tt2 ? (tt1 < tt3 ? tt1 : tt3) : (tt2 < tt3 ? tt2 : tt3);
				s3 = tt1 > tt2 ? (tt1 > tt3 ? tt1 : tt3) : (tt2 > tt3 ? tt2 : tt3);
				s2 = tt1 < tt2 ? (tt1 > tt3 ? tt1 : (tt2 < tt3 ? tt2 : tt3)) : (tt1 < tt3 ? tt1 : (tt2 < tt3 ? tt3: tt2));

				bool matched = false;
				unsigned int ass_vert;

				for (unsigned int j = 0; j < q; j++) {
					if (s1 == triangle_analyzed[4*j] && s2 == triangle_analyzed[4*j+1] && s3 == triangle_analyzed[4*j+2]) {
						matched = true;
						ass_vert = triangle_analyzed[4*j+3];
						break;
					}
				}

				if (matched) {
					new_vert_ind[i] = ass_vert;
				} else {
					triangle_analyzed[4*q] = s1;
					triangle_analyzed[4*q+1] = s2;
					triangle_analyzed[4*q+2] = s3;
					triangle_analyzed[4*q+3] = real_vert_count;
					q++;
					real_vert_list[2*real_vert_count+1] = fmod(2*xmax+vert_list[2*i+1], xmax);
					real_vert_list[2*real_vert_count] = fmod(2*ymax+vert_list[2*i], ymax);
					real_dual_list[3*real_vert_count] = t1 % num;
					real_dual_list[3*real_vert_count+1] = t2 % num;
					real_dual_list[3*real_vert_count+2] = t3 % num;
					new_vert_ind[i] = real_vert_count;
					real_vert_count++;
				}
			}
		}
		for (unsigned int i = 0; i < edge_count; i++) {
			unsigned int v1, v2;
			v1 = edge_list[2 * i];
			v2 = edge_list[2 * i+1];
			if (v1 != UINT_MAX && v2 != UINT_MAX) {
				if (new_vert_ind[v1] >= 0 && new_vert_ind[v2] >=0) {
					unsigned int nv1 = new_vert_ind[v1];
					unsigned int nv2 = new_vert_ind[v2];
					real_edge_list[2*real_edge_count] = nv1 < nv2 ? nv1 : nv2;
					real_edge_list[2*real_edge_count+1] = nv1 < nv2 ? nv2 : nv1;
					real_edge_count++;
				}
			}
		}
		unsigned int new_edge_count = delete_duplicates(real_edge_count, real_edge_list);
		real_edge_count = new_edge_count;
		free(triangle_analyzed);
		free(new_vert_ind);
		free(dual_list);
		free(edge_list);
		free(vert_list);
	}


	intptr_t *output = (intptr_t *)malloc(6 * sizeof(intptr_t));
	output[0] = (intptr_t)malloc(4 * sizeof(unsigned int));
	((unsigned int *)output[0])[0] = real_vert_count;
	((unsigned int *)output[0])[1] = real_edge_count;
	((unsigned int *)output[0])[2] = line_count;
	((unsigned int *)output[0])[3] = real_vert_count;
	output[1] = (intptr_t)site_list;
	output[2] = (intptr_t)real_vert_list;
	output[3] = (intptr_t)real_edge_list;
	output[4] = (intptr_t)line_list;
	output[5] = (intptr_t)real_dual_list;

	free(ELhash);
	free(PQhash);
	free(sites);
	for (unsigned int i = 0; i < allocnum; i++) {
		free(alloclist[i]);
	}
	free(alloclist);

	return output;
}

/* return a single in-storage site */
struct Site *nextone() {
	for (; siteidx < nsites; siteidx += 1) {
		if (siteidx == 0 || sites[siteidx].coord.x != sites[siteidx - 1].coord.x ||
				sites[siteidx].coord.y != sites[siteidx - 1].coord.y) {
			siteidx += 1;
			return (&sites[siteidx - 1]);
		};
	};
	return ((struct Site *)NULL);
}