#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);

void execute_division(quad_arc_t e, char type);

quad_arc_list make_edge_list(quad_arc_t e, int nsites);

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);
    list = make_edge_list(e, nsites);
    plot_triangulation_debug(e, st, nsites, "out/debug", eps, list);
    fprintf(stderr, "Plotting delaunay...\n");
    plot_triangulation(e, st, nsites, "out/delrandom", eps);


//  TESTE - MOD6
    e = quad_onext(quad_dprev(e));
    execute_division(e, 'd');
    e = quad_sym(quad_dnext(e));
    execute_division(e, 'd');
    e = quad_sym(quad_dnext(e));
    execute_division(e, 'd');
    execute_division(quad_sym(e), 'r');
//

    plot_triangulation(e, st, nsites, "out/delrandom2", 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; 
	}
}

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 (g == quad_lnext(quad_lnext(quad_lnext(g))))
	add_quad_arc_l(&front_list_begin, g);
	if (h == quad_lnext(quad_lnext(quad_lnext(h))))
	add_quad_arc_l(&front_list_begin, h);
//fprintf(stderr, "%d %d %d\n",(*ORG(front_list_begin->e)).pt.c.c[1], (*ORG(front_list_begin->next->e)).pt.c.c[1], (*ORG(front_list_begin->next->next->e)).pt.c.c[1]);
	while(front_list_begin!=NULL){
		front_list = front_list_begin;
		do{
			f = front_list->e; 
fprintf(stderr, "\naresta do front a processar: %d %d\t%d %d\n", (*ORG(f)).pt.c.c[1], (*ORG(f)).pt.c.c[2], (*DST(f)).pt.c.c[1], (*DST(f)).pt.c.c[2]);
			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);				fprintf(stderr, "caso1\n");
				add_quad_arc_l(&edge_list, quad_sym(g));
fprintf(stderr, "adiciona na lista de arestas: %d %d\t%d %d\n", (*ORG(quad_sym(g))).pt.c.c[1], (*ORG(quad_sym(g))).pt.c.c[2], (*DST(quad_sym(g))).pt.c.c[1], (*DST(quad_sym(g))).pt.c.c[2]);
				remove_quad_arc_l(&front_list_begin,&front_list);
fprintf(stderr, "remove do front: %d %d\t%d %d\n", (*ORG(right->e)).pt.c.c[1], (*ORG(right->e)).pt.c.c[2], (*DST(right->e)).pt.c.c[1], (*DST(right->e)).pt.c.c[2]);
				remove_quad_arc_l(&front_list_begin,&right);
				if (g == quad_lnext(quad_lnext(quad_lnext(g)))){ //se g nao for aresta de borda
fprintf(stderr, "adiciona no front: %d %d\t%d %d\n", (*ORG(g)).pt.c.c[1], (*ORG(g)).pt.c.c[2], (*DST(g)).pt.c.c[1], (*DST(g)).pt.c.c[2]);
					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);				fprintf(stderr, "caso2\n");				
				add_quad_arc_l(&edge_list, quad_sym(g));
fprintf(stderr, "adiciona na lista de arestas: %d %d\t%d %d\n", (*ORG(quad_sym(g))).pt.c.c[1], (*ORG(quad_sym(g))).pt.c.c[2], (*DST(quad_sym(g))).pt.c.c[1], (*DST(quad_sym(g))).pt.c.c[2]);
				remove_quad_arc_l(&front_list_begin,&front_list);
fprintf(stderr, "remove do front: %d %d\t%d %d\n", (*ORG(left->e)).pt.c.c[1], (*ORG(left->e)).pt.c.c[2], (*DST(left->e)).pt.c.c[1], (*DST(left->e)).pt.c.c[2]);
				remove_quad_arc_l(&front_list_begin,&left);
				if (g == quad_lnext(quad_lnext(quad_lnext(g)))){ //se g nao for aresta  de borda
fprintf(stderr, "adiciona no front: %d %d\t%d %d\n", (*ORG(g)).pt.c.c[1], (*ORG(g)).pt.c.c[2], (*DST(g)).pt.c.c[1], (*DST(g)).pt.c.c[2]);
					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
//fprintf(stderr, "%d %d\t%d %d\n", (*ORG(f)).pt.c.c[1], (*ORG(f)).pt.c.c[2], (*DST(f)).pt.c.c[1], (*DST(f)).pt.c.c[2]);
				(*ORG(quad_dprev(f))).index = 0;
				g = quad_onext(f);				fprintf(stderr, "caso3\n");
				h = quad_dprev(f);
				add_quad_arc_l(&edge_list, f);
fprintf(stderr, "adiciona na lista de arestas: %d %d\t%d %d\n", (*ORG(f)).pt.c.c[1], (*ORG(f)).pt.c.c[2], (*DST(f)).pt.c.c[1], (*DST(f)).pt.c.c[2]);
				remove_quad_arc_l(&front_list_begin,&front_list);
				if (g == quad_lnext(quad_lnext(quad_lnext(g)))){ //se g nao for aresta de borda
fprintf(stderr, "adiciona no front: %d %d\t%d %d\n", (*ORG(g)).pt.c.c[1], (*ORG(g)).pt.c.c[2], (*DST(g)).pt.c.c[1], (*DST(g)).pt.c.c[2]);
					add_quad_arc_l(&front_list_begin, g);}	
				if (h == quad_lnext(quad_lnext(quad_lnext(h)))){ //se h nao for aresta  de borda
fprintf(stderr, "adiciona no front: %d %d\t%d %d\n", (*ORG(h)).pt.c.c[1], (*ORG(h)).pt.c.c[2], (*DST(h)).pt.c.c[1], (*DST(h)).pt.c.c[2]);
					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;				fprintf(stderr, "caso4\n");
			}
		}while (front_list!=front_list_begin);
        }
	return edge_list;
}

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));
}

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

	switch(type){
		case 'a': case 'c': case 'n':
		break;
		case 'b': case 'o':
			divide_triangle2(e);
		break;
		case 'd':
			divide_triangle1(e);
		break;
		case 'e':
			f = quad_sym(quad_dprev(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 'g':
			divide_triangle1(e);
			f = quad_sym(quad_onext(e));
			divide_triangle2(f);
			f = quad_dprev(f);
			divide_triangle2(f);
		break;
		case 'h': 
			divide_triangle1(e);
			divide_triangle2(e);
			f = quad_sym(quad_dprev(quad_dprev(e)));
			divide_triangle2(f);
		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 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_delaunay_edges(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);
   }