#include <stdio.h>
#include <stdlib.h>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <list>
#include <iterator>
#include <queue>
#include <deque>
#include <stack>
#include <algorithm>
#include <set>
#include <math.h>
#include <assert.h>
#include <limits>
#include <climits>
#include <string>
#include <time.h>
#include <malloc.h>
#include <getopt.h>
#include "tri.h"
#include "vor.h"
#include "background.h"
#include "quad.h"

#define tprintf(...) printf("### %lf ### ", clock ()*1e-6); printf(__VA_ARGS__)

// void init_inside(Triangulation *tri)
// {
//     int i0, i1;
//     int rs = tri->getVertexCount()-4;
//     Vertex *pv, *pn;
// 
//     Triangle *t, *s;
//     int j;
//     int ret;
// 
//     printf("*** Init START ***\n");
// 
//     for (t = tri->getTri(); t < tri->getTri()+ tri->getTriangleCount(); ++t)
//         t->set_flag(-1);
// 
//     for (int i = 0; i < rs; i++)
//     {
//         i0 = i+4;
//         i1 = (i+1 < rs) ? i0+1 : (i+1)%rs + 4;
//         pv = tri->get_vertex(i0);
//         pn = tri->get_vertex(i1);
// 
//         // flag de triangles
//         ret = tri->get_tri(pv, pn, t, j=0);
//         assert(ret != 0);
// 
//         s = t->get_adjacent(j);
//         t->set_flag(1);
//         s->set_flag(0);
//     }
// 
//     int cont = 1;
//     while (cont)
//     {
// 	cont = 0;
//         for (t = tri->getTri(); t < tri->getTri()+ tri->getTriangleCount(); ++t)
//             for (int i = 0; i < 3; i++)
//                 if (t->get_adjacent(i) != NULL && t->get_adjacent(i)->get_flag() == -1)
//                 {
//                     t->get_adjacent(i)->set_flag(t->get_flag());
//                     cont = 1;
//                 }
//     }

//     printf("*** Init END ***\n");
// }

void flag_border(Triangulation *tri)
{
    Triangle *t;
    for (t = tri->getTri(); t < tri->getTri()+ tri->getTriangleCount(); ++t)
    {
        tri->get_vertex(t->get_vertex(0))->flag = 1;
        tri->get_vertex(t->get_vertex(1))->flag = 1;
        tri->get_vertex(t->get_vertex(2))->flag = 1;
        t->set_flag(1);
    }
    for (Vertex* v = tri->getVertices(); v < tri->getVertices()+ tri->getVertexCount(); ++v)
    {
        if ((v->p[0] < 0.0 || v->p[0] > 10.0) || (v->p[1] < 0.0 || v->p[1] > 10.0))
        {
            v->flag = -1;
        }
    }
}

void insert_center(Triangulation *tri)
{
    for (int loop = 0; loop < 2500; loop++)
    {
        double max_vol = 0.0;
        Triangle *max_t = NULL;
        for (Triangle *t = tri->getTri(); t < tri->getTri()+ tri->getTriangleCount(); ++t)
        {
            if (t != NULL && t->get_flag() == 1)
            {
                double v = tri->volume(t);
// 		double v = tri->get_center(t)[2];
                if (v > max_vol)
                {
                    max_vol = v;
                    max_t = t;
                }
            }
        }
        assert(max_t != NULL);
        double max_vol_a = 0.0;
        Triangle *max_t_a = NULL;
        for (int i = 0; i < 3; i++)
        {
            Triangle *at = max_t->get_adjacent(i);
            if (at != NULL && at->get_flag() == 1)
            {
                double va = tri->volume(at);
// 		double va = tri->get_center(at)[2];
                if (va > max_vol_a)
                {
                    max_vol_a = va;
                    max_t_a = at;
                }
            }
        }
        const double *ct, *cta;
        double vx[2];
        ct = tri->get_center(max_t);
        if (max_t_a != NULL)
        {
            cta = tri->get_center(max_t_a);
            vx[0] = (ct[0]+cta[0])/2.0;
            vx[1] = (ct[1]+cta[1])/2.0;
        }
        else
        {
            vx[0] = ct[0];
            vx[1] = ct[1];
        }
        Vertex *pv = tri->add_vertex(vx);
        tri->insert_lawson(pv);
    }
}

void classic_tri(double *vx, int nb_vx)
{
    Triangulation *tri = new Triangulation;
    Background *bg = new Background(0.0, 0.0, 10.0, 10.0);
    tri->set_background(bg);
    
    tprintf("*** Triangulation START ***\n");
    tri->init(vx, nb_vx, 1.0, false);
    tprintf("*** Triangulation END ***\n");

    //     init_inside(&m_trg);
//     insert_center(&m_trg);
    
//     assert(m_trg.check_tri() == 0);
    tri->export_vtk();
    VoronoiDiagram m_vor(tri, bg);
    m_vor.export_vtk();
//     m_vor.print();
}

void lloyd(double *vx, int nb_vx, double relax, int max_loop)
{
    tprintf("*** Lloyd START ***\n");
    
    int dump_rate = 10;
    
//     Background *bg = new Background("portrait.png", 2, 0.0, 0.0, 10.0, 10.0);
//     Background *bg = new Background("vincent.png", 0, 0.0, 0.0, 10.0, 10.0);
    Background *bg = new Background(0.0, 0.0, 10.0, 10.0);
    
    progress_bar pg;    
    pg.start(max_loop);   
        
    for (int i = 0; i < max_loop; i++)
    {
	pg.progress(i);
        Triangulation *tri = new Triangulation;
        tri->init(vx, nb_vx, 1.0, false);
// 	flag_border(tri);
//      tri->export_vtk();

        VoronoiDiagram *vor = new VoronoiDiagram(tri, bg);
	flag_border(tri);
	vor->compute_centroid();
        if (i%dump_rate == 0)
            vor->export_vtk();
	for (int j = 0; j < nb_vx; j++)
	{
            VoronoiCell *vc = vor->get_cell(j);
#if defined (POWER)
                int index = 3*j;
#else
                int index = 2*j;
#endif
            if (vc->_center->flag == 1)
            {
                vx[index] = vc->_center->p[0] + (vc->_centroid->p[0]-vc->_center->p[0])*relax;
                vx[index+1] = vc->_center->p[1] + (vc->_centroid->p[1]-vc->_center->p[1])*relax;
            }
            else
            {
                vx[index] = vc->_center->p[0];
                vx[index+1] = vc->_center->p[1];
            }
            // calcul de la ponderation
#if defined (POWER)
	double size = bg->get_size(&vx[index]);
	assert(size < 1.0);
	vx[index+2] = size/12.;
#endif
	}
	if (i == max_loop-1)
	{
	  pg.end();  
#if defined (QUAD)
	  QuadMesh *qmesh = new QuadMesh(tri);
	  qmesh->export_vtk();
#endif
	  tri->export_vtk();
	}
	delete tri;
	delete vor;
    }
    
    
    tprintf("*** Lloyd END ***\n");
}

void test_orient(double *vx, int nb_vx)
{
    tprintf("*** Triangulation START ***\n");
    
    Background *bg = new Background;
    
    double *vx_save = new double[2 * nb_vx];
    double angle = 0.0;
    for (int i = 0; i < 360; i++)
    {	
        for (int i = 0; i < 2*nb_vx; i++)
            vx_save[i] = vx[i];
	
        Triangulation *tri = new Triangulation;
	tri->set_angle(angle);
        tri->init(vx_save, nb_vx, 1.0, false);
        VoronoiDiagram m_vor(tri, bg);
	m_vor.compute_centroid();
        m_vor.export_vtk();
	delete tri;
	angle += M_PI/180.0;
    }
    tprintf("*** Triangulation END ***\n");
}

int main(int argc, char** argv)
{
      int nb = 50000;
      srand(0);
      
//     srand(0);
//     for (int i = 0; i < 100000; i++)
//     {
//       printf("%d %d\n", rand()%1000, rand()%1000);
//     }
//     exit(0);
//     
      
//     printf("%d\n", nb);
//     for (int i = 0; i < nb; i++)
//     {
//       printf("%d %lf %lf %lf\n", i, (double)rand()/(double)RAND_MAX, (double)rand()/(double)RAND_MAX, 0.0);
//     }
//     exit(0);
    
//         srand(0);
//     printf("%d\n", nb);
//     for (int i = 0; i < nb; i++)
//     {
//       printf("%lf %lf %lf\n", 10.0*(double)rand()/(double)RAND_MAX, 10.0*(double)rand()/(double)RAND_MAX, 0.1*(double)rand()/(double)RAND_MAX);
//     }
//     exit(0);
// 
//     int num = 100;
//     for (int i = 0; i < num; i++)
//       for (int j = 0; j < num; j++)
//     {
//       printf("%d %d\n", i, j);
//     }
//     exit(0);

//     Triangulation m_trg;
    char filename[100];
    int error = 0;
    sprintf(filename, "%s.dat", argv[1]);
    printf("Reading %s ....\n", filename);
    FILE *fp = fopen(filename, "r");
    int nb_vx;
    if (fscanf(fp, "%d", &nb_vx) != 1)
        error++;
    printf("--> %d vertices\n", nb_vx);
#if defined (POWER)  
        double *vx = new double[3 * nb_vx];
// 	    srand(0);
    for (int i = 0; i < nb_vx; i++)
    {
        double x, y, z, w;
	int index;
// 	if (fscanf(fp, "%lf %lf", &x, &y) != 2)
     if (fscanf(fp, "%lf %lf %lf", &x, &y, &w) != 3)
// 	  if (fscanf(fp, "%lf %lf %lf", &x, &y, &z) != 3)
// 	if (fscanf(fp, "%d %lf %lf %lf", &index, &x, &y, &z) != 4)
            error++; 
        vx[3*i] = x;
        vx[3*i+1] = y;
// 	vx[3*i+2] = w/2000.;
	vx[3*i+2] = (x)/8.;
// 	vx[3*i+2] = 0.0;
    }
#else
    double *vx = new double[2 * nb_vx];
    for (int i = 0; i < nb_vx; i++)
    {
        double x, y, z;
        int index;
//         if (fscanf(fp, "%lf %lf", &x, &y) != 2)
        if (fscanf(fp, "%lf %lf %lf", &x, &y, &z) != 3)
//         if (fscanf(fp, "%d %lf %lf %lf", &index, &x, &y, &z) != 4)
                error++;
        vx[2*i] = (double)x;
        vx[2*i+1] = (double)y;
    }
#endif
        
    fclose(fp);

#if defined(CLASSIC)       
     classic_tri(vx, nb_vx);
#endif     
//     test_orient(vx, nb_vx); 
#if defined(LLOYD)
    lloyd(vx, nb_vx, 1.0, 500);
#endif    
    
    

}