#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 MC (delaunay_MAX_COORD)

void divide_triangle1(quad_arc_t e);

void divide_triangle2(quad_arc_t e);

void execute_division(quad_arc_t e, char type);

delaunay_site_t *makesites(int nsites, bool_t normal);

int main(int argc, char **argv)
  { 
    int nsites = atoi(argv[1]);
    bool_t normal = 1;
    bool_t eps = 1;
    quad_arc_t e;
    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, "Plotting delaunay...\n");
    plot_delaunay(e, st, nsites, "out/delrandom", eps);

/*  MOD1
    e = quad_dprev(e);
    execute_division(e, 'd');
    e = quad_sym(e);
    execute_division(e, 'b');
*/

/*  MOD2
    quad_arc_t f = quad_dprev(e);
    execute_division(f, 'd');
    quad_arc_t g = quad_onext(quad_dnext(f));
    execute_division(g, 'd');
    execute_division(e, 'g');
*/

/*  MOD3
    quad_arc_t f = quad_dprev(e);
    execute_division(f, 'd');
    quad_arc_t g = quad_onext(quad_dnext(f));
    execute_division(g, 'd');
    execute_division(e, 'h');
*/

    e = quad_dprev(e);
    execute_division(e, 'd');
    e = quad_dnext(e);
    execute_division(e, 'f');




    plot_delaunay(e, st, nsites, "out/delrandom2", eps);
    return(0);
  }

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 = notnull(malloc(sizeof(delaunay_site_t)), "no mem");

	e_odata->pt.c.c[0] = (*ORG(e)).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] + (*DST(e)).pt.c.c[1] * (*ORG(e)).pt.c.c[0];   
	e_odata->pt.c.c[2] = (*ORG(e)).pt.c.c[2] * (*DST(e)).pt.c.c[0] + (*DST(e)).pt.c.c[2] * (*ORG(e)).pt.c.c[0];
	e_odata->index = -1; 

	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;

	switch(type){
		case 'a':

		break;
		case 'b':

			divide_triangle2(e);
		break;
		case 'c':

		break;
		case 'd':

			divide_triangle1(e);
		break;
		case 'e':

		break;
		case 'f':
			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;
	}
}

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 = k;
        /* 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 = delaunay_hi2_from_r2(&p);
      }
    if (VERBOSE) 
      { for (k = 0; k < nsites; k++)
          { assert(st[k].index == k);
            deldebug_print_site ("given", &(st[k]));  fprintf(stderr, "\n");
            assert(abs(st[k].pt.c.c[0]) <= MC);
            assert(abs(st[k].pt.c.c[1]) <= MC);
            assert(abs(st[k].pt.c.c[2]) <= MC);
          }
      }
    return st;
  }