#include <delaunay.h>
#include <delaunay_plot.h>
#include <delaunay_debug.h>

#include <quad.h>
#include <affirm.h>
#include <pswr.h>
#include <jsrandom.h>

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

#define NORMAL 1
#define NTOSS 9
#define VERBOSE 1
#define RBIG 10.0
#define PLOTDIM 10000

#define MC (delaunay_MAX_COORD)

typedef struct quadArcList{
	quad_arc_t e;
	struct quadArcList *next;
	struct quadArcList *prev;
}*quad_arc_list;

delaunay_site_t* divide_edge(quad_arc_t e, int color);

void divide_triangle1(quad_arc_t e);

void divide_triangle2(quad_arc_t e);

char execute_match(quad_arc_t e);

void execute_division(quad_arc_t e, char type);

quad_arc_list make_edge_list(quad_arc_t e, int nsites);

void make_subdivision(quad_arc_list l);

delaunay_site_t *makesites(int nsites, bool_t normal);

void plot_triangulation (quad_arc_t e, delaunay_site_t *st, int nsites, char *prefix, bool_t eps);

void plot_triangulation_debug (quad_arc_t e, delaunay_site_t *st, int nsites, char *prefix, bool_t eps, quad_arc_list list);

int main(int argc, char **argv)
  { 
    int nsites = atoi(argv[1]);
    assert(nsites > 2);
    bool_t normal = 1;
    bool_t eps = 1;
    quad_arc_t e;
    quad_arc_list list;
    fprintf(stderr, "Creating %d sites...\n", nsites);
    delaunay_site_t *st = makesites(nsites, normal);
    fprintf(stderr, "Building delaunay...\n");
    e = delaunay_build (st, nsites);
    fprintf(stderr, "Building the order of the triangles that will be processed...\n");
    list = make_edge_list(e, nsites);
    fprintf(stderr, "Building subdivision...\n");
    //plot_triangulation_debug(e, st, nsites, "out/debug", eps, list);
    make_subdivision(list);   
    fprintf(stderr, "Plotting delaunay...\n");
    plot_triangulation(e, st, nsites, "out/delrandom", eps);
    return(0);
  }


quad_arc_list find_quad_arc_l(quad_arc_list list, quad_arc_t e){

    if (list==NULL) return NULL;
    quad_arc_list l;
    for (l=list;l->next!=list;l=l->next){ 
	if (l->e == e) return l;
    }
    if (l->e == e) return l;
    return NULL;
}


void add_quad_arc_l(quad_arc_list *list, quad_arc_t e){

	if(*list == NULL){
		*list = notnull(malloc(sizeof(struct quadArcList)), "no mem");
		(*list)->next = (*list)->prev = (*list);
		(*list)->e = e;
	}
	else{
		quad_arc_list aux;
		aux = notnull(malloc(sizeof(struct quadArcList)), "no mem");
		aux->e = e;
		aux->next = (*list);
		aux->prev = (*list)->prev;
		(*list)->prev->next = aux;
		(*list)->prev = aux;
	}
}

void remove_quad_arc_l(quad_arc_list *begin, quad_arc_list *list){
	if(*list != NULL){
		if((*begin)==(*list)) (*begin) = (*begin)->next;
		if((*begin)==(*list)) (*begin) = NULL;
		(*list)->next->prev = (*list)->prev;
		(*list)->prev->next = (*list)->next;
		free(*list);
		*list=NULL; 
	}
}

int is_interior_edge(quad_arc_t e, quad_arc_t f){
	if (f == quad_lnext(quad_lnext(quad_lnext(f))) && f!=quad_sym(e) && f!=quad_dnext(e) && f!=quad_oprev(e))
		return 1;
	else
		return 0;
}

quad_arc_list make_edge_list(quad_arc_t e, int nsites){

	//Make a list of edges in the order that they will be called by the match function. Also set the colors of the vertices of the first triangle.

	quad_arc_list edge_list, front_list, front_list_begin, left, right;
	edge_list = front_list = front_list_begin = NULL;

	int i=0;
	quad_arc_t f,g,h;

	// set the colors of the vertices of the first triangle.
	(*ORG(e)).index = 1;
	(*DST(e)).index = 2;
	(*ORG(quad_dprev(e))).index = 3;

	//Inicialize edges from the front
	g = quad_onext(e);
	h = quad_dprev(e);
	if (is_interior_edge(e, g))
	add_quad_arc_l(&front_list_begin, g);
	if (is_interior_edge(e, h))
	add_quad_arc_l(&front_list_begin, h);
	while(front_list_begin!=NULL){
		front_list = front_list_begin;
		do{
			f = front_list->e; 
			right = find_quad_arc_l(front_list_begin,quad_sym(quad_dprev(f)));
			left = find_quad_arc_l(front_list_begin,quad_sym(quad_onext(f)));		
			if(right!=NULL){ //duas arestas do triangulo no front (na triangulacao de delaunay nao pode haver tres arestas de um triangulo no front)
				g = quad_onext(f);
				add_quad_arc_l(&edge_list, quad_sym(g));
				remove_quad_arc_l(&front_list_begin,&front_list);
				remove_quad_arc_l(&front_list_begin,&right);
				if (is_interior_edge(e, g)){
					add_quad_arc_l(&front_list_begin, g);}
				break;
			}			
			else if (left!=NULL){ //duas arestas do triangulo no front (pelo outro lado)
				g = quad_dprev(f);				
				add_quad_arc_l(&edge_list, quad_sym(g));
				remove_quad_arc_l(&front_list_begin,&front_list);
				remove_quad_arc_l(&front_list_begin,&left);
				if (is_interior_edge(e, g)){ 
					add_quad_arc_l(&front_list_begin, g);}	
				break;	
			}
			else if ((*ORG(quad_dprev(f))).index == -1){ //1 aresta do triangulo no front, vértice oposto nao pintado
				(*ORG(quad_dprev(f))).index = 0;
				g = quad_onext(f);
				h = quad_dprev(f);
				add_quad_arc_l(&edge_list, f);
				remove_quad_arc_l(&front_list_begin,&front_list);
				if (is_interior_edge(e, g)){ 
					add_quad_arc_l(&front_list_begin, g);}	
				if (is_interior_edge(e, h)){ 
					add_quad_arc_l(&front_list_begin, h);}	
				break;				
			}
			else{	//1 arestas do triangulo no front, vértice oposto ja pintado
				front_list = front_list->next;
			}
		}while (front_list!=front_list_begin);
        }
	return edge_list;
}

void make_subdivision(quad_arc_list l){
	quad_arc_list begin = l;
	char type = 'a';
	while (l!=NULL){
		type = execute_match(l->e);
		execute_division(l->e, type);
		remove_quad_arc_l(&begin,&l);
		l = begin;
	}
} 

int is_colinear(hi2_point_t a, hi2_point_t b, hi2_point_t c){
	if(a.c.c[1]*b.c.c[2]*c.c.c[0] + a.c.c[2]*b.c.c[0]*c.c.c[1] + a.c.c[0]*b.c.c[1]*c.c.c[2] == 0)
		return 1;
	else
		return 0;
}

char execute_match(quad_arc_t e){
	quad_arc_t f;
	char state[7]="000000";
	int color1 = (*ORG(e)).index;
	int color2 = -1;
	int i=1;
	if((*DST(e)).index != color1){
		color2 = (*DST(e)).index ;
		state[0]='2';
	}
	else
		state[0]='1';
	for(f=quad_lnext(e);f!=e;f=quad_lnext(f)){
		if((*DST(f)).index == color1)
			state[i]='1';
		else if((*DST(f)).index == color2)
			state[i]='2';
		else if((*DST(f)).index != 0){
			if(color2== -1){
				color2 = (*DST(f)).index;
				state[i]='2';
			}
			else
				state[i]='3';
		}
		i++;
	}
	if(state[0] == '1') //state[1]=='2'; //state[2]=='1'
		if (state[3]=='0')
			return 'd';
		else if (state[3]=='2')
			return 'h';		
		else	//state[3]='3'
			return 'g';
	else  //state[0]=='2'
		if (state[1]=='3')
			if(state[3]=='1')
				return 'e';
			else
				return 'c';
		else if (state[1]=='1')
			if(state[2]=='0')
				return 'b';
			else if(state[2]=='3')
				return 'E';
			else if(is_colinear( (*DST(e)).pt, (*DST(quad_lnext(e))).pt, (*DST(quad_lnext(quad_lnext(e)))).pt ))
				return 'f';
			else
				return 'F';
	return 'a';
}

delaunay_site_t* divide_edge(quad_arc_t e, int color){

	delaunay_site_t *e_odata = notnull(malloc(sizeof(delaunay_site_t)), "no mem");

	if((*DST(e)).pt.c.c[0] > (*ORG(e)).pt.c.c[0]){
		e_odata->pt.c.c[0] = (*DST(e)).pt.c.c[0] * 2; 
		e_odata->pt.c.c[1] = (*ORG(e)).pt.c.c[1] * (*DST(e)).pt.c.c[0] / (*ORG(e)).pt.c.c[0] + (*DST(e)).pt.c.c[1];
		e_odata->pt.c.c[2] = (*ORG(e)).pt.c.c[2] * (*DST(e)).pt.c.c[0] / (*ORG(e)).pt.c.c[0] + (*DST(e)).pt.c.c[2];
	}
	else{ 
		e_odata->pt.c.c[0] = (*ORG(e)).pt.c.c[0] * 2; 
		e_odata->pt.c.c[1] = (*ORG(e)).pt.c.c[1] + (*DST(e)).pt.c.c[1] * (*ORG(e)).pt.c.c[0] / (*DST(e)).pt.c.c[0];
		e_odata->pt.c.c[2] = (*ORG(e)).pt.c.c[2] + (*DST(e)).pt.c.c[2] * (*ORG(e)).pt.c.c[0] / (*DST(e)).pt.c.c[0];
	}

	e_odata->index = color; 

	return e_odata;

}

void divide_triangle1 (quad_arc_t e)
{
	quad_arc_t a = quad_make_edge();
	quad_arc_t b = quad_sym(quad_onext(e));
	quad_splice(b,a);
	quad_arc_t c = quad_oprev(e);
	quad_splice(c,e);
	quad_arc_t d = quad_make_edge();
	quad_splice(c,d);
	quad_splice(quad_sym(d),e);
	quad_splice(e,quad_sym(a));

	delaunay_site_t *e_odata = divide_edge(e, 6 - (*ORG(b)).index - (*DST(e)).index);

	SET_ORG(e,e_odata);	
	SET_ORG(a,ORG(b));
	SET_ORG(d,DST(b));
	SET_DST(a,ORG(e));
	SET_DST(d,ORG(e));

}

void divide_triangle2 (quad_arc_t e)
{
	quad_arc_t a = quad_make_edge();
	quad_arc_t b = quad_sym(quad_onext(e));
	quad_splice(b,a);
	quad_arc_t c = quad_sym(quad_dprev(e));
	quad_splice(c,quad_sym(a));

	SET_ORG(a,ORG(b));
	SET_DST(a,DST(e));
}
int global_c = 0;

void execute_division(quad_arc_t e, char type)
{
	quad_arc_t f,a,b,c;				global_c++;

	switch(type){
		case 'a': 
			(*ORG(quad_dprev(e))).index = 6 - (*ORG(e)).index - (*DST(e)).index; fprintf(stdout,"cor vertice %d: %d\n", global_c,(*ORG(quad_dprev(e))).index);
		break;		
		case 'b': 
			divide_triangle2(e);
			(*ORG(quad_dprev(e))).index = 6 - (*ORG(e)).index - (*DST(e)).index; fprintf(stdout,"cor vertice %d: %d\n", global_c,(*ORG(quad_dprev(e))).index);
		break;
		case 'c': case 'n':
		break;
		case 'd':
			divide_triangle1(e);
		break;
		case 'e':
			f = quad_sym(quad_dprev(e));
			divide_triangle2(f);
		break;
		case 'E':
			f = quad_dnext(quad_onext(e));
			divide_triangle2(f);
		break;
		case 'f': case 'p':
			f = quad_sym(quad_dprev(e));
			divide_triangle2(f);
			f  = quad_dprev(f);
			divide_triangle1(f);
			f  = quad_sym(f);
			divide_triangle2(f);
		break;
		case 'F': case 'P':
			f = quad_dnext(quad_onext(e));
			divide_triangle2(f);
			f  = quad_dprev(f);
			divide_triangle1(f);
			f  = quad_sym(f);
			divide_triangle2(f);	
		break;
		case 'g': 
			divide_triangle1(e);
			divide_triangle2(e);
			f = quad_sym(quad_dprev(quad_dprev(e)));
			divide_triangle2(f);
		break;	
		case 'h':
			divide_triangle1(e);
			f = quad_sym(quad_onext(e));
			divide_triangle2(f);
			f = quad_dprev(f);
			divide_triangle2(f);
		break;
		case 'o':
			divide_triangle2(e);
		break;
		case 'q':
			divide_triangle2(e);
			f = quad_dprev(e);
			divide_triangle2(f);
			f = quad_dprev(f);
			divide_triangle2(f);
		break;
		case 'r':
			a = quad_make_edge();
			b = quad_dnext(quad_onext(e));
			quad_splice(b,a);
			c = quad_sym(quad_dprev(e));
			quad_splice(c,quad_sym(a));
			SET_ORG(a,ORG(b));
			SET_DST(a,DST(e));
			divide_triangle1(a);
			divide_triangle2(quad_sym(quad_onext(a)));
			divide_triangle2(quad_sym(a));
			divide_triangle2(quad_sym(quad_oprev(a)));
		break;
	}
}

delaunay_site_t *makesites(int nsites, bool_t normal)
  { delaunay_site_t *st = notnull(malloc(nsites*sizeof(delaunay_site_t)), "no mem");
    
    srandom(1000);

    int k;
    for (k=0; k < nsites; k++) 
      { st[k].index = -1;
        /* Generate a random point {p} in the square {[-1 _ =1]^2}: */
        r2_t p;
        if (normal) 
          { /* Central limit approx to normal distribuition: */
            double s = 0, t = 0;
            int j;
            for (j=0; j<NTOSS; j++)
              { s += 2 * (drandom() - 0.5);
                t += 2 * (drandom() - 0.5);
              }
            p.c[0] = s/sqrt(NTOSS);
            p.c[1] = t/sqrt(NTOSS);
          }
        else 
          { /* Uniform distribution: */
            p.c[0] = 2 * (drandom() - 0.5);
            p.c[1] = 2 * (drandom() - 0.5);
          }
        /* Convert {p} to integer homogeneous coordinates: */
        st[k].pt.c.c[0] = PLOTDIM;
        st[k].pt.c.c[1] = (int)(p.c[0]*PLOTDIM);
        st[k].pt.c.c[2] = (int)(p.c[1]*PLOTDIM);
	fprintf(stderr, "site %d: \t[\t%d\t%d\t%d\t]\n",k,st[k].pt.c.c[0],st[k].pt.c.c[1],st[k].pt.c.c[2]);
      }
    return st;
  }

void draw_colored_triangulation(PSStream *fps, quad_arc_t e)
  {
    auto void draw_colored_edge(quad_arc_t e);
    
    quad_arc_vec_t root = quad_arc_vec_make_desc(&e, 1);
    quad_enum(&root, draw_colored_edge);
    
    void draw_colored_edge(quad_arc_t e)
      { char text[2] = "0";
	delaunay_site_t *op = ORG(e);
        delaunay_site_t *dp = DST(e);
        r2_t rop = delaunay_r2_from_hi2(&(op->pt));
        r2_t rdp = delaunay_r2_from_hi2(&(dp->pt));
        pswr_segment(fps, rop.c[0], rop.c[1], rdp.c[0], rdp.c[1]);
	text[0] = (char)op->index + 48; 
	pswr_label(fps,text, rop.c[0], rop.c[1], 0, 0);
	text[0] = (char)dp->index + 48; 
	pswr_label(fps,text, rdp.c[0], rdp.c[1], 0, 0);
      }
 }

void plot_triangulation (quad_arc_t e, delaunay_site_t *st, int nsites, char *prefix, bool_t eps)
  { 
    double r = max_site_radius(st, nsites);
  
    /* Create Postscript document or EPS figure stream. */
    PSStream *fps = (PSStream *)open_ps_stream(r, prefix, eps);
    
    /* Start a new picture: */
    double wm = 0.2*r;
    double xmin = -wm/2; double xmax = +wm/2;
    double ymin = -wm/2; double ymax = +wm/2;
    pswr_new_picture(fps, xmin,xmax, ymin, ymax);
    
    /* Plot Delaunay edges, solid: */
    float penwd = (eps ? 0.20 : 0.10);
    pswr_set_pen(fps, 0,0,0, penwd, 0.0, 0.0);
    draw_colored_triangulation(fps, e);

    if (! eps)
      { /* Add caption and frame: */
        pswr_set_pen(fps, 0,0,0, 0.10, 0.0, 0.0);
        pswr_add_caption(fps, "Voronoi/Delaunay diagram", 0.5);
        pswr_set_pen(fps, 0,0,0, 0.20, 0.0, 0.0);
        pswr_frame(fps);
      }
    /* We are done: */
    pswr_close_stream(fps);
  }

void plot_triangulation_debug(quad_arc_t e, delaunay_site_t *st, int nsites, char *prefix, bool_t eps, quad_arc_list list)
   {
    double r = max_site_radius(st, nsites);
  
    /* Create Postscript document or EPS figure stream. */
    PSStream *fps = (PSStream *)open_ps_stream(r, prefix, eps);
    
    /* Start a new picture: */
    double wm = 0.2*r;
    double xmin = -wm/2; double xmax = +wm/2;
    double ymin = -wm/2; double ymax = +wm/2;
    pswr_new_picture(fps, xmin,xmax, ymin, ymax);
    
    /* Plot Delaunay edges, solid: */
    float penwd = (eps ? 0.20 : 0.10);
    pswr_set_pen(fps, 0,0,0, penwd, 0.0, 0.0);
    draw_delaunay_edges(fps, e);

    /*draw edges too debug*/
    char text[2] = "a";
    int i; double x,y;
    quad_arc_list l;
    x =  (((double)(*DST(e)).pt.c.c[1]) / ((double)(*DST(e)).pt.c.c[0]) + ((double)(*DST(quad_onext(e))).pt.c.c[1]) / ((double)(*DST(quad_onext(e))).pt.c.c[0]) +
	((double)(*ORG(e)).pt.c.c[1]) / ((double)(*ORG(e)).pt.c.c[0]))/3;
    y =  (((double)(*DST(e)).pt.c.c[2]) / ((double)(*DST(e)).pt.c.c[0]) + ((double)(*DST(quad_onext(e))).pt.c.c[2]) / ((double)(*DST(quad_onext(e))).pt.c.c[0]) +
	((double)(*ORG(e)).pt.c.c[2]) / ((double)(*ORG(e)).pt.c.c[0]))/3;
    pswr_label(fps,"-", x, y, 0, 0);
    for (l=list;l->next!=list;l=l->next){
	x =  (((double)(*DST(l->e)).pt.c.c[1]) / ((double)(*DST(l->e)).pt.c.c[0]) + 
		((double)(*DST(quad_onext(l->e))).pt.c.c[1]) / ((double)(*DST(quad_onext(l->e))).pt.c.c[0]) +
		((double)(*ORG(l->e)).pt.c.c[1]) / ((double)(*ORG(l->e)).pt.c.c[0]))/3;
	y =  (((double)(*DST(l->e)).pt.c.c[2]) / ((double)(*DST(l->e)).pt.c.c[0]) + 
		((double)(*DST(quad_onext(l->e))).pt.c.c[2]) / ((double)(*DST(quad_onext(l->e))).pt.c.c[0]) +
		((double)(*ORG(l->e)).pt.c.c[2]) / ((double)(*ORG(l->e)).pt.c.c[0]))/3;
	pswr_label(fps,text, x, y, 0, 0);
	text[0]++; if(text[0]=='z'+1) text[0]='A';
    }
	x =  (((double)(*DST(l->e)).pt.c.c[1]) / ((double)(*DST(l->e)).pt.c.c[0]) + 
		((double)(*DST(quad_onext(l->e))).pt.c.c[1]) / ((double)(*DST(quad_onext(l->e))).pt.c.c[0]) +
		((double)(*ORG(l->e)).pt.c.c[1]) / ((double)(*ORG(l->e)).pt.c.c[0]))/3;
	y =  (((double)(*DST(l->e)).pt.c.c[2]) / ((double)(*DST(l->e)).pt.c.c[0]) + 
		((double)(*DST(quad_onext(l->e))).pt.c.c[2]) / ((double)(*DST(quad_onext(l->e))).pt.c.c[0]) +
		((double)(*ORG(l->e)).pt.c.c[2]) / ((double)(*ORG(l->e)).pt.c.c[0]))/3;
    pswr_label(fps,text, x, y, 0, 0);

    if (! eps)
      { /* Add caption and frame: */
        pswr_set_pen(fps, 0,0,0, 0.10, 0.0, 0.0);
        pswr_add_caption(fps, "Voronoi/Delaunay diagram", 0.5);
        pswr_set_pen(fps, 0,0,0, 0.20, 0.0, 0.0);
        pswr_frame(fps);
      }
    /* We are done: */
    pswr_close_stream(fps);
   }