#include <cstring>
#include <cstdio>
#include "tri.h"
#include "background.h"
#include <stdlib.h>
#include <queue>



// Algorithm.
// The initial triangulation is the two triangles that triagulate
// the bounding box.  Then we:
// First strategy: insert each site of |s| one by one randomly.
// Second strategy: insert each site in an enumeration order of buckets such that two consecutive buckets are neighbors.

/**
 *    The triangulation constructor.
 *    @param *s the vertices array (s[0] = vertex[0].x & s[1] = vertex[0].y).
 *    @param n number of vertices.
 *    @param epsilon how much enlarge the bounding box.
 *    @param regular triangulation is or not regular.
 */
Triangulation::Triangulation ():
        nv (0),
        ntri (0),
        last_tri (0),
        leftoncount (0)
{
    cav = new Cavity(this);
    ifree_t = 0;
}


Triangulation::~Triangulation()
{
    delete []v;
    delete []tri;
    delete []free_t;
    delete cav;
}

int dump_steps = 0;

void dump_conf(Triangulation *tri, Triangle * t0, Triangle * t1, Triangle *t2)
{
    static int count = 0;
    char filename[100];
    sprintf(filename, "dump_conf_%d.vtk", count++);
    printf(" write file --> %s\n", filename);

    FILE *fp_rw = fopen (filename, "w");

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET POLYDATA\n");
    int num = 6; 
    if (t0->get_adjacent(0) == NULL)
        num--;
    if (t1->get_adjacent(1) == NULL)
        num--;
    if (t2->get_adjacent(2) == NULL)
        num--;
    fprintf (fp_rw, "POINTS %d float\n", 3*num);

    Vertex *pv;

    for (int i = 0; i < 3; i++)
    {
        pv = tri->get_vertex(t0->get_vertex(i));
        fprintf(fp_rw, "%lf %lf %lf\n", pv->p[0], pv->p[1], 0.0);
    }
    
    if (t0->get_adjacent(0) != NULL)
    for (int i = 0; i < 3; i++)
    {
        pv = tri->get_vertex(t0->get_adjacent(0)->get_vertex(i));
        fprintf(fp_rw, "%lf %lf %lf\n", pv->p[0], pv->p[1], 0.0);
    }
    
    for (int i = 0; i < 3; i++)
    {
        pv = tri->get_vertex(t1->get_vertex(i));
        fprintf(fp_rw, "%lf %lf %lf\n", pv->p[0], pv->p[1], 0.0);
    }
    if (t1->get_adjacent(1) != NULL)
        for (int i = 0; i < 3; i++)
    {
        pv = tri->get_vertex(t1->get_adjacent(1)->get_vertex(i));
        fprintf(fp_rw, "%lf %lf %lf\n", pv->p[0], pv->p[1], 0.0);
    }
    
    for (int i = 0; i < 3; i++)
    {
        pv = tri->get_vertex(t2->get_vertex(i));
        fprintf(fp_rw, "%lf %lf %lf\n", pv->p[0], pv->p[1], 0.0);
    }
    if (t2->get_adjacent(2) != NULL)
        for (int i = 0; i < 3; i++)
    {
        pv = tri->get_vertex(t2->get_adjacent(2)->get_vertex(i));
        fprintf(fp_rw, "%lf %lf %lf\n", pv->p[0], pv->p[1], 0.0);
    }

        fprintf (fp_rw, "POLYGONS %d %d\n", num, 4*num);
        fprintf (fp_rw, "%d %d %d %d\n", 3, 0, 1, 2);
        fprintf (fp_rw, "%d %d %d %d\n", 3, 3, 4, 5);
	fprintf (fp_rw, "%d %d %d %d\n", 3, 6, 7, 8);
    int index = 9;
        if (t0->get_adjacent(0) != NULL)
        fprintf (fp_rw, "%d %d %d %d\n", 3, index++, index++, index++);
	    if (t1->get_adjacent(1) != NULL)
        fprintf (fp_rw, "%d %d %d %d\n", 3, index++, index++, index++);
    if (t2->get_adjacent(2) != NULL)
        fprintf (fp_rw, "%d %d %d %d\n", 3, index++, index++, index++);


    fclose(fp_rw);
}

void Triangulation::debug_step(Triangle *t, Triangle *s)
{
    if (dump_steps)
    {
        t->pos2 = 1;
        if (s)
            s->pos2 = 2;
	check_delaunay();
        export_vtk();
        t->pos2 = 0;
        if (s)
            s->pos2 = 0;
    }
}

void Triangulation::check_power(double *vx, int size)
{
    ANNpointArray nodes = annAllocPts(size, 3);

    for (int i = 0; i < size; i++)
    {
        nodes[i][0] = vx[3*i];
        nodes[i][1] = vx[3*i+1];
        nodes[i][2] = vx[3*i+2]*vx[3*i+2];
    }
    ANNkd_tree *kdtree = new ANNkd_tree(nodes, size, 2);

    srand(0);
    ANNidxArray index = new ANNidx;
    ANNdistArray dist = new ANNdist;
    ANNpoint xyz = new ANNcoord[2];
    
    for (int l = 0; l < 10*size; l++)      
    {
	int i = rand()%(size-1);
//     for (int i = 0; i < size; i++)      
//     {
// 	int i = rand()%(size-1);
        xyz[0] = vx[3*i];
        xyz[1] = vx[3*i+1];
        double w0 = vx[3*i+2]*vx[3*i+2];
        kdtree->annkSearch(xyz, 2, index, dist);
        ANNpointArray pts = kdtree->thePoints();
        ANNpoint cl = pts[index[1]];
        double w1 = cl[2];
        double d = (xyz[0]-cl[0])*(xyz[0]-cl[0])+(xyz[1]-cl[1])*(xyz[1]-cl[1]);
// 	printf("d = %lf, w0 = %lf, w1 = %lf\n", d, w0, w1);
        if (d <= fabs (w0-w1))
        {
// 	  printf("Problem with weight p %d -- p %d\n", i, index[1]);
//             double tmin = (d-w1)/(2.0*d);
// // 	  printf("tmin = %lf\n", tmin);
//             if (tmin < 0.0)
//                 tmin = 0.0;
//             double t = tmin + (1.0-tmin)*(double)rand()/(double)RAND_MAX; // 0 < t < 1
// // 	  printf("t = %lf\n", t);
//             w0 = 2.0*t*d-d+w1;
// // 	  printf("w0 = %lf\n", w0);
//             assert(w0 >= 0.0);
// 	    w0 = w1;	    
//             vx[3*i+2] = sqrt(w0);
// 	    vx[3*i+2] = 0.0;
        }
    }
    delete kdtree;
    delete index;
    delete dist;
    delete xyz;
}

void Triangulation::export_vtk()
{
    static int count = 0;
    char filename[100];
    int nb_tri;
    
    printf("*** Export mesh START ***\n");

    sprintf(filename, "test_%d.vtk", count);
    printf(" write file --> %s\n", filename);

    FILE *fp_rw = fopen (filename, "w");

    printf(" %d vertices, %d triangles\n", getVertexCount(), getTriangleCount());

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET POLYDATA\n");
    fprintf (fp_rw, "POINTS %d float\n", getVertexCount());

    Vertex *pv;

    for (int i = 0; i < getVertexCount(); i++)
    {
        pv = get_vertex(i);
	double p[2] = {pv->p[0], pv->p[1]};
#if defined(ALIGN)
	disalign(p);
#endif
// #if defined(ANGLE)
// 	
//         double x = p[0];
//         double y = p[1];
//         p[0] = x*cos(_angle) - y*sin(_angle);
//         p[1] = x*sin(_angle) + y*cos(_angle);
// #endif
        fprintf(fp_rw, "%lf %lf %lf\n", p[0], p[1], 0.0);
    }
    
    Triangle *tri = getTri();
    Triangle *tend = tri + getTriangleCount();
    Triangle *t;
    
    int tri_count = 0;
    for (t = tri; t < tend; ++t)
      if (t->get_flag() != FREED_TRIANGLE_FLAG)
	tri_count++;
      else
	dprintf("tetra %d is erased !\n", (int)(t-getTri()));

    fprintf (fp_rw, "POLYGONS %d %d\n", tri_count, 4*tri_count);
    for (t = tri; t < tend; ++t)
      if (t->get_flag() != FREED_TRIANGLE_FLAG)
        fprintf (fp_rw, "%d %d %d %d\n", 3, get_vertex (t, 0), get_vertex (t, 1), get_vertex (t, 2));

    fprintf (fp_rw, "CELL_DATA %d\n", tri_count);

//     for (t = tri; t < tend; ++t)
//         if ((get_center(t))[0] != INFINITY && (get_center(t))[1] != INFINITY)
//             fprintf (fp_rw, "%d\n", 1);
//         else
//             fprintf (fp_rw, "%d\n", 0);
//         for (t = tri; t < tend; ++t)
// 	  if (t->get_flag() != FREED_TRIANGLE_FLAG)
//             fprintf (fp_rw, "%d\n", t->get_flag());
	  
    fprintf (fp_rw, "SCALARS index int 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");
    int i = 0;
    for (t = tri; t < tend; ++t)
        if (t->get_flag() != FREED_TRIANGLE_FLAG)
            fprintf (fp_rw, "%d\n", i++);
	
    fprintf (fp_rw, "SCALARS flag int 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");
    for (t = tri; t < tend; ++t)
        if (t->get_flag() != FREED_TRIANGLE_FLAG)
            fprintf (fp_rw, "%d\n", t->get_flag());
    
    fprintf (fp_rw, "SCALARS angle float 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");
    for (t = tri; t < tend; ++t)
        if (t->get_flag() != FREED_TRIANGLE_FLAG)
            fprintf (fp_rw, "%f\n", t->get_angle());

    fclose(fp_rw);
   
#if defined(Linf) && defined(POWER)
    sprintf(filename, "test_power_%d.vtk", count);
    printf(" write file --> %s\n", filename);

    fp_rw = fopen (filename, "w");

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET POLYDATA\n");
    fprintf (fp_rw, "POINTS %d float\n", 4*getVertexCount());
    
    for (int j = 0; j < getVertexCount(); j++)
    {
        pv = get_vertex(j);
        double p0[2] = {pv->p[0]-pv->w, pv->p[1]-pv->w};
        double p1[2] = {pv->p[0]+pv->w, pv->p[1]-pv->w};
        double p2[2] = {pv->p[0]+pv->w, pv->p[1]+pv->w};
        double p3[2] = {pv->p[0]-pv->w, pv->p[1]+pv->w};
	fprintf(fp_rw, "%lf %lf %lf\n", p0[0], p0[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p1[0], p1[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p2[0], p2[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p3[0], p3[1], 0.0);
    }
    
    fprintf (fp_rw, "POLYGONS %d %d\n", getVertexCount(), 5*getVertexCount());
   
    for (int j = 0; j < getVertexCount(); j++)
    {
        fprintf (fp_rw, "%d %d %d %d %d\n", 4, 4*j, 4*j+1, 4*j+2, 4*j+3);
    }
    
    fprintf (fp_rw, "CELL_DATA %d\n", getVertexCount());
    fprintf (fp_rw, "SCALARS index int 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");      
    for (int j = 0; j < getVertexCount(); j++)
      fprintf (fp_rw, "%d\n", j);
    
    fclose(fp_rw);    
#endif  
    
#if defined(CIRCLE_DUMP)   
    
#if defined(Linf)

    sprintf(filename, "test_polygon_%d.vtk", count);
    printf(" write file --> %s\n", filename);

    fp_rw = fopen (filename, "w");

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET POLYDATA\n");
    fprintf (fp_rw, "POINTS %d float\n", 4*getTriangleCount());

    for (t = tri; t < tend; ++t)
    {
        
#if defined(POWER)
        double p0[2] = {(get_center(t))[0]-((get_center(t))[2]+get_vertex(t->get_vertex(0))->w), (get_center(t))[1]-((get_center(t))[3]+get_vertex(t->get_vertex(0))->w)};
        double p1[2] = {(get_center(t))[0]+((get_center(t))[2]+get_vertex(t->get_vertex(0))->w), (get_center(t))[1]-((get_center(t))[3]+get_vertex(t->get_vertex(0))->w)};
        double p2[2] = {(get_center(t))[0]+((get_center(t))[2]+get_vertex(t->get_vertex(0))->w), (get_center(t))[1]+((get_center(t))[3]+get_vertex(t->get_vertex(0))->w)};
        double p3[2] = {(get_center(t))[0]-((get_center(t))[2]+get_vertex(t->get_vertex(0))->w), (get_center(t))[1]+((get_center(t))[3]+get_vertex(t->get_vertex(0))->w)};
#else
        double p0[2] = {(get_center(t))[0]-(get_center(t))[2], (get_center(t))[1]-(get_center(t))[3]};
        double p1[2] = {(get_center(t))[0]+(get_center(t))[2], (get_center(t))[1]-(get_center(t))[3]};
        double p2[2] = {(get_center(t))[0]+(get_center(t))[2], (get_center(t))[1]+(get_center(t))[3]};
        double p3[2] = {(get_center(t))[0]-(get_center(t))[2], (get_center(t))[1]+(get_center(t))[3]};
#endif
#if defined(ANGLE)
	double angle = t->get_angle();
	
	
	
        double x, y;
	x = p0[0];
	y = p0[1];
        p0[0] = x*cos(angle) - y*sin(angle);
        p0[1] = x*sin(angle) + y*cos(angle);
        x = p1[0];
        y = p1[1];
        p1[0] = x*cos(angle) - y*sin(angle);
        p1[1] = x*sin(angle) + y*cos(angle);
        x = p2[0];
        y = p2[1];
        p2[0] = x*cos(angle) - y*sin(angle);
        p2[1] = x*sin(angle) + y*cos(angle);
        x = p3[0];
        y = p3[1];
        p3[0] = x*cos(angle) - y*sin(angle);
        p3[1] = x*sin(angle) + y*cos(angle);
#endif

        fprintf(fp_rw, "%lf %lf %lf\n", p0[0], p0[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p1[0], p1[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p2[0], p2[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p3[0], p3[1], 0.0);
    }

    fprintf (fp_rw, "POLYGONS %d %d\n", getTriangleCount(), 5*getTriangleCount());
    i = 0;
    for (t = tri; t < tend; ++t)
    {
        fprintf (fp_rw, "%d %d %d %d %d\n", 4, 4*i, 4*i+1, 4*i+2, 4*i+3);
        i++;
    }
    
    fprintf (fp_rw, "CELL_DATA %d\n", getTriangleCount());
    fprintf (fp_rw, "SCALARS index int 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");      
    i = 0;
    for (t = tri; t < tend; ++t)
      fprintf (fp_rw, "%d\n", i++);
    
    fclose(fp_rw);
#endif 
    
#if defined(L2)

    int nbcirclepts = 40;
    sprintf(filename, "test_circle_%d.vtk", count);
    printf(" write file --> %s\n", filename);

    fp_rw = fopen (filename, "w");

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET POLYDATA\n");
    fprintf (fp_rw, "POINTS %d float\n", nbcirclepts*getTriangleCount());

    for (t = tri; t < tend; ++t)
    {
#if defined(POWER)
	double c[3] = {get_center(t)[0], get_center(t)[1], get_center(t)[2]+get_vertex(t->get_vertex(0))->w};
#else
	double c[3] = {get_center(t)[0], get_center(t)[1], get_center(t)[2]};
#endif
// 	printf("--> %lf %lf %lf\n", c[0], c[1], c[2]);
	for (int i = 0; i < nbcirclepts; i++)
        {
            fprintf(fp_rw, "%lf %lf %lf\n", c[0]+sqrt(c[2])*cos((double)i*2.0*M_PI/(double)nbcirclepts), c[1]+sqrt(c[2])*sin((double)i*2.0*M_PI/(double)nbcirclepts), 0.0);
        }
    }

    fprintf (fp_rw, "POLYGONS %d %d\n", getTriangleCount(), (nbcirclepts+1)*getTriangleCount());
    int circ = 0;
    for (t = tri; t < tend; ++t)
    {
      fprintf (fp_rw, "%d", nbcirclepts);
      for (int i = 0; i < nbcirclepts; i++)
        fprintf (fp_rw, " %d", circ*nbcirclepts + i);
      fprintf (fp_rw, "\n");
      circ++;
    }
    
    fclose(fp_rw);
#endif 
    
#endif    
    
    count++;

    printf("*** Export mesh END ***\n");
}

void Triangulation::align (double *vx, int size, double &xmin, double &xmax, double &ymin, double &ymax)
{
    // attention alignement --> arrondis possible !
    double *s, *send; // pointer used to travel on |s|
    
    double dx = xmax-xmin;
    double dy = ymax-ymin;
    double factor = 1.0;
    
    if (dx > dy)
    {
      _factor = 1e15/dx;
    }
    else
    {
      _factor = 1e15/dy;
    }
    
    _xmin = xmin;
    _ymin = ymin;
    
    xmin = 0.0;
    xmax = (xmax-_xmin)*_factor;
    ymin = 0.0;
    ymax = (ymax-_ymin)*_factor;
    
    s = vx;
    
#if defined (POWER)
    send = vx + 3*size; // send est un pointeur sur la derniere case du tableau de pointeur
#else
    send = vx + 2*size;
#endif
    
    while (s < send)
    {
        *s = round((*s-_xmin)*_factor);
        ++s;
	*s = round((*s-_ymin)*_factor);
        ++s;
#if defined (POWER)
	++s;
#endif	
    }    
}

void Triangulation::disalign()
{
    for (int i = 0; i < getVertexCount(); i++)
    {
        Vertex *pv = get_vertex(i);
        pv->p[0] = (pv->p[0]/_factor)+_xmin;
        pv->p[1] = (pv->p[1]/_factor)+_ymin;
    }

    Triangle *tri = getTri();
    Triangle *tend = tri + getTriangleCount();
    Triangle *t;

    for (t = tri; t < tend; ++t)
    {
        double c[3] = {get_center(t)[0], get_center(t)[1], get_center(t)[2]};

#if  defined(Linf)
        t->set_square_center((c[0]/_factor)+_xmin, (c[1]/_factor)+_ymin, c[2]/_factor, c[2]/_factor);
#endif
#if defined(L2)
        t->set_circle_center((c[0]/_factor)+_xmin, (c[1]/_factor)+_ymin, c[2]/_factor);
#endif
    }
}

void Triangulation::disalign(double *p)
{
    p[0] = (p[0]/_factor)+_xmin;
    p[1] = (p[1]/_factor)+_ymin;
}

void Triangulation::init (double *vx, int n, double epsilon, bool regular)
{
    /*!
     *  Implementation:
     *     - find the bounding box enlarged by |epsilon|
     *     - initialize vertex and triangle arrays.
     *     - initialize bucket (each point is put into a bucket).
     *     - initialize the triangulation by the two triangles of the bounding box.
     *     - insert vertices into triangulation.
     */   
    
// #if defined(ANGLE)
//     for (int i = 0; i < n; i++)
//     {
// 	double x = vx[2*i];
// 	double y = vx[2*i+1];
// 	vx[2*i] = x*cos(_angle) + y*sin(_angle);
// 	vx[2*i+1] = -x*sin(_angle) + y*cos(_angle);
//     }
// #endif

    double *s, *send; // pointer used to travel on |s|
    double xmin, xmax, ymin, ymax; // bounding box
#if defined (POWER)
    double wmax;
#endif
    double dx, dy; // extra size enlarged by |epsilon|
    
    srand(0);

    s = vx; // t est donc un pointeur sur le tableau de vertices
#if defined (POWER)
//     check_power(vx, n);
    send = vx + 3*n /* previously: + ptind(n)*/; // send est un pointeur sur la derniere case du tableau de pointeur
    xmin = *s++;
    ymin = *s++;
    wmax = *s++;
    xmax = xmin;
    ymax = ymin;
    
#else
    send = vx + 2*n /* previously: + ptind(n)*/;
    xmin = *s++;
    ymin = *s++;
    xmax = xmin;
    ymax = ymin;
#endif

    while (s < send) // creation de la boite englobante, c'est a dire que boucle sur tous les sommets (t -> t+2*n)
    {
        if (*s < xmin)
            xmin = *s;
        else if (*s > xmax)
            xmax = *s;
        ++s;
        if (*s < ymin)
            ymin = *s;
        else if (*s > ymax)
            ymax = *s;
        ++s;
#if defined (POWER)
        if (*s > wmax)
            wmax = *s;
	++s;
#endif
    }
    
#if defined (ALIGN) 
     align (vx, n, xmin, xmax, ymin, ymax);    
#endif
     
//   printf("Build Bounding Box index=%d epsilon=%f\n", index, epsilon);

    dx = (xmax - xmin) * epsilon; // ajout du facteur epsilon a la taille de la boite englobante
    xmin -= dx;
    xmax += dx;
    dy = (ymax - ymin) * epsilon;
    ymin -= dy;
    ymax += dy;

//   printf("Build dx,dy %lf %lf %lf %lf %lf %lf\n", xmin, xmax, ymin, ymax, dx, dy);

    v_capacity = n + 4; // le tableau de vertex a maintenant 4 vertices de plus
//     v_capacity *= 1000;
    v = new Vertex[v_capacity];
    tri_capacity = (n + 1) * 2 /*<<1*/;
//     tri_capacity *= 1000;
    tri = new Triangle[tri_capacity];
    
    free_t = new Triangle*[FREE_CAPACITY];

//   printf("Build Vertex and Triangle Array\n");

    llx = xmin - dx; // construction du bucket
    lly = ymin - dy;
    wr = BUCKET_SIZE / (xmax + dx - llx);
    hr = BUCKET_SIZE / (ymax + dy - lly);
    bzero (bucket_array, BUCKET_SIZE * BUCKET_SIZE * sizeof (Vertex *));

#if defined (POWER)
//     printf("wmax = %lf\n", wmax);
    new_vertex (xmin, ymin, wmax); // construction des 4 sommets supplementaires (mais ne les place pas dans le bucket)
    new_vertex (xmax, ymin, wmax);
    new_vertex (xmax, ymax, wmax);
    new_vertex (xmin, ymax, wmax);

    //   printf("Build VertexBucket\n");

    send = vx + 3*n /* previously: + ptind(n)*/;
    for (s = vx; s < send; ++s, ++s, ++s)
        new_vertex (s); // ajoute tous les sommets du nuage dans le bucket
#else
    new_vertex (xmin, ymin); // construction des 4 sommets supplementaires (mais ne les place pas dans le bucket)
    new_vertex (xmax, ymin);
    new_vertex (xmax, ymax);
    new_vertex (xmin, ymax);

    //   printf("Build VertexBucket\n");

    send = vx + 2*n /* previously: + ptind(n)*/;
    for (s = vx; s < send; ++s, ++s)
        new_vertex (s); // ajoute tous les sommets du nuage dans le bucket
#endif

//   printf("Build Vertices\n");

    Triangle *t0, *t1;
    t0 = new_triangle ();
    t1 = new_triangle ();

//   printf("Build Two First Triangles\n");

    t0->set_members (0, 2, 3, NULL, NULL, t1, -1, -1, 1); // v0, v1, v2, adj0, adj1, adj2, opp0, opp1, opp2
    t1->set_members (0, 1, 2, NULL, t0, NULL, -1, 2, -1);
    
//     ////////////
//     double a0 = (get_vertex(0)->a + get_vertex(2)->a + get_vertex(3)->a)/3.0;
//     double a1 = (get_vertex(0)->a + get_vertex(1)->a + get_vertex(2)->a)/3.0;
//     
//     t0->set_angle(a0, cos(a0), sin(a0));
//     t1->set_angle(a1, cos(a1), sin(a1));
//     ///////////
    
    
    get_vertex(0)->t = t0;
    get_vertex(1)->t = t1;
    get_vertex(2)->t = t0;
    get_vertex(3)->t = t0;

    t0->set_modified ();
    t1->set_modified ();

//   printf("Modify Triangles\n");

    last_tri = tri;
    last_i = 0;
    Vertex **b = bucket_array;
    int x = 1;
    //////////////////////////////////////////////////////////////////////////////////////////////////////////
    //                                    IMPORTANT                                                         //
    //////////////////////////////////////////////////////////////////////////////////////////////////////////
    // Boucle sur le bucket ET NON SUR LE TABLEAU DE POINTS. Permet de reutiliser le triangle precedemment
    // trouve pour l'insertion du nouveau sommet (variables : last_tri et last_i). la boucle sur le bucket
    // est telle que l'on passe continuellement d'une case voisine a une autre.
    //////////////////////////////////////////////////////////////////////////////////////////////////////////

    int count = 0;
    for (int j = 0; j < BUCKET_SIZE - 1; ++j, x = -x, b += BUCKET_SIZE)
    {
        for (int i = 0; i < BUCKET_SIZE; ++i, b += x)
        {
            for (Vertex * pv = *b; pv != NULL; pv = pv->n)
            {
		set_angle(pv->a);
                ////////////////////////////////////////////////
                ///                Insertion                   //
                ////////////////////////////////////////////////
                dprintf("insert number --> %d\n", count);
//                 if (count % 10000 == 0)
// //                     if (count >= 701930)
// 		    {
// 		        printf("insert number --> %d\n", count);
//                         check_tri();
// // 			getchar();
// 		    }

#if defined(STEP_DEBUG)
                int step_debug = 1;
                int range_debug[2] = {-1, -1};

		 step_debug = 1517;
		 range_debug[0] = step_debug-1;
		 range_debug[1] = step_debug+1;
// 		 range_debug[0] = 9470;
// 		 range_debug[1] = 9480;
// 		 range_debug[0] = 0;
// 		 range_debug[1] = 10;
// 		 step_debug = count;
                if (count >= range_debug[0] && count <= range_debug[1])
                {
                    printf("insert number --> %d\n", count);
                    export_vtk();
                    if (count == step_debug)
                    {
// 		  export_vtk();
                        dump_steps = 1;
                    }
                    else
                        dump_steps = 0;
                }
                else
                    dump_steps = 0;
#endif
		
#if defined(LAWSON_ALGORITHM)
                if (!insert_lawson (pv)) // main part, we insert each vertex
#endif
#if defined(BOWYERWATSON_ALGORITHM)
                if (!insert_bowyer_watson (pv)) // main part, we insert each vertex
#endif	
                {
                    printf("Exiting at %d !! \n", count);
                    return;
                }
                export_vtk();

//                 check_delaunay();
                ////////////////////////////////////////////////
#if defined(STEP_DEBUG)                
                if (count == range_debug[1])
                    return;
#endif		
// 		if (count == 3)
// 		  return;
		count++;
            }
        }
    }
    dprintf("free size --> %d\n", ifree_t);
    
    check_delaunay();
    
// #if defined (ALIGN)    
//      disalign();
// #endif
    return;
}

const int Triangle::p[3] = { 2, 0, 1 };
const int Triangle::n[3] = { 1, 2, 0 };

bool Triangulation::insert_lawson (Vertex * pv) // insert a point into an existing triangulation
{
    Triangle *t, *s;
    int i, j, ni;
    int stack = 0;

    t = last_tri; // pointeur sur le dernier triangle visite
    i = last_i; // dernier indice trouve dans last_tri
    
//     printf("pv --> %lf %lf : bbox (%lf %lf, %lf %lf)\n", pv->p[0], pv->p[1], get_vertex(0)->p[0], get_vertex(0)->p[1], get_vertex(2)->p[0], get_vertex(2)->p[1]);
    assert(pv->p[0] > get_vertex(0)->p[0]);
    assert(pv->p[1] > get_vertex(0)->p[1]);
    assert(pv->p[0] < get_vertex(2)->p[0]);
    assert(pv->p[1] < get_vertex(2)->p[1]);

    int p = locate (pv->p, t, i);
    /////////////////////////////////////////////////////////////////////////////
    // pv->p : les deux coordonnees du vertex, p est l'indice du triangle trouve
    // (mais t et i sont modifies)
    // valeurs de retour :
    // -1 : le sommet est sur une arete
    // 0 : le sommet est confondu avec un sommet existant
    // >0 : l'indice du triangle contenant le sommet
    /////////////////////////////////////////////////////////////////////////////

    if (p == 0) // cas d'un sommet confondu
    {
// 	printf("--> insert vertex on vertex\n");
        return true;
    }
    else if (p < 0) // cas d'un sommet sur une arete
    {
// 	printf("--> insert vertex on edge\n");
        // pv est le pointeur vers le sommet, v est le pointeur du tableau de sommets de la triangulation
        // pv - v retourne donc l'indice entier du sommet
        p = pv - v;

        // Renvoie s le triangle adjacent par la ieme arete de t. j est l'indice de cette meme arete dans s.
        get_nb (t, i, s, j); // "get_nb" signifie "get_neighbor".

        ///////////////////////////////////////////////////////////////////////
        //                            Split                                 //
        ///////////////////////////////////////////////////////////////////////
        // If |p| is on the symmetric edges (t,i) and (s,j)
        // We first make a copy |t1| of |t| and a copy |s1| of |s|, then make
        // proper adjustments for each triangle. Finally, |i| is changed to be
        // the triangle index for |p| in |t|.
        ///////////////////////////////////////////////////////////////////////

        int pi = Triangle::prev (i); // on recupere just un indice (0, 1 ou 2)
        int ni = Triangle::next (i);
        int pj = Triangle::prev (j);
        int nj = Triangle::next (j);

        Triangle *t1 = copy (t); // cree une copie de t
        Triangle *s1 = copy (s);	
        change_id (t1, ni); // integre t1 comme adjacence de la ni-ieme adjacence de t
        change_id (s1, pj); // integre s1 comme adjacence de la pj-ieme adjacence de s
	t1->set_flag(1);
	s1->set_flag(1);

#if 0
        Triangle *&tmp = v[get_vertex (t, pi)].t;
        if (tmp == t || tmp == s)
            tmp = t1;
#endif

        /////////////////////////////////////////////////////////////////
        // rappel : i est l'indice de l'arete dans t
        // nous avons maintenant 4 triangles dans la "cavite"
        /////////////////////////////////////////////////////////////////

        set_vertex (t, pi, p); // place p comme le pi-ieme sommet de t
        set_vertex (t1, ni, p); // place p comme le ni-ieme sommet de t1
        set_vertex (s, nj, p); // place p comme le nj-ieme sommet de s
        set_vertex (s1, pj, p); // place p comme le pj-ieme sommet de s1

        set_nb (t, ni, t1, pi); // la ni-ieme adjacence de t est t1 && la pi-ieme adjacence de t1 est t
        set_nb (t1, i, s1, j); // la i-ieme adjacence de t1 est s1 && la j-ieme adjacence de s1 est t1
        set_nb (s, pj, s1, nj); // la pj-ieme adjacence de s est s1 && la nj-ieme adjacence de s1 est s

        t->set_modified ();
        s->set_modified ();
        t1->set_modified ();
        s1->set_modified ();
	
        i = pi;
        stack = 4; // taille de la stack contenant les aretes de la cavite
    }
    else // cas general d'un sommet dans un triangle
    {
//         printf("--> insert vertex on triangle\n");
        p = pv - v;
	
        // Splitting a triangle is similar to splitting an edge.
        // On etoile le triangle t a partir du point insere p
        Triangle *t1 = copy (t);
        change_id (t1, 2);
	t1->set_flag(1);
        Triangle *t2 = copy (t);
        change_id (t2, 1);
	t2->set_flag(1);
#if 0
        Triangle *&tmp = v[get_vertex (t, 0)].t;
        if (tmp == t)
            tmp = t1;
#endif

        set_vertex (t, 0, p);
        set_vertex (t2, 1, p);
        set_vertex (t1, 2, p);


        set_nb (t, 2, t1, 0);
        set_nb (t, 1, t2, 0);
        set_nb (t1, 1, t2, 2);

        t->set_modified ();
        t1->set_modified ();
        t2->set_modified ();
		
//         dump_conf(this, t, t2, t1);
//         assert(volume(t)>0);
//         assert(volume(t1)>0);
//         assert(volume(t2)>0);
	
        i = 0;
        stack = 3;
    }

    ////////////////////////////////////////////////////////////////////////////////////
    //                                Corrections                                     //
    ////////////////////////////////////////////////////////////////////////////////////
    // A ce stade, nous disposons d'un ensemble de triangles flagges, d'un nombre d'aretes
    // dont la delaunay admissibilite doit etre verifiee et eventuellement corrigee.
    ////////////////////////////////////////////////////////////////////////////////////

    ni = Triangle::next (i);

#if defined(L2)
    do
    {
        get_nb (t, i, s, j); // recupere le triangle voisin de t
// 	printf("stack = %d\n", stack);
	
	debug_step(t, s);
        if ((s != NULL) && (incircle (s, p) > 0) && flip (t, i, s, j, stack))
            {
// 	      printf("flip !!\n");
            }
            else
            {
                --stack;
                s = get_adjacent (t, ni);
                i = Triangle::next (get_index (t, ni));
                ni = Triangle::next (i);
                t = s;
            }
        assert(stack < 10000);
    }
    while (stack > 0);
    return true;
    
#endif
    
#if  defined(Linf)
    do
    {
// 	printf("p --> %d\n", p);
        get_nb (t, i, s, j); // recupere le triangle voisin de t
// 	printf("stack = %d\n", stack);

        debug_step(t, s);
//         if ((s != NULL) && (insquare (s, p) >= 0 || insquare (t, s->get_vertex(j)) >= 0) && flip (t, i, s, j, stack))
	if ((s != NULL) && mutual_insquare(t, i, s, j) && flip (t, i, s, j, stack))
        {
            while (check_previous(t, i, stack));
            assert(t->get_vertex(i) == p);
        }
        else
        {
            while (check_previous(t, i, stack));
            assert(t->get_vertex(i) == p);
            --stack;
            s = get_adjacent (t, ni);
            i = Triangle::next (get_index (t, ni));
            ni = Triangle::next (i);
            t = s;
        }
        assert(stack < 10000);

        if (stack == 0)
            while (check_previous(t, i, stack));
    }
    while (stack > 0);            
    return true;
#endif
    
}

inline bool Triangulation::check_previous(Triangle * t, int i, int &stack)
{
//     printf("check previous\n");
    Triangle *c;
    int ii = Triangle::prev (i);
    int jj;
    get_nb (t, ii, c, jj); // recupere le triangle voisin de t
    
    debug_step(t, c);
/*    std::cerr << "triangle 0 (";
    print_triangle (std::cerr, c) << ")\n";*/ 
    
//     if ((c != NULL) && (insquare (t, c->get_vertex(jj)) >= 0 || insquare (c, t->get_vertex(ii))  >= 0 ) && flip (t, ii, c, jj, stack))
    if ((c != NULL) && mutual_insquare(t, ii, c, jj) && flip (t, ii, c, jj, stack))
    {
// 	    std::cerr << "triangle 1 (";
//     print_triangle (std::cerr, c) << ")\n";
        check_neigh (t, i, stack);
        check_neigh (t, Triangle::next(i), stack);
        return true;
    }
    return false;
}

inline bool Triangulation::flip (Triangle * t, int i, Triangle * s, int j, int &stack)
{   
 
//   printf("check classical\n");
// #if defined(Linf)
    if (!checkflip2to2(t, i, s, j))
      return false;
// #endif    
    ++stack;
    
    if (dump_steps)
      printf("--> flip !!\n");
    
 #ifdef DEBUG    
    printf("tri before flip : v0(%lf, %lf) v1(%lf, %lf) v2(%lf, %lf)\n", v[t->get_vertex(0)].p[0], v[t->get_vertex(0)].p[1],
           v[t->get_vertex(1)].p[0], v[t->get_vertex(1)].p[1],
           v[t->get_vertex(2)].p[0], v[t->get_vertex(2)].p[1]);
    printf("tri before flip : v0(%lf, %lf) v1(%lf, %lf) v2(%lf, %lf)\n", v[s->get_vertex(0)].p[0], v[s->get_vertex(0)].p[1],
           v[s->get_vertex(1)].p[0], v[s->get_vertex(1)].p[1],
           v[s->get_vertex(2)].p[0], v[s->get_vertex(2)].p[1]);
#endif   
    
    flip2to2 (t, i, s, j);
    
#ifdef DEBUG    
    printf("tri after flip : v0(%lf, %lf) v1(%lf, %lf) v2(%lf, %lf)\n", v[t->get_vertex(0)].p[0], v[t->get_vertex(0)].p[1],
           v[t->get_vertex(1)].p[0], v[t->get_vertex(1)].p[1],
           v[t->get_vertex(2)].p[0], v[t->get_vertex(2)].p[1]);    
    printf("tri after flip : v0(%lf, %lf) v1(%lf, %lf) v2(%lf, %lf)\n", v[s->get_vertex(0)].p[0], v[s->get_vertex(0)].p[1],
           v[s->get_vertex(1)].p[0], v[s->get_vertex(1)].p[1],
           v[s->get_vertex(2)].p[0], v[s->get_vertex(2)].p[1]);
#endif
        
    assert(volume(s)>=0);
    assert(volume(t)>=0);
    
    return true;
}

bool Triangulation::check_neigh (Triangle * t, int i, int &stack)
{
    Triangle *tt;
    int ii;
            
    get_nb (t, i, tt, ii);
    
    if (tt == NULL)
      return false;
        
    int i0 = ii;
    ii = Triangle::next(ii);

    while (ii != i0)
    {
        Triangle *ss = NULL;
        int jj;

        get_nb (tt, ii, ss, jj);
	
	debug_step(tt, ss);
	
//         printf("check cocircular !\n");
//         if ((ss != NULL) && (insquare (ss, tt->get_vertex(ii)) >= 0 || insquare (tt, ss->get_vertex(jj)) >= 0) && flip (tt, ii, ss, jj, stack))
	if ((ss != NULL) && mutual_insquare(tt, ii, ss, jj) && flip (tt, ii, ss, jj, stack))
        {
		check_neigh (tt, ii, stack);
		check_neigh (tt, Triangle::next(ii), stack);
        }
        else
            ii = Triangle::next(ii);
    }
    return true;
}

#if defined(ANGLE)
bool Triangulation::checkflip2to2 (Triangle * t, int i, int p, double angle)
#else
bool Triangulation::checkflip2to2 (Triangle * t, int i, int p)
#endif
{
  //     printf("check\n");
    Triangle *s;
    s = new Triangle;
    s->set_members(t->get_vertex(Triangle::prev (i)), t->get_vertex(Triangle::next (i)), p, NULL, NULL, NULL, -1, -1, -1);
    s->set_modified();
#if defined(Linf)
    Triangle *t1;
    t1 = new Triangle;
    Triangle *s1;
    s1 = new Triangle; 
    t1->set_members(t->get_vertex(i), t->get_vertex(Triangle::next (i)), p, NULL, NULL, NULL, -1, -1, -1);
    s1->set_members(p, t->get_vertex(Triangle::prev (i)), t->get_vertex (i), NULL, NULL, NULL, -1, -1, -1);
//     t1->set_modified();
//     s1->set_modified();
//           
//     // check square sizes   
//     const double *ct0 = get_center(t);
//     const double *cs0 = get_center(s);
//     const double *ct1 = get_center(t1);
//     const double *cs1 = get_center(s1);
    
    double ct0[4];
    double cs0[4];
    double ct1[4];
    double cs1[4];
    
#if defined(ANGLE)
    circumsquare(t, ct0, angle);
    circumsquare(s, cs0, angle);
    circumsquare(t1, ct1, angle);
    circumsquare(s1, cs1, angle);
#else
    circumsquare(t, ct0);
    circumsquare(s, cs0);
    circumsquare(t1, ct1);
    circumsquare(s1, cs1);
#endif    

    if (ct0[0] != INFINITY && cs0[0] != INFINITY)
    {
        if ((ct0[0] == cs0[0]) && (ct0[1] == cs0[1]) && (ct1[0] == cs1[0]) && (ct1[1] == cs1[1])) // cocyclicity : same center
        {
            delete s;
            delete t1;
            delete s1;
            return false;
        }

        if ((ct0[2] == ct1[2]) && (cs0[2] == cs1[2])) // cocyclicity : same radius
        {
            delete s;
            delete t1;
            delete s1;
            return false;
        }

        if ((ct0[2] == cs1[2]) && (cs0[2] == ct1[2])) // cocyclicity : same radius
        { 
            delete s;
            delete t1;
            delete s1;
            return false;
        }
    }
    
    double at0 = ct0[2]*ct0[3];
    double as0 = cs0[2]*cs0[3];
    double at1 = ct1[2]*ct1[3];
    double as1 = cs1[2]*cs1[3];
    
    double sa0 = at0 + as0;
    double sa1 = at1 + as1;    
        
    if (dump_steps)
    {
        printf("old : (%lf %lf) %lf --> %lf | (%lf %lf) %lf --> %lf\n", ct0[0], ct0[1], ct0[2], at0, cs0[0], cs0[1], cs0[2], as0);
        printf("new : (%lf %lf) %lf --> %lf | (%lf %lf) %lf --> %lf\n", ct1[0], ct1[1], ct1[2], at1, cs1[0], cs1[1], cs1[2], as1);
        printf("at0 --> %lf | as0 --> %lf\n", at0, as0);
        printf("at1 --> %lf | as1 --> %lf\n", at1, as1);
	printf("sa0 --> %lf | sa1 --> %lf ==> diff = %lf\n", sa0, sa1, sa1-sa0); // diff can be truncated using alignement
    }
    
    if (sa1 > sa0) // wrong swap
    {
        delete s; 
	delete t1;
        delete s1;
        return false;
    }

    if ((ct0[0] == INFINITY) && (ct1[0] == INFINITY) && (as1 > as0))
    {
        delete s; 
	delete t1;
        delete s1;
        return false;
    }
    if ((ct0[0] == INFINITY) && (cs1[0] == INFINITY) && (at1 > as0))
    {
        delete s; 
	delete t1;
        delete s1;
        return false;
    }
    if ((cs0[0] == INFINITY) && (cs1[0] == INFINITY) && (at1 > at0))
    {
        delete s; 
	delete t1;
        delete s1;
        return false;
    }
    if ((cs0[0] == INFINITY) && (ct1[0] == INFINITY) && (as1 > at0))
    {
        delete s; 
	delete t1;
        delete s1;
        return false;
    }
    
    delete t1;
    delete s1;
        
    // check volumes
#endif
    
    Vertex *v[3];
    v[0] = get_vertex(t->get_vertex(i));
    v[1] = get_vertex(t->get_vertex(Triangle::next (i)));
    v[2] = get_vertex(p);    
    if ( ((v[1]->p[0]-v[0]->p[0])*(v[2]->p[1]-v[0]->p[1]) - (v[1]->p[1]-v[0]->p[1])*(v[2]->p[0]-v[0]->p[0])) <= 0.0)
    {
	delete s; 
        return false;
    }

    v[0] = get_vertex(p);
    v[1] = get_vertex(t->get_vertex(Triangle::prev (i)));
    v[2] = get_vertex(t->get_vertex (i));
    if ( ((v[1]->p[0]-v[0]->p[0])*(v[2]->p[1]-v[0]->p[1]) - (v[1]->p[1]-v[0]->p[1])*(v[2]->p[0]-v[0]->p[0])) <= 0.0)
    {
        delete s;
        return false;
    }
    
    delete s;
    return true;
}

bool Triangulation::checkflip2to2 (Triangle * t, int i, Triangle * s, int j)
{
//     printf("check\n");
#if defined(Linf)
    Triangle *t1;
    t1 = new Triangle;
    Triangle *s1;
    s1 = new Triangle; 
    t1->set_members(t->get_vertex(i), t->get_vertex(Triangle::next (i)), s->get_vertex (j), NULL, NULL, NULL, -1, -1, -1);
    s1->set_members(s->get_vertex(j), s->get_vertex(Triangle::next (j)), t->get_vertex (i), NULL, NULL, NULL, -1, -1, -1);
    t1->set_modified();
    s1->set_modified();
    
      
    // check square sizes
   
    const double *ct0 = get_center(t);
    const double *cs0 = get_center(s);
    const double *ct1 = get_center(t1);
    const double *cs1 = get_center(s1);

    if (ct0[0] != INFINITY && cs0[0] != INFINITY)
    {
        if ((ct0[0] == cs0[0]) && (ct0[1] == cs0[1]) && (ct1[0] == cs1[0]) && (ct1[1] == cs1[1])) // cocyclicity : same center
        {
	    delete t1; delete s1;
            return false;
        }

        if ((ct0[2] == ct1[2]) && (cs0[2] == cs1[2])) // cocyclicity : same radius
        {
	    delete t1; delete s1;
            return false;
        }

        if ((ct0[2] == cs1[2]) && (cs0[2] == ct1[2])) // cocyclicity : same radius
        { 
	    delete t1; delete s1;
            return false;
        }
    }
    
    double at0 = ct0[2]*ct0[3];
    double as0 = cs0[2]*cs0[3];
    double at1 = ct1[2]*ct1[3];
    double as1 = cs1[2]*cs1[3];
    
    double sa0 = at0 + as0;
    double sa1 = at1 + as1;    
        
    if (dump_steps)
    {
        printf("old : (%lf %lf) %lf --> %lf | (%lf %lf) %lf --> %lf\n", ct0[0], ct0[1], ct0[2], at0, cs0[0], cs0[1], cs0[2], as0);
        printf("new : (%lf %lf) %lf --> %lf | (%lf %lf) %lf --> %lf\n", ct1[0], ct1[1], ct1[2], at1, cs1[0], cs1[1], cs1[2], as1);
        printf("at0 --> %lf | as0 --> %lf\n", at0, as0);
        printf("at1 --> %lf | as1 --> %lf\n", at1, as1);
	printf("sa0 --> %lf | sa1 --> %lf ==> diff = %lf\n", sa0, sa1, sa1-sa0); // diff can be truncated using alignement
    }
    
    if (sa1 > sa0) // wrong swap
    {
        delete t1;
        delete s1;
        return false;
    }

    if ((ct0[0] == INFINITY) && (ct1[0] == INFINITY) && (as1 > as0))
    {
        delete t1;
        delete s1;
        return false;
    }
    if ((ct0[0] == INFINITY) && (cs1[0] == INFINITY) && (at1 > as0))
    {
        delete t1;
        delete s1;
        return false;
    }
    if ((cs0[0] == INFINITY) && (cs1[0] == INFINITY) && (at1 > at0))
    {
        delete t1;
        delete s1;
        return false;
    }
    if ((cs0[0] == INFINITY) && (ct1[0] == INFINITY) && (as1 > at0))
    {
        delete t1;
        delete s1;
        return false;
    }
    
    delete t1;
    delete s1;
        
    // check volumes
#endif
    
    Vertex *v[3];
    v[0] = get_vertex(t->get_vertex(i));
    v[1] = get_vertex(t->get_vertex(Triangle::next (i)));
    v[2] = get_vertex(s->get_vertex (j));
    if ( ((v[1]->p[0]-v[0]->p[0])*(v[2]->p[1]-v[0]->p[1]) - (v[1]->p[1]-v[0]->p[1])*(v[2]->p[0]-v[0]->p[0])) <= 0.0)
        return false;

    v[0] = get_vertex(s->get_vertex(j));
    v[1] = get_vertex(s->get_vertex(Triangle::next (j)));
    v[2] = get_vertex(t->get_vertex (i));
    if ( ((v[1]->p[0]-v[0]->p[0])*(v[2]->p[1]-v[0]->p[1]) - (v[1]->p[1]-v[0]->p[1])*(v[2]->p[0]-v[0]->p[0])) <= 0.0)
        return false;

    return true;
}

void Triangulation::flip2to2 (Triangle * t, int i, Triangle * s, int j)
{
//     printf("flip !\n");
    Triangle *u, *v;
    int k, l;
    int ni, pi, nj, pj;
    nj = Triangle::next (j);
    pj = Triangle::prev (j);
    ni = Triangle::next (i);
    pi = Triangle::prev (i);
    get_nb (t, ni, u, k);
    get_nb (s, nj, v, l);
    set_nb (t, i, v, l);
    set_nb (s, j, u, k);
    set_nb (t, ni, s, nj);

#if 0
    Triangle *&tmpt = Triangulation::v[get_vertex (t, pi)].t;
    if (tmpt == t)
        tmpt = s;
    Triangle *&tmps = Triangulation::v[get_vertex (s, pj)].t;
    if (tmps == s)
        tmps = t;
#endif

    set_vertex (t, pi, get_vertex (s, j));
    set_vertex (s, pj, get_vertex (t, i));
    t->set_modified ();
    s->set_modified ();
}

int Triangulation::check_delaunay()
{
#if !defined(CHECK)
    return 0;
#endif
    printf("## --> check tri ##\n");
    Triangle *tri = getTri();
    Triangle *tend = tri + getTriangleCount();
    Triangle *t, *s;
    int error = 0;

    std::vector<Triangle *> tf;
    
    for (t = tri; t < tend; ++t)
    {
        assert(t->get_flag() != FREED_TRIANGLE_FLAG);
        t->set_flag(DEFAULT_TRIANGLE_FLAG);
    }

    for (t = tri; t < tend; ++t)
    {
        for (int i = 0; i < 3; i++)
        {
            int j;
            get_nb (t, i, s, j); // recupere le triangle voisin de t
#if defined(L2)
            if ((s != NULL) && (incircle (t, s->get_vertex(j)) > 0))
	    {
	      			dprintf("--> %.29lf\n", incircle (s, t->get_vertex(i)));

			t->set_flag(1);
                        error++;
// 			getchar();
	    }
#endif
#if defined(Linf)
            if ((s != NULL) && (insquare (s, t->get_vertex(i)) > 0 || insquare (t, s->get_vertex(j)) > 0))
            {
                if (checkflip2to2(t, i, s, j))
                {
                    dprintf("--> %.29lf ## %.29lf\n", insquare (s, t->get_vertex(i)), insquare (t, s->get_vertex(j)));

                    t->set_flag(1);
                    error++;
// 			getchar();
                }
            }
#endif
            
        }
    }
    printf("### Found %d delaunay errors\n", error);
//     assert(error == 0);
    return error; 
}

double Triangulation::volume(Triangle *t)
{
  Vertex *v[3];
  v[0] = get_vertex(t->get_vertex(0));
  v[1] = get_vertex(t->get_vertex(1));
  v[2] = get_vertex(t->get_vertex(2));
  return ((v[1]->p[0]-v[0]->p[0])*(v[2]->p[1]-v[0]->p[1]) - (v[1]->p[1]-v[0]->p[1])*(v[2]->p[0]-v[0]->p[0]))/2.0;
}


bool Triangulation::insert_bowyer_watson (Vertex * pv) // insert a point into an existing triangulation
{
  static int fv[3][2] = {{1,2},{2,0},{0,1}};
    Triangle *t, *s;
    int i;

    t = last_tri; // pointeur sur le dernier Triangle visite
    i = last_i; // dernier indice trouve dans last_tri

    dprintf("pv --> %lf %lf : bbox (%lf %lf, %lf %lf)\n", pv->p[0], pv->p[1], get_vertex(0)->p[0], get_vertex(0)->p[1], get_vertex(2)->p[0], get_vertex(2)->p[1]);
    
    int p = locate (pv->p, t, i);
    /////////////////////////////////////////////////////////////////////////////
    // valeurs de retour :
    // 0 : le sommet est confondu avec un sommet existant
    // -1 : le sommet est sur une arete
    // 1 : le sommet est dans le triangle
    /////////////////////////////////////////////////////////////////////////////
    
    // pv - v retourne donc l'indice entier du sommet
    int pvi = pv - v;    
    dprintf("t found = %d, cas %d\n", (int) (t - getTri()), p);
    
    if (p == 0) // cas d'un sommet confondu
    {
// 	printf("--> insert vertex on vertex\n");
        return true;
    }
    else
    {
        Triangle *base = t;

        cav->set_old_flag(DEFAULT_TRIANGLE_FLAG);
	if (p == -1) // arete
	  cav->set_base(base, i, pvi);
	else // triangle
	  cav->set_base(base, -1, pvi);

        cav->init();
        if (dump_steps)
            cav->export_vtk();
	
        build_initial_cavity(cav, pvi);

        cav->dump();
	cav->check();
        if (dump_steps)
            cav->export_vtk();
	cav->export_vtk();

        while (correct_cavity(cav, pvi))
        {
            cav->inc_old_flag();
            cav->init();
            if (dump_steps)
                cav->export_vtk();
            build_final_cavity(cav, pvi);
	    
	    cav->dump();
	    cav->check();
            if (dump_steps)
                cav->export_vtk();
	    cav->export_vtk();
        }
        
	finalize(cav);
        
        cav->dump();
        cav->check();
        if (dump_steps)
            cav->export_vtk();

        star_cavity(cav, pvi);

        // verif

//         export_vtk();
        check_adjacencies();

        return true;
    }
}

void Triangulation::finalize(Cavity *cav)
{    
    cav->init_read_vertex_cavity();
    Vertex *v = cav->next_vertex_cavity();
    while (v != NULL)
    {
	v->flag = 0;
        v = cav->next_vertex_cavity();
    }
    cav->clear_vertex_cavity();
    
    cav->init_read_triangle_cavity();
    Triangle *t = cav->next_triangle_cavity();
    while (t != NULL)
    {
	dprintf("--> erasing tetra %d from cavity\n", get_triangle_index(t));
	if (t->get_flag() == cav->get_new_flag())
	{
	  t->set_flag(FREED_TRIANGLE_FLAG);
	  free_triangle(t);
	}
	else
	{
	  t->set_flag(DEFAULT_TRIANGLE_FLAG);
	}
        t = cav->next_triangle_cavity();
    }
    cav->clear_triangle_cavity();
}

bool Triangulation::mutual_insquare(Triangle *t, int i, Triangle *s, int j)
{
    if ((insquare (t, s->get_vertex(j)) >= 0.0) || (insquare(s, t->get_vertex(i)) >= 0.0))
      return true;
    else
      return false;
}

#if  defined(ANGLE)
bool Triangulation::mutual_insquare(Triangle *t, int i, int p)
{
    Triangle *s;
    s = new Triangle;
    s->set_members(t->get_vertex(Triangle::prev (i)), t->get_vertex(Triangle::next (i)), p, NULL, NULL, NULL, -1, -1, -1);
    s->set_modified();
    if ((insquare (t, p, t->get_angle()) >= 0.0) || (insquare(s, t->get_vertex(i), _angle) >= 0.0))
//     if ((insquare (t, p) >= 0.0) || (insquare(s, t->get_vertex(i)) >= 0.0))
    {
      delete s;
      return true;
    }
    else
      delete s;
      return false;
}
#else
bool Triangulation::mutual_insquare(Triangle *t, int i, int p)
{
    Triangle *s;
    s = new Triangle;
    s->set_members(t->get_vertex(Triangle::prev (i)), t->get_vertex(Triangle::next (i)), p, NULL, NULL, NULL, -1, -1, -1);
    s->set_modified();
    if ((insquare (t, p) >= 0.0) || (insquare(s, t->get_vertex(i)) >= 0.0))
    {
      delete s;
      return true;
    }
    else
      delete s;
      return false;
}
#endif

bool Triangulation::build_initial_cavity(Cavity *cav, int p)
{
//     dprintf("Build cavity\n");    
    printf("Build cavity\n");    

    cav->init_read_facet();
    Facet *ft = cav->get_next_facet();

    while (ft != NULL)
    {
        if (ft->get_flag() == ACTIVE_FACET)
        {
            Triangle *t = ft->get_ext_triangle();
            if (t != NULL)
            {
		int face = ft->get_ext_face();
                dprintf("checking facet %d --> %d %d : t int %d, t ext %d\n", cav->get_index(ft),
                       ft->get_vertex(0), ft->get_vertex(1), get_triangle_index(ft->get_int_triangle()), get_triangle_index(ft->get_ext_triangle()));

                if (t->get_flag() == cav->get_old_flag())
                {
#if  defined(L2)
                    if ((incircle (t, p) >= 0.0))
                    {
		      if (checkflip2to2 (t, ft->get_ext_face(), p))
		      {
                        cav->star_facet(ft);
                        cav->add_triangle_cavity(t);
		      }
                    }
#endif 
#if  defined(Linf)
#if  defined(ANGLE)
		    if (mutual_insquare(t, face, p))
// 		    if (insquare(t,p, t->get_angle()) >= 0.0)
#else
		    if (insquare (t, p) >= 0.0)
#endif
                    {
// 		      if (checkflip2to2 (t, face, p))
		      {
                        cav->star_facet(ft);
                        cav->add_triangle_cavity(t);
		      }
                    }
#endif                  
                }
                else
                {
		  dprintf("int flag %d <---> ext flag %d\n", ft->get_int_triangle()->get_flag(), ft->get_ext_triangle()->get_flag());
		  if (ft->get_ext_triangle()->get_flag() == cav->get_new_flag())
		  cav->remove_facet(ft);
                }
                if (dump_steps)
		  cav->export_vtk();
            }
        }
        ft = cav->get_next_facet();
    }
}

bool Triangulation::build_final_cavity(Cavity *cav, int p)
{
    dprintf("Build cavity\n");    

    cav->init_read_facet();
    Facet *ft = cav->get_next_facet();

    while (ft != NULL)
    {
        if (ft->get_flag() == ACTIVE_FACET)
        {
            Triangle *t = ft->get_ext_triangle();
            if (t != NULL)
            {
                dprintf("checking facet %d --> %d %d : t int %d, t ext %d\n", cav->get_index(ft),
                       ft->get_vertex(0), ft->get_vertex(1), get_triangle_index(ft->get_int_triangle()), get_triangle_index(ft->get_ext_triangle()));

                if (t->get_flag() == cav->get_old_flag())
                {
                    cav->star_facet(ft);
                }
                else
                {
                    dprintf("int flag %d <---> ext flag %d\n", ft->get_int_triangle()->get_flag(), ft->get_ext_triangle()->get_flag());
                    if (ft->get_ext_triangle()->get_flag() == cav->get_new_flag())
                        cav->remove_facet(ft);
                }
                if (dump_steps)
		  cav->export_vtk();
            }
        }
        ft = cav->get_next_facet();
    }  
}

void Triangulation::ball(Vertex *v, int &nb, Triangle **tlist)
{
   nb = 0;
   Triangle *t = v->t;
   Triangle *s = NULL;
   tlist[nb++] = t;
   
    int f = -1;
    for (int i = 0; i < 3 && f == -1; i++)
        if (get_vertex(t->get_vertex(i)) == v)
            f = i;
    assert(f != -1);
    
    f = (f+1)%3;
    
    s = t->get_adjacent(f);
    f = t->get_opp(f);
    while (s != tlist[0])
    {
        f = (f+2)%3;
        tlist[nb++] = s;
	t = s;
	assert (get_vertex(t->get_vertex((f+2)%3)) == v);
        s = t->get_adjacent(f);
        f = t->get_opp(f);
    }
}

bool Triangulation::correct_cavity(Cavity *cav, int p)
{
   static int fv[3][2] = {{1,2},{2,0},{0,1}};
    dprintf("--> Correcting cavity\n");
    
    int correction = 0;
    Facet *ft = NULL;
    
    // alone vertex    
    cav->init_read_vertex_cavity();
    Vertex *v = cav->next_vertex_cavity();
    while (v != NULL)
    {
        v->flag = 0;
	v = cav->next_vertex_cavity();
    }
    dprintf("--> %d vertex in cavity\n", cav->vertex_cavity_size());
    
    cav->init_read_facet();    
    ft = cav->get_next_facet();
    while (ft != NULL)
    {
        if (ft->get_flag() == ACTIVE_FACET)
        {
            cav->inc_vertex_cavity_flag(ft->get_vertex(0));
            cav->inc_vertex_cavity_flag(ft->get_vertex(1));  
        }
        ft = cav->get_next_facet();
    }
    
    cav->init_read_vertex_cavity();
    v = cav->next_vertex_cavity();
    while (v != NULL)
    {
      dprintf("vertex %d --> flag = %d\n", get_vertex_index(v), v->flag);
        if (v->flag == 0)
        {
	   Triangle *tball[1000];
	    int nb;
	    ball(v, nb, tball);
	    Triangle *t = NULL;
	    for (int i = 0; i < nb; i++)
	    {
	      t = tball[i];
	      if (!cav->is_base(t) && t->get_flag() == cav->get_new_flag())
		break;
	    }
	    assert(t != NULL);

            if (t->get_flag() != LOCKED_TRIANGLE_FLAG)
            {
		assert( !cav->is_base(t));
		assert( t->get_flag() == cav->get_new_flag());		
		dprintf("--> locking Triangle %d\n", get_triangle_index(t));
                t->set_flag(LOCKED_TRIANGLE_FLAG);
//                 cav->add_unmodified_triangle(t);
                correction++;
                dprintf("correction needed for vertex %d --> flag = %d\n", get_vertex_index(v), v->flag);
            }
        }
        v = cav->next_vertex_cavity();
    }

// facet orientation    
    cav->init_read_facet();    
    ft = cav->get_next_facet();
    while (ft != NULL)
    {
        if (ft->get_flag() == ACTIVE_FACET)
        {
            Triangle *intt = ft->get_int_triangle();
	    assert(intt != NULL);
            if (intt->get_flag() == cav->get_new_flag())
            {
                int face = ft->get_int_face();
                dprintf("check orient facet %d :: Triangle %d --> face %d = %d %d\n", cav->get_index(ft), (int)(intt-getTri()), face, ft->get_vertex(0), ft->get_vertex(1));
                double test = lefton_test(get_vertex(p)->p, intt, face);
                if (test <= 0.0)
                {
		    dprintf("--> locking Triangle %d\n", get_triangle_index(intt));
                    intt->set_flag(LOCKED_TRIANGLE_FLAG);
		    correction++;
                }
            }
        }
        ft = cav->get_next_facet();
    }   
 
    printf("--> %d corrections\n", correction);
    return correction;
}

bool Triangulation::star_cavity(Cavity *cav, int p)
{
    dprintf("Staring cavity\n");

    static int fv[3][2] = {{1,2},{2,0},{0,1}};

    Triangle *tcav[CAVITY_CAPACITY];

    for (int i = 0; i < cav->get_nb_facet(); i++)
    {
        Facet *ft = cav->get_facet(i);
        if (ft->get_flag() == ACTIVE_FACET)
        {
	    dprintf("using facet %d\n", cav->get_index(ft));
            tcav[i] = new_triangle();
            Triangle *nt =  tcav[i];
            int v0, v1;
            Triangle *t;
            int f;

            ft->get_vertices(v0, v1);

            t = ft->get_ext_triangle();
            f = ft->get_ext_face();

            nt->set_vertex(0, v0);
            nt->set_vertex(1, v1);
            nt->set_vertex(2, p);
            nt->set_modified();
            nt->set_flag(DEFAULT_TRIANGLE_FLAG);
            v[v0].t = v[v1].t = v[p].t = nt;
	    
#if defined(ANGLE)
	   nt->xor_modified();
	   update_circumsquare(nt);
	   assert(nt->get_center()[0] != INFINITY);
#endif

            dprintf("created triangle %d --> %d %d %d\n", (int)(nt-getTri()),  nt->get_vertex(0), nt->get_vertex(1), nt->get_vertex(2));

	    assert(volume(nt) > 0.0);

            if (t != NULL)
            {
// 		printf("set adj : Triangle %d face %d (%d %d %d) <---> Triangle %d face %d (%d %d %d)\n", (int)(t-tet), f, t->get_vertex(fv[f][0]), t->get_vertex(fv[f][1]), t->get_vertex(fv[f][2]),
// 		       (int)(nt-tet), 3, nt->get_vertex(fv[3][0]), nt->get_vertex(fv[3][1]), nt->get_vertex(fv[3][2]));
                t->set_adjacent(f, nt);
                t->set_opp(f, 2);
                nt->set_adjacent(2, t);
                nt->set_opp(2, f);
		
		assert(t->get_vertex(fv[f][0]) + t->get_vertex(fv[f][1]) == nt->get_vertex(fv[2][0]) + nt->get_vertex(fv[2][1]));
            }
        }
    }

    for (int i = 0; i < cav->get_nb_facet(); i++)
    {
        Facet *ft = cav->get_facet(i);
        if (ft->get_flag() == ACTIVE_FACET)
        {
            Triangle *nt =  tcav[i];
            Triangle *t;
	    
	    int id = cav->get_index(ft->get_adj(0));
	    t = tcav[id];
	    set_nb(nt, 0, t, 1);
			
            last_tri = nt;
	   }
    }
    
//         for (int i = 0; i < cav->get_nb_facet(); i++)
//     {
//         Facet *ft = cav->get_facet(i);
//         if (ft->get_flag() == ACTIVE_FACET)
//         {
//             Triangle *nt =  tcav[i];
// 	    dprintf("created triangle %d --> %d %d %d, opp %d %d %d\n", (int)(nt-getTri()),  nt->get_vertex(0), nt->get_vertex(1), nt->get_vertex(2),  nt->get_opp(0), nt->get_opp(1), nt->get_opp(2));
// 	}
//     }
}

int Triangulation::check_adjacencies()
{
#if !defined(DEBUG)
    return 0;
#endif
    static int fv[3][2] = {{1,2},{2,0},{0,1}};
    Triangle *tinit = getTri();
    Triangle *tend = tinit + getTriangleCount();
    Triangle *t, *at;

    int pb = 0;

    for (t = tinit; t < tend; ++t)
    {
        if (t->get_flag() == DEFAULT_TRIANGLE_FLAG)
        {
//             std::cerr<<"Triangle "<<(int)(t-tet)<<" (";
//             print_triangle(std::cerr,t)<<")\n";
            for (int i = 0; i < 3; i++)
            {
	        at = t->get_adjacent(i);
                if (at != NULL)
                {
                    int j = t->get_opp(i);

                    if (at->get_adjacent(j) != t)
                    {
		      printf("(Triangle %d) %d %d <---> (Triangle %d) %d %d\n", (int)(t-getTri()), t->get_vertex(fv[i][0]), t->get_vertex(fv[i][1]),
                           (int)(at-getTri()), at->get_vertex(fv[j][0]), at->get_vertex(fv[j][1]) );
                        printf("problem !!\n");
                        pb++;
                    }
                }
            }
        }
    }
    assert(pb == 0);
}

void Cavity::init()
{    
//     dprintf("Setting base cavity --> tetra %d\n", tri->get_triangle_index(get_base()));

//     clear_triangle_cavity();
    clear_vertex_cavity();
    clear_cavity_facet();

    static int fv[3][2] = {{1,2},{2,0},{0,1}};

    Triangle *t = tbase;
    assert(t->get_flag() != LOCKED_TRIANGLE_FLAG);
    t->set_flag(get_new_flag());

    Facet *pf[3];

    add_vertex_cavity(t->get_vertex(0));
    add_vertex_cavity(t->get_vertex(1));
    add_vertex_cavity(t->get_vertex(2));

    for (int i = 0; i < 3; i++)
    {
        pf[i] = new_facet();
	pf[i]->set_children(NULL, NULL);
	pf[i]->set_parent(NULL);
	pf[i]->set_int_triangle(NULL, -1);
	pf[i]->set_ext_triangle(NULL, -1);
	pf[i]->set_flag(ACTIVE_FACET);
    }

    pf[0]->set_members(t->get_vertex(fv[0][0]), t->get_vertex(fv[0][1]), pf[fv[0][0]], pf[fv[0][1]], 1, 0);
    pf[1]->set_members(t->get_vertex(fv[1][0]), t->get_vertex(fv[1][1]), pf[fv[1][0]], pf[fv[1][1]], 1, 0);
    pf[2]->set_members(t->get_vertex(fv[2][0]), t->get_vertex(fv[2][1]), pf[fv[2][0]], pf[fv[2][1]], 1, 0);
    
    for (int i = 0; i < 3; i++)
    {
        Triangle *s = t->get_adjacent(i);
	pf[i]->set_int_triangle(t, i);
        if (s != NULL)
            pf[i]->set_ext_triangle(s, t->get_opp(i));
    }
    
    if (fbase != -1)
      star_facet(pf[fbase]);
    
    dump();
    check();
}

void Cavity::remove_facet(Facet *fti)
{
    dprintf("--> remove facet %d --> %d %d\n", get_index(fti), fti->get_vertex(0), fti->get_vertex(1));
    for (int i = 0; i < 2; i++)
    {
	Facet *adj = fti->get_adj(i);
	int opp = fti->get_opp(i);
	
        if (fti->get_vertex(i) == adj->get_vertex(opp))
        {
	  Facet *aa0 = fti->get_adj((i+1)%2);
	  int oo0 = fti->get_opp((i+1)%2);
	  
	  Facet *aa1 = adj->get_adj((opp+1)%2);
	  int oo1 = adj->get_opp((opp+1)%2);
	  
	  aa0->set_nb(oo0, aa1, oo1);
	  aa1->set_nb(oo1, aa0, oo0); 
	  
	  fti->set_flag(INACTIVE_FACET);
	  fti->set_members(-1, -1, NULL, NULL, -1, -1);
	  free_facet(fti);
	  
	  dprintf("--> remove facet %d\n", get_index(adj));
	  adj->set_flag(INACTIVE_FACET);
	  adj->set_members(-1, -1, NULL, NULL, -1, -1);
	  free_facet(adj);
	  break;
        }
    }

    if (fti->get_flag() == ACTIVE_FACET)
    {
      dump();
      export_vtk();
      check();
    }
    assert(fti->get_flag() == INACTIVE_FACET);
}

void Cavity::star_facet(Facet *fti) // used to enlarge cavity
{
    static int fv[3][2] = {{1,2},{2,0},{0,1}};
    static int v2[3] = {0, 1, 0};
    
    dprintf("Staring facet %d --> tri int %d, tri ext %d\n", get_index(fti), tri->get_triangle_index(fti->get_int_triangle()), tri->get_triangle_index(fti->get_ext_triangle()));
    
    Facet *pf[2];
    int vft[2];
    Facet *nft[2];
    int oft[2];
    
    assert(fti->get_int_triangle()->get_flag() != LOCKED_TRIANGLE_FLAG);
    assert(fti->get_int_triangle()->get_flag() != FREED_TRIANGLE_FLAG);
    assert(fti->get_ext_triangle()->get_flag() != LOCKED_TRIANGLE_FLAG);
    assert(fti->get_ext_triangle()->get_flag() != FREED_TRIANGLE_FLAG);
    
    fti->get_ext_triangle()->set_flag(get_new_flag());
        
    for (int i = 0; i < 2; i++)
    {
        dprintf("creating facet %d\n", icavfw);
        pf[i] = new_facet();
        pf[i]->set_children(NULL, NULL);
        pf[i]->set_parent(fti);
	pf[i]->set_int_triangle(NULL, -1);
	pf[i]->set_ext_triangle(NULL, -1);
        pf[i]->set_flag(ACTIVE_FACET);
    }
    
    fti->set_children(pf[0], pf[1]);
    fti->set_flag(INACTIVE_FACET);
    
    fti->get_vertices(vft[0], vft[1]);
    fti->get_neigh(nft[0], nft[1], oft[0], oft[1]);
    
    Triangle *t = fti->get_ext_triangle();
    int face = fti->get_ext_face();
    int p = t->get_vertex(face);
    
    add_vertex_cavity(p);
    
    dprintf("staring facet %d --> %d %d, adj %d %d\n", get_index(fti), fti->get_vertex(0), fti->get_vertex(1), get_index(fti->get_adj(0)), get_index(fti->get_adj(1)));
    
    Triangle *s = NULL;

    pf[0]->set_members(fti->get_vertex(0), p, pf[1], nft[1], 1, oft[1]);
    pf[0]->set_int_triangle(t, (face+1)%3);
    s = t->get_adjacent((face+1)%3);
    if (s != NULL)
        pf[0]->set_ext_triangle(s, t->get_opp((face+1)%3));
    nft[1]->set_nb(oft[1], pf[0], 1);
    
    dprintf("facet %d --> %d %d\n", get_index(pf[0]), pf[0]->get_vertex(0), pf[0]->get_vertex(1));
    dprintf("set int triangle %d, face %d --> %d %d\n", tri->get_triangle_index(t), (face+1)%3, t->get_vertex(fv[(face+1)%3][0]), t->get_vertex(fv[(face+1)%3][1]));

    pf[1]->set_members(p, fti->get_vertex(1), nft[0], pf[0], oft[0], 0);
    pf[1]->set_int_triangle(t, (face+2)%3);
    s = t->get_adjacent((face+2)%3);
    if (s != NULL)
        pf[1]->set_ext_triangle(s, t->get_opp((face+2)%3));
    nft[0]->set_nb(oft[0], pf[1], 0);

    dprintf("facet %d --> %d %d\n", get_index(pf[1]), pf[1]->get_vertex(0), pf[1]->get_vertex(1));
    dprintf("set int triangle %d, face %d --> %d %d\n", tri->get_triangle_index(t), (face+2)%3, t->get_vertex(fv[(face+2)%3][0]), t->get_vertex(fv[(face+2)%3][1]));
}

void Cavity::dump()
{
#if !defined(DEBUG)
    return;
#endif
    printf("Cavity --> %d facet\n", icavfw);

    for (int i = 0; i < get_nb_facet(); i++)
    {
        Facet *ft = get_facet(i);
        if (ft->get_flag() == ACTIVE_FACET)
        {
            int vft[2];
            Facet *nft[2];
            int oft[2];
            ft->get_vertices(vft[0], vft[1]);
            ft->get_neigh(nft[0], nft[1], oft[0], oft[1]);

            printf("facet %d : v0 = %d, adj0 = %d, opp0 = %d || v1 = %d, adj1 = %d, opp1 = %d\n", (int)(ft-cav_f),
                   vft[0], (int)(nft[0]-cav_f), oft[0], vft[1], (int)(nft[1]-cav_f), oft[1]);
        }
  }
  
  init_read_triangle_cavity();
  Triangle *t = next_triangle_cavity();
    while (t != NULL)
    {
      printf("Triangle %d has flag %d\n", tri->get_triangle_index(t), t->get_flag());
      t = next_triangle_cavity();
    }
}

void Cavity::check()
{
#if !defined(DEBUG)
    return;
#endif
    int pb = 0;
    for (int i = 0; i < get_nb_facet(); i++)
    {
        Facet *ft = get_facet(i);
        if (ft->get_flag() == ACTIVE_FACET)
        {
            for (int j = 0; j < 2; j++)
            {
                Facet *adj = ft->get_adj(j);
                int opp = ft->get_opp(j);
                if (ft->get_vertex(j) == adj->get_vertex(opp))
                {
                    printf("## /!\\ ## facet %d has same opposite vertex %d --> facet %d\n", get_index(ft), ft->get_vertex(j), get_index(adj));
                    pb++;
                }
            }
        }
    }
	  
    for (int i = 0; i < get_nb_facet(); i++)
    {
        Facet *ft = get_facet(i);
        if (ft->get_flag() == ACTIVE_FACET)
        {
            for (int j = 0; j < 2; j++)
            {
                Facet *aft = ft->get_adj(j);
                if (aft != NULL)
		{
                    if (aft->get_adj(ft->get_opp(j)) != ft)
		    {
                        printf("## /!\\ ## facet %d --> adj %d = %d, opp %d = %d --> %d\n", get_index(ft), j, get_index(aft), j, ft->get_opp(j), get_index(aft->get_adj(ft->get_opp(j))));
			pb++;
		    }
		}

            }
        }
    }
    if (pb)
      exit(1);
}

void Cavity::export_vtk()
{    
//     getchar();
    static int count = 0;
    char filename[100];
    int nb_tri;
        
    sprintf(filename, "cavity_%d.vtk", count);
    printf(" write file --> %s\n", filename);

    FILE *fp_rw = fopen (filename, "w");
    
    int tcount = 0;
    for (int i = 0; i < get_nb_facet(); i++)
    {
        Facet *ft = get_facet(i);
	if (ft->get_flag() == ACTIVE_FACET)
	  tcount++;
    }

//     printf(" %d vertices, %d triangles\n", tri->getVertexCount(), tcount);

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET UNSTRUCTURED_GRID\n");
//     fprintf (fp_rw, "DATASET POLYDATA\n");
    fprintf (fp_rw, "POINTS %d float\n", tri->getVertexCount()+1);

    Vertex *pv;

    for (int i = 0; i < tri->getVertexCount(); i++)
    {
        pv = tri->get_vertex(i);
        double p[2] = {pv->p[0], pv->p[1]};
// #if defined(ALIGN)
//         disalign(p);
// #endif
// #if defined(ANGLE)
// 
//         double x = p[0];
//         double y = p[1];
//         p[0] = x*cos(_angle) - y*sin(_angle);
//         p[1] = x*sin(_angle) + y*cos(_angle);
// #endif
        fprintf(fp_rw, "%lf %lf %lf\n", p[0], p[1], 0.0);
    }
    
    Vertex *c = tri->get_vertex(center);
    fprintf(fp_rw, "%lf %lf %lf\n", c->p[0], c->p[1], 0.0);
    
//     fprintf (fp_rw, "VERTICES %d %d\n", 1, 2);
//     fprintf(fp_rw, "%d %d\n", 1, tri->getVertexCount());

//     fprintf (fp_rw, "POLYGONS %d %d\n", tcount, 4*tcount);
//     for (int i = 0; i < get_nb_facet(); i++)
//     {
//         Facet *ft = get_facet(i);
//         if (ft->get_flag() == ACTIVE_FACET)
// 	{
// 	  fprintf (fp_rw, "%d %d %d %d\n", 3, ft->get_vertex(0), ft->get_vertex(1), ft->get_vertex(2));
// 	}
//     }

    fprintf (fp_rw, "CELLS %d %d\n", triangle_cavity_size()+1, 4*triangle_cavity_size()+2);
    fprintf (fp_rw, "%d %d\n", 1, tri->getVertexCount());
    init_read_triangle_cavity();
    Triangle *t = next_triangle_cavity();
    while (t != NULL)
    {
        fprintf (fp_rw, "%d %d %d %d\n", 3, t->get_vertex (0), t->get_vertex (1), t->get_vertex (2));
	t = next_triangle_cavity();
    }

    fprintf (fp_rw, "CELL_TYPES %d\n", triangle_cavity_size()+1);
    fprintf (fp_rw, "%d\n", 1);
    for (int i = 0; i < triangle_cavity_size(); i++)
      fprintf (fp_rw, "%d\n", 5);

    fprintf (fp_rw, "CELL_DATA %d\n", triangle_cavity_size()+1);
    fprintf (fp_rw, "SCALARS flag int 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");
    fprintf (fp_rw, "%d\n", 0);
     init_read_triangle_cavity();
     t = next_triangle_cavity();
    while (t != NULL)
    {
        fprintf (fp_rw, "%d\n", t->get_flag ());
	t = next_triangle_cavity();
    }

    fclose(fp_rw);
    
#if defined(CIRCLE_DUMP)   
    
#if defined(Linf)

    sprintf(filename, "cavity_polygon_%d.vtk", count);
    printf(" write file --> %s\n", filename);

    fp_rw = fopen (filename, "w");

    fprintf (fp_rw, "# vtk DataFile Version 2.0\n");
    fprintf (fp_rw, "Really cool data\n");
    fprintf (fp_rw, "ASCII\n");
    fprintf (fp_rw, "DATASET POLYDATA\n");
    fprintf (fp_rw, "POINTS %d float\n", 4*triangle_cavity_size());

    init_read_triangle_cavity();
    t = next_triangle_cavity();
    while (t != NULL)
    {        
	double p0[2] = {(tri->get_center(t))[0]-(tri->get_center(t))[2], (tri->get_center(t))[1]-(tri->get_center(t))[3]};
        double p1[2] = {(tri->get_center(t))[0]+(tri->get_center(t))[2], (tri->get_center(t))[1]-(tri->get_center(t))[3]};
        double p2[2] = {(tri->get_center(t))[0]+(tri->get_center(t))[2], (tri->get_center(t))[1]+(tri->get_center(t))[3]};
        double p3[2] = {(tri->get_center(t))[0]-(tri->get_center(t))[2], (tri->get_center(t))[1]+(tri->get_center(t))[3]};
	
	double angle = t->get_angle();
        double x, y;
	x = p0[0];
	y = p0[1];
        p0[0] = x*cos(angle) - y*sin(angle);
        p0[1] = x*sin(angle) + y*cos(angle);
        x = p1[0];
        y = p1[1];
        p1[0] = x*cos(angle) - y*sin(angle);
        p1[1] = x*sin(angle) + y*cos(angle);
        x = p2[0];
        y = p2[1];
        p2[0] = x*cos(angle) - y*sin(angle);
        p2[1] = x*sin(angle) + y*cos(angle);
        x = p3[0];
        y = p3[1];
        p3[0] = x*cos(angle) - y*sin(angle);
        p3[1] = x*sin(angle) + y*cos(angle);
	
	fprintf(fp_rw, "%lf %lf %lf\n", p0[0], p0[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p1[0], p1[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p2[0], p2[1], 0.0);
        fprintf(fp_rw, "%lf %lf %lf\n", p3[0], p3[1], 0.0);
	
	t = next_triangle_cavity();
    }
    
    fprintf (fp_rw, "POLYGONS %d %d\n", triangle_cavity_size(), 5*triangle_cavity_size());
    for (int i = 0; i < triangle_cavity_size(); i++)
    {
        fprintf (fp_rw, "%d %d %d %d %d\n", 4, 4*i, 4*i+1, 4*i+2, 4*i+3);
    }
    
    fprintf (fp_rw, "CELL_DATA %d\n", triangle_cavity_size());
    fprintf (fp_rw, "SCALARS index int 1\n");
    fprintf (fp_rw, "LOOKUP_TABLE default\n");      
    for (int i = 0; i < triangle_cavity_size(); i++)
      fprintf (fp_rw, "%d\n", i);
    
    fclose(fp_rw); 
#endif
#endif 
    
    count++;
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                              Circumcenter                                           //
/////////////////////////////////////////////////////////////////////////////////////////////////////////
/**
 *    The compute circumcenter fonction.
 *    @param *q the two coordinates of the vertex.
 *    @param *&t the triangle which wil be find.
 *    @param *&i the edge of the triangle from which we start the test.
 *    @return the index of the found triangle or an error code.
 */
// Given 3 points a, b, c. Circumcenter of these three points is define by its center p and its radius R :
//
// p.x = (((a.x - c.x) * (a.x + c.x) + (a.y - c.y) * (a.y + c.y)) / 2 * (b.y - c.y)
//    -  ((b.x - c.x) * (b.x + c.x) + (b.y - c.y) * (b.y + c.y)) / 2 * (a.y - c.y))
//     / D
//
// p.y = (((b.x - c.x) * (b.x + c.x) + (b.y - c.y) * (b.y + c.y)) / 2 * (a.x - c.x)
//     -  ((a.x - c.x) * (a.x + c.x) + (a.y - c.y) * (a.y + c.y)) / 2 * (b.x - c.x))
//     / D
//
// where D = (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y)
//
// The _squared_ circumradius is then:
//
// r^2 = (c.x - p.x)^2 + (c.y - p.y)^2
//
// This approach uses Cramer's Rule to find the intersection of two
// perpendicular bisectors of triangle edges.
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                          Cercle circonscrit                                                 //
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// Given 3 points a, b, c. Circumcenter of these three points is define by its center p and its radius R :
///
/// \f$ p_x = \frac{((a_x - c_x)(a_x + c_x) + (a_y - c_y)(a_y + c_y)) * (b_y - c_y)
///     -  ((b_x - c_x)(b_x + c_x) + (b_y - c_y)(b_y + c_y)) * (a_y - c_y)}{2D} \f$
///
/// \f$ p_y = \frac{((b_x - c_x)(b_x + c_x) + (b_y - c_y)(b_y + c_y)) * (a_x - c_x)
///     -  ((a_x - c_x)(a_x + c_x) + (a_y - c_y)(a_y + c_y)) * (b_x - c_x)}{2D} \f$
///
/// where \f$ D = (a_x - c_x)(b_y - c_y) - (b_x - c_x)(a_y - c_y) \f$
///
/// The _squared_ circumradius is then:
///
/// \f$ r^2 = (c_x - p_x)^2 + (c_y - p_y)^2 \f$
///
/// This approach uses Cramer's Rule to find the intersection of two perpendicular bisectors of triangle edges.
///
/// Let :
/// \f[ \left\{\begin{matrix}
/// ax+by = e\\
/// cx+dy = f\end{matrix}\right \f]
/// So :
/// \f[ x = { \begin{vmatrix}e&b\\f&d\end{vmatrix} \over \begin{vmatrix}a&b\\c&d\end{vmatrix} } = { ed - bf \over ad - bc} \f]
/// and
///\f[ y = { \begin{vmatrix}a&e\\c&f\end{vmatrix} \over \begin{vmatrix}a&b\\c&d\end{vmatrix} } = { af - ec \over ad - bc} \f]
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void Triangulation::update_circumcircle (Triangle * t) const
{
  double circle[3];
#if defined(POWER)
    power_circumcircle (t, circle);
#else
    circumcircle (t, circle);
#endif
    t->set_circle_center (circle[0], circle[1], circle[2]);
}

void Triangulation::circumcircle (Triangle * t, double *circle) const
{
    int v0, v1, v2;
    v0 = t->get_vertex (0);
    v1 = t->get_vertex (1);
    v2 = t->get_vertex (2);
    double x0, y0, x1, y1, x2, y2;
    double *ptr;
    ptr = v[v0].p; // ptr pointe maintenant sur les coordonnees du sommet d'indice v0, premier sommet du triangle *t
    x0 = *ptr++; // on recupere la coordonnee x, puis on incremente le pointeur (attention a l'ordre des operations : *ptr++)
    y0 = *ptr;
    ptr = v[v1].p; // meme chose pour v1
    x1 = *ptr++;
    y1 = *ptr;
    ptr = v[v2].p; // meme chose pour v2
    x2 = *ptr++;
    y2 = *ptr;
    double a11, a21, a12, a22, b1, b2, d; // calcul du centre du cercle circonscrit
    d = x0 * x0 + y0 * y0;
    a11 = x1 - x0;
    a12 = y1 - y0;
    b1 = .5 * (x1 * x1 + y1 * y1 - d);
    a21 = x2 - x0;
    a22 = y2 - y0;
    b2 = .5 * (x2 * x2 + y2 * y2 - d);
    d = 1. / (a11 * a22 - a21 * a12);
    
    circle[0] = (b1 * a22 - b2 * a12) * d;
    circle[1] = (b2 * a11 - b1 * a21) * d;
    circle[2] = sqrsum (circle[0] - x0, circle[1] - y0);
//     t->set_circle_center (x, y, d);
#ifdef DEBUG

    double u0, u1, u2;
    u0 = sqrsum (x - x0, y - y0);
    u1 = sqrsum (x - x1, y - y1);
    u2 = sqrsum (x - x2, y - y2);
    std::cerr << u0 << " = " << u1 << " = " << u2 << "\n";
#endif
}

void Triangulation::power_circumcircle (Triangle * t, double *circle) const
{
    int v0, v1, v2;
    v0 = t->get_vertex (0);
    v1 = t->get_vertex (1);
    v2 = t->get_vertex (2);
    double x0, y0, x1, y1, x2, y2;
    double w0, w1, w2;
    double *ptr;
    ptr = v[v0].p; // ptr pointe maintenant sur les coordonnees du sommet d'indice v0, premier sommet du triangle *t
    x0 = *ptr++; // on recupere la coordonnee x, puis on incremente le pointeur (attention a l'ordre des operations : *ptr++)
    y0 = *ptr;
    w0 = v[v0].w;
//     w0 *= w0;
    ptr = v[v1].p; // meme chose pour v1
    x1 = *ptr++;
    y1 = *ptr;
     w1 = v[v1].w;
//      w1 *= w1;
    ptr = v[v2].p; // meme chose pour v2
    x2 = *ptr++;
    y2 = *ptr;
     w2 = v[v2].w;
//      w2 *= w2;
     
    double a11, a21, a12, a22, b1, b2, d; // calcul du centre du cercle circonscrit
    d = x0 * x0 + y0 * y0;
    a11 = x1 - x0;
    a12 = y1 - y0;
    a21 = x2 - x0;
    a22 = y2 - y0;
         // compute parameters;
    double ab2 = a11*a11+a12*a12;
    double ac2 = a21*a21+a22*a22;
    
//     printf("w0 = %lf, w1 = %lf, w2 = %lf\n", w0, w1, w2);
    
     double t0 = (ab2+w0-w1)/(2.0*ab2);
     double t1 = (ac2+w0-w2)/(2.0*ac2);  
     
//      printf("w0 = %lf, w1 = %lf, t0 = %lf (%lf %lf)\n", w0, w1, t0, t0*ab2-w0, (1.-t0)*ab2-w1);
//      printf("w0 = %lf, w2 = %lf, t1 = %lf (%lf %lf)\n", w0, w2, t1, t1*ac2-w0, (1.-t1)*ac2-w2);
    
    b1 = t0*x1*x1 + t0*y1*y1 - (1.0-t0)*x0*x0 - (1.0-t0)*y0*y0 + (1.0-2.0*t0)*x0*x1 + (1.0-2.0*t0)*y0*y1;
    b2 = t1*x2*x2 + t1*y2*y2 - (1.0-t1)*x0*x0 - (1.0-t1)*y0*y0 + (1.0-2.0*t1)*x0*x2 + (1.0-2.0*t1)*y0*y2;

    d = 1. / (a11 * a22 - a21 * a12);
    double x, y;
//     x = (b1 * a22 - b2 * a12) * d;
//     y = (b2 * a11 - b1 * a21) * d;
//     d = sqrsum (x - x0, y - y0) - w0;
    
    circle[0] = (b1 * a22 - b2 * a12) * d;
    circle[1] = (b2 * a11 - b1 * a21) * d;
    circle[2] = sqrsum (circle[0] - x0, circle[1] - y0) - w0;
    
//     assert(!isnan(x));
//     assert(!isnan(y));
//     assert(!isnan(d));
    
//     t->set_circle_center (x, y, d);

#ifdef DEBUG
    double u0, u1, u2;
    u0 = sqrsum (x - x0, y - y0) - w0;
    u1 = sqrsum (x - x1, y - y1) - w1;
    u2 = sqrsum (x - x2, y - y2) - w2;
    std::cerr << u0 << " = " << u1 << " = " << u2 << "\n";
#endif
}

void Triangulation::update_circumsquare (Triangle * t) const
{
    //     printf("--> tetra %d %d %d %d\n", t->get_vertex (0), t->get_vertex (1), t->get_vertex (2), t->get_vertex (3));
    double square[4];
#if defined(ANGLE)    
    circumsquare(t, square, _angle);
    t->set_square_center (square[0], square[1], square[2], square[3]);
//     t->set_angle(_angle, _cosa, _sina);
    double angle2 = 0.0;
    for (int i = 0; i < 3; i++)
    {
        Vertex *pv = v+t->get_vertex (i);
	angle2 += pv->a;
    }
    angle2 /= 3.0;
    t->set_angle(angle2, cos(angle2), sin(angle2));
#else
    circumsquare(t, square);
    t->set_square_center (square[0], square[1], square[2], square[3]);
    t->set_angle(_angle, _cosa, _sina);
#endif
}

#if defined(ANGLE)
void Triangulation::circumsquare (Triangle * t, double *square, double angle) const
#else
void Triangulation::circumsquare (Triangle * t, double *square) const
#endif
{
#if defined(POWER)
    double tc[3][3];
#else
    double tc[3][2];
#endif
    double *ptr = NULL;
    double minmax[2][2] = {{INFINITY, -INFINITY}, {INFINITY, -INFINITY}};
    
    double bound_values[2][2];
    int bound_vertices[2][2][3];
    int bound_nbv[2][2];
    int vertices_nb_bound[3];
    int score_enlarge_bound[2];
    
    // initialize values
    for (int i = 0; i < 2; i++)
    {
        bound_values[i][0] = INFINITY;
        bound_values[i][1] = -INFINITY;
        bound_nbv[i][0] = 0;
        bound_nbv[i][1] = 0;
        for (int j = 0; j < 3; j++)
        {
            bound_vertices[i][0][j] = -1;
            bound_vertices[i][1][j] = -1;
        }
    }
    
#if defined(POWER)
    // get vertices coord and compute box
    int bound_index[2][2] = {{-1,-1},{-1,-1}};
    for (int i = 0; i < 3; i++)
    {
        vertices_nb_bound[i] = 0;
        ptr = v[t->get_vertex (i)].p;
	tc[i][2] = v[t->get_vertex (i)].w;
        for (int j = 0; j < 2; j++)
        {
            tc[i][j] = *ptr++;
            if (tc[i][j]-tc[i][2] < bound_values[j][0])
	    {
                bound_values[j][0] = tc[i][j]-tc[i][2];
		bound_index[j][0] = i;
	    }
            if (tc[i][j]+tc[i][2] > bound_values[j][1])
	    {
                bound_values[j][1] = tc[i][j]+tc[i][2];
		bound_index[j][1] = i;
	    }
        }
    }
    
    bound_values[0][0] += 2*tc[bound_index[0][0]][2];
    bound_values[0][1] -= 2*tc[bound_index[0][1]][2];
    bound_values[1][0] += 2*tc[bound_index[1][0]][2];
    bound_values[1][1] -= 2*tc[bound_index[1][1]][2]; 

//     printf("triangle : v0(%lf, %lf) w %lf - v1(%lf, %lf) w %lf - v2(%lf, %lf) w %lf\n",
// 	   tc[0][0], tc[0][1], tc[0][2],
// 	   tc[1][0], tc[1][1], tc[1][2],
// 	   tc[2][0], tc[2][1], tc[2][2]
// 	  );
//     printf("bbox --> min = %lf %lf, max = %lf %lf\n", bound_values[0][0], bound_values[1][0], bound_values[0][1], bound_values[1][1]);

#else
    
#if defined(ANGLE)
    double angle2 = 0.0;
    for (int i = 0; i < 3; i++)
    {
        Vertex *pv = v+t->get_vertex (i);
	angle2 += pv->a;
    }
    angle2 /= 3.0;
#endif
    

// get vertices coord and compute box
    for (int i = 0; i < 3; i++)
    {
        vertices_nb_bound[i] = 0;
        Vertex *pv = v+t->get_vertex (i);
#if defined(ANGLE)	
// 	tc[i][0] = pv->p[0]*cos(angle) + pv->p[1]*sin(angle);
//         tc[i][1] = -pv->p[0]*sin(angle) + pv->p[1]*cos(angle);
	tc[i][0] = pv->p[0]*cos(angle2) + pv->p[1]*sin(angle2);
        tc[i][1] = -pv->p[0]*sin(angle2) + pv->p[1]*cos(angle2);
#else
	tc[i][0] = pv->p[0];
	tc[i][1] = pv->p[1];
#endif
        for (int j = 0; j < 2; j++)
        {
            if (tc[i][j] < bound_values[j][0])
                bound_values[j][0] = tc[i][j];
            if (tc[i][j] > bound_values[j][1])
                bound_values[j][1] = tc[i][j];  
	}
    }

//     printf("triangle : v0(%lf, %lf) v1(%lf, %lf) v2(%lf, %lf)\n", tc[0][0], tc[0][1], tc[1][0], tc[1][1], tc[2][0], tc[2][1]);
//     printf("bbox --> min = %lf %lf, max = %lf %lf\n", bound_values[0][0], bound_values[1][0], bound_values[0][1], bound_values[1][1]);
#endif

    // infinite radius if flat box and exit
    for (int i = 0; i < 2; i++)
    {
        if (bound_values[i][0] == bound_values[i][1])
        {
//             t->set_square_center (INFINITY, INFINITY, INFINITY, INFINITY);
	    square[0] = square[1] = square[2] = square[3] = INFINITY;
//             printf("--> infinite circumcube : flat bbox\n");
            return;
        }
    }

    // compute max direction
    double dl[2];
    for (int i = 0; i < 2; i++)
        dl[i] = bound_values[i][1]-bound_values[i][0];

//     printf("dlx = %lf, dly = %lf\n", dl[0], dl[1]);

    int max_dir = 0;
    if (dl[0] < dl[1])
        max_dir = 1; 

//     printf("max dir --> %d\n", max_dir);

    // return if equal
    if ((dl[0] == dl[1]))
    {
//         assert(max_dir == -1);
        t->set_square_center ((bound_values[0][0]+bound_values[0][1])/2.0, (bound_values[1][0]+bound_values[1][1])/2.0,
                            (bound_values[0][1]-bound_values[0][0])/2.0, (bound_values[1][1]-bound_values[1][0])/2.0);
        return;
    }

    assert(max_dir != -1);
    
    // fill bound counters

#if defined(POWER)
    printf("checking ...\n");
    for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 2; j++)
        {
	    printf("testing v %d %lf --> %lf\n", i, tc[i][j]+tc[i][2], bound_values[j][0]);
	    printf("testing v %d %lf --> %lf\n", i, tc[i][j]-tc[i][2], bound_values[j][1]);
//             if (tc[i][j]+tc[i][2] == bound_values[j][0])
	    if (fabs(tc[i][j]+tc[i][2] - bound_values[j][0]) < 1e-12)
            {
		printf("--> accepted\n");
                vertices_nb_bound[i]++;
		bound_nbv[j][0]++;
                bound_vertices[j][0][i] = i;
            }
//             else if (tc[i][j]-tc[i][2] == bound_values[j][1])
	    else if (fabs(tc[i][j]-tc[i][2] - bound_values[j][1]) < 1e-12) 
            {
	        printf("--> accepted\n");
                vertices_nb_bound[i]++;
		bound_nbv[j][1]++;
                bound_vertices[j][1][i] = i;
            }
        }

        if (vertices_nb_bound[i] == 0)
        {
//             t->set_square_center (INFINITY, INFINITY, INFINITY, INFINITY);
	    square[0] = square[1] = square[2] = square[3] = INFINITY;
// 	    printf("--> infinite circumcube : no vertex bound\n");
            return;
        }
    }

#else

for (int i = 0; i < 3; i++)
    {
        for (int j = 0; j < 2; j++)
        {
            if (tc[i][j] == bound_values[j][0])
            {
                vertices_nb_bound[i]++;
		bound_nbv[j][0]++;
                bound_vertices[j][0][i] = i;
            }
            else if (tc[i][j] == bound_values[j][1])
            {
                vertices_nb_bound[i]++;
		bound_nbv[j][1]++;
                bound_vertices[j][1][i] = i;
            }
        }

        if (vertices_nb_bound[i] == 0)
        {
            square[0] = square[1] = square[2] = square[3] = INFINITY;
            return;
        }
    }
#endif
    
#if defined(DEBUG)
    printf("## bounds : ");
    for (int i = 0; i < 2; i++)
    {
        printf("bound %d -->\n", i);
        printf("dir - : %d vertices --> ", bound_nbv[i][0]);
        for (int j = 0; j < 3; j++)
            if (bound_vertices[i][0][j] != -1)
                printf(" %d", bound_vertices[i][0][j]);
        printf("\n");
        printf("dir + : %d vertices --> ", bound_nbv[i][1]);
        for (int j = 0; j < 3; j++)
            if (bound_vertices[i][1][j] != -1)
                printf(" %d", bound_vertices[i][1][j]);
        printf("\n");
    }
#endif
    
    int extdir = (max_dir+1)%2;

    // choose extension
    int extent[2] = {1, 1};
    for (int i = 0; i < 3; i++)
    {
        int vx = bound_vertices[extdir][0][i];
        if (vx != -1)
        {
            if (vertices_nb_bound[vx] == 1)
                extent[0] = 0;
        }
    }

    for (int i = 0; i < 3; i++)
    {
        int vx = bound_vertices[extdir][1][i];
        if (vx != -1)
        {
            if (vertices_nb_bound[vx] == 1)
                extent[1] = 0;
        }
    }

    // cleanup structures
    if (extent[0])
    {
        for (int i = 0; i < 3; i++)
        {
// 		assert(bound_nbv[cdir][0] >= 2);
            int vx = bound_vertices[extdir][0][i];
            if (vx != -1)
            {
                bound_vertices[extdir][0][i] = -1;
                vertices_nb_bound[vx]--;
                bound_nbv[extdir][0]--;
            }
        }
    }

    if (extent[1])
    {
        for (int i = 0; i < 3; i++)
        {
//                 assert(bound_nbv[cdir][1] >= 2);
            int vx = bound_vertices[extdir][1][i];
            if (vx != -1)
            {
                bound_vertices[extdir][1][i] = -1;
                vertices_nb_bound[vx]--;
                bound_nbv[extdir][1]--;
            }
        }
    }

//     printf("extensions --> %d %d\n", extent[0], extent[1]);

    double enlarge = dl[max_dir]-dl[extdir];
    double factor = 0.5;
    
//         printf("--> enlarging dir %d of %lf\n", cdir, enlarge);
//         printf("dir %d --> min = %lf, max = %lf\n", cdir, bound_values[cdir][0], bound_values[cdir][1]);

    if (extent[0] && extent[1])
    {
        bound_values[extdir][0] -= enlarge*factor;
        bound_values[extdir][1] += enlarge*factor;
    }
    else if (extent[0])
    {
        bound_values[extdir][0] -= enlarge;
    }
    else if (extent[1])
    {
        bound_values[extdir][1] += enlarge;
    }
    else
    {
	square[0] = square[1] = square[2] = square[3] = INFINITY;
        return;
    }
//     printf(" --> min = %lf, max = %lf\n", bound_values[extdir][0], bound_values[extdir][1]);

    square[0] = (bound_values[0][0]+bound_values[0][1])/2.0;
    square[1] = (bound_values[1][0]+bound_values[1][1])/2.0;
    
#if defined(POWER)
    double rx, ry;
    
//     if (fabs(tc[0][0]-x) >= fabs(tc[0][1]-y))
//       rx = ry = fabs(tc[0][0]-x)-tc[0][2];
//     else
//       rx = ry = fabs(tc[0][1]-y)-tc[0][2];
    
    double pmin = tc[0][2];
    int imin = 0;
    if (pmin > tc[1][2])
    {
      pmin > tc[1][2];
      imin = 1;
    }
    if (pmin > tc[2][2])
    {
      pmin > tc[2][2];
      imin = 2;
    }
        
    square[2] = square[3] = std::max(fabs(square[0]-tc[imin][0])-pmin, fabs(square[1]-tc[imin][1])-pmin);
    
    // control
      double lc[3];
      printf("--> new triangle\n");
    for (int i = 0; i < 3; i++)
    {
	double ll[2];
	ll[0] = fabs(square[0]-tc[i][0]);
	ll[1] = fabs(square[1]-tc[i][1]);
        double l = std::max(ll[0], ll[1]);
        lc[i] = l - tc[i][2];
	printf("--> P = %lf\n", lc[i]);
    }
//     for (int i = 0; i < 3; i++)
//     {
//       assert(fabs(lc[i] - lc[(i+1)%3]) < 1e-8);
//     }
#else
    square[2] = (bound_values[0][1]-bound_values[0][0])/2.0;
    square[3] = (bound_values[1][1]-bound_values[1][0])/2.0;
    assert(square[2] > 0);
    assert(square[3] > 0);
#endif

//     printf("square --> xmin = %lf, xmax = %lf (sum %lf), ymin = %lf, ymax = %lf (sum %lf)\n", bound_values[0][0], bound_values[0][1], rx,
// 	    bound_values[1][0], bound_values[1][1], ry);
    
//     t->set_square_center (x, y, rx, ry);
}

// static long onvertex, onedge;

/////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                              Locate                                                 //
/////////////////////////////////////////////////////////////////////////////////////////////////////////
/**
 *    The locate fonction.
 *    @param *q the two coordinates of the vertex.
 *    @param *&t the triangle which will be find.
 *    @param *&i the edge of the triangle from which we start the test.
 *    @return the index of the found triangle or an error code.
 */
/// Given a point |p| with coordinates |(q[0],q[1])|, |locate(q,t,i)| locates the triangle that encloses |p|, searching from
/// the given edge (t,i), and returns the vertex index p. Since all sites are enclosed by an epsilon-enlarged bounding box,
/// |locate| will never fail. When return, it sets |t| by the triangle. It returns 0 if |p| coincides with an existing vertex
/// in which case no new vertex is created; it returns |-p| if |p| is on edge |i| of |t|;
/// it returns |p| if |p| is strictly inside |t|.
/// To locate the triangle, we continue to test if |q| is on the left of or on every edge of |t|, and try to find the triangle
/// of which each edge passes the test;
/// if one edge fails the test, we immediately change |t| to its adjacent triangle and start over the tests.
////////////////////////////////////////////////////////////////////////////////////////////////////////

int Triangulation::locate (double *q, Triangle * &t, int &i)
{
    /**
      Il faut declarer un compteur.
      A chaque fois qu'on fait get_adjacent on incremente le compteur.
      On veut stocker le max et la moyenne de ce compteur.
    */
    int ni, pi;
    double fi, fni, fpi;
    Triangle *s;
    if (q[0] == v[get_vertex (t, 0)].p[0] && q[1] == v[get_vertex (t, 0)].p[1])
        return 0; // cas ou le sommet insere est confondu avec un sommet existant
    if (q[0] == v[get_vertex (t, 1)].p[0] && q[1] == v[get_vertex (t, 1)].p[1])
        return 0; // cas ou le sommet insere est confondu avec un sommet existant
    if (q[0] == v[get_vertex (t, 2)].p[0] && q[1] == v[get_vertex (t, 2)].p[1])
        return 0; // cas ou le sommet insere est confondu avec un sommet existant

    fi = lefton_test (q, t, i);// The |lefton_test| is a counterclockwise |ccw| test.
	
    
    if (fi < 0) // on passe dans le triangle voisin (ieme adjacence)
    {
        s = get_adjacent (t, i);
        i = get_index (t, i);
        t = s;
        fi = -fi;
    }

    int stop = 0;
    while (!stop)
    {
	int choice = rand()%2;
        switch (choice)
        {
        case 0:
            ni = Triangle::next (i); // retourne le ieme+1 vertex (c'est a dire 0, 1 ou 2)
            if ((fni = lefton_test (q, t, ni)) < 0.0) // on passe dans le triangle voisin (ieme+1 adjacence)
            {
                fi = -fni;
                i = get_index (t, ni);
                t = get_adjacent (t, ni);
            }
            else
            {
                pi = Triangle::prev (i); // retourne le ieme-1 vertex (c'est a dire 0, 1 ou 2)
                if ((fpi = lefton_test (q, t, pi)) < 0.0) // on passe dans le triangle voisin (ieme-1 adjacence)
                {
                    fi = -fpi;
                    i = get_index (t, pi);
                    t = get_adjacent (t, pi);
                }
                else
                    stop = 1; // on est arrive dans le bon triangle
            }
            break;
        case 1 :
            pi = Triangle::prev (i); // retourne le ieme-1 vertex (c'est a dire 0, 1 ou 2)
            if ((fpi = lefton_test (q, t, pi)) < 0.0) // on passe dans le triangle voisin (ieme-1 adjacence)
            {
                fi = -fpi;
                i = get_index (t, pi);
                t = get_adjacent (t, pi);
            }
            else
            {
                ni = Triangle::next (i); // retourne le ieme+1 vertex (c'est a dire 0, 1 ou 2)
                if ((fni = lefton_test (q, t, ni)) < 0.0) // on passe dans le triangle voisin (ieme+1 adjacence)
                {
                    fi = -fni;
                    i = get_index (t, ni);
                    t = get_adjacent (t, ni);
                }
                else
                    stop = 1; // on est arrive dans le bon triangle
            }
            break;
        }
    }

    last_tri = t; // on met a jour ces deux variables globales
    last_i = i;
            
    if (fni == 0.0 && fpi == 0.0)
        return 0;

    // test cas particuliers (voir les valeurs de retour dans Triangulation::insert)
    if (fi == 0.0)
        return -1;
    
    if (fni == 0.0)
    {
        i = ni;
        return -1;
    }

    if (fpi == 0.0)
    {
        i = pi;
        return -1;
    }
    
    return 1;
}

long ccwcount;

void Triangulation::print (const char *file) const
{
    char node[20], ele[20];
    std::cerr << "number of calls of ccw: " << ccwcount << "\n";
    std::cerr << "number of on vertex: " << onvertex << "\n";
    std::cerr << "number of on edge: " << onedge << "\n";
    std::cerr << "number of on lefton test: " << leftoncount << "\n";


    sprintf (node, "%s.node", file);
    sprintf (ele, "%s.ele", file);
    FILE *fnode = fopen (node, "w");
    FILE *fele = fopen (ele, "w");

    int i, j;

    fprintf (fnode, "%d 2 %d 0\n", nv, 0);
    for (i = 0, j = 0; i < nv; ++i)
    {
        fprintf (fnode, "%d\t%f  %f\n", i + 1, v[i].p[0], v[i].p[1]);
    }

    fprintf (fele, "%d 3 0\n", ntri);
    Triangle *t, *tend = tri + ntri;
    for (i = 0, t = tri; t < tend; ++t, ++i)
    {
        fprintf (fele, "%d\t%d\t%d\t%d\n", i + 1, get_vertex (t, 0) + 1,
                 get_vertex (t, 1) + 1, get_vertex (t, 2) + 1);
    }
    fclose (fnode);
    fclose (fele);
}

std::ostream & Triangulation::print_vertex (std::ostream & os, int i) const
{
    const double *
    p = v[i].p;
    os << p[0] << " " << p[1];
    return os;
}

std::ostream & Triangulation::print_triangle (std::ostream & os, Triangle * t,
        bool shift) const
{
    if (shift)
    {
        os << get_vertex (t, 0) + 1 << " "
        << get_vertex (t, 1) + 1 << " " << get_vertex (t, 2) + 1;
    }
    else
    {
        os << get_vertex (t, 0) << " "
        << get_vertex (t, 1) << " " << get_vertex (t, 2);
    }
    return os;
}

static const double THRESHOLD = 2.98023232758737586699e-8;

bool Triangulation::is_valid () const
{
    bool result = true;
    int i;
    Triangle *t, *tend = tri + ntri;
    for (t = tri; t < tend; ++t) // test tous les triangles ...
    {
        for (i = 0; i < nv; ++i) // pour chacun des sommets
        {
//             double d = power_distance (t, i);
	    double d = max_distance (t, i);
            if (d < -THRESHOLD)
            {
                result = false;
                std::cerr << "triangulation error: vertex " << i << " = (";
                print_vertex (std::cerr, i) << ") is in triangle (";
                print_triangle (std::cerr, t) << ")\n";
//                 std::cerr << "\tpower distance = " << d << "\n";
		std::cerr << "\tmax distance = " << d << "\n";
            }
        }
    }
    return result;
}

/// |check_nb() checks if the neighboring relationships are well maintained.

bool Triangulation::check_nb () const
{
    bool result = true;
    int i, j, i1;
    Triangle *t, *tend = tri + ntri, *t1, *s;
    for (t = tri; t < tend; ++t) // test tous les triangles ...
    {
        for (i = 0; i < 3; ++i) // pour chacun de ses trois voisins
        {
            get_nb (t, i, s, j);
            get_nb (s, j, t1, i1);
            if (t != t1 || i != i1)
            {
                std::cerr << "neighbor error: triangle (";
                print_triangle (std::cerr, t) << ") " << i << "-neighbor is (";
                print_triangle (std::cerr,
                                s) << ") with " << j << "-neighbor\n";
                std::cerr << "But, its neighbor is triangle (";
                print_triangle (std::cerr, t) << ") " << i1 << "-neighbor\n";
                result = false;
            }
        }
    }
    std::cerr << "check neighbor = " << std::boolalpha << result << "\n";
    return result;
}

// The implementation of |find| is obvious but tedious.
// We search the nonempty bucket level by level starting with level 0 at bucket |b|.
// At each level, the buckets form a round strip; and we start to search from the lower left
// bucket, then we continue along the strip counterclockwisely, i.e., move to the right, move
// upwardly, move to the left, move downwardly. Before searching at the current level, the
// bucket |b| is at the lower left corner for the last level; after searching at current
// level, |b| automatically comes back to the lower left corner.

Vertex* Triangulation::find (Vertex ** b, int x, int y) const
{
    if (*b != NULL)
        return *b;
    int xmin, xmax, ymin, ymax;
    xmin = xmax = x;
    ymin = ymax = y;
    for (int level = 0; level < SEARCH_LEVEL; ++level) // boucle jusqu'au niveau de recherche maximal
    {
        if (xmin > 0) // determine xmin
        {
            --xmin; // le decremente si possible
            --b; //
        }
        if (xmax < BUCKET_SIZE - 1)
        {
            ++xmax;
        }
        if (ymin > 0)
        {
            --ymin;
            b -= BUCKET_SIZE;
        }
        if (ymax < BUCKET_SIZE - 1)
        {
            ++ymax;
        }
        for (x = xmin; x < xmax; ++x, ++b)
            if (*b != NULL)
                return *b;
        for (y = ymin; y < ymax; ++y, b += BUCKET_SIZE)
            if (*b != NULL)
                return *b;
        for (; x > xmin; --x, --b)
            if (*b != NULL)
                return *b;
        for (; y > ymin; --y, b -= BUCKET_SIZE)
            if (*b != NULL)
                return *b;
    }
    return NULL;
}

// void Triangulation::update_circumsquare (Triangle * t) const
// {
//     int v0, v1, v2;
//     v0 = t->get_vertex (0);
//     v1 = t->get_vertex (1);
//     v2 = t->get_vertex (2);
//     double x0, y0, x1, y1, x2, y2;
//     double *ptr;
//     ptr = v[v0].p; // ptr pointe maintenant sur les coordonnees du sommet d'indice v0, premier sommet du triangle *t
//     x0 = *ptr++; // on recupere la coordonnee x, puis on incremente le pointeur (attention a l'ordre des operations : *ptr++)
//     y0 = *ptr;
//     ptr = v[v1].p; // meme chose pour v1
//     x1 = *ptr++;
//     y1 = *ptr;
//     ptr = v[v2].p; // meme chose pour v2
//     x2 = *ptr++;
//     y2 = *ptr;
//         
//     double xmin, xmax, ymin, ymax;
//     int ixmin[2], ixmax[2], iymin[2], iymax[2];
//     int c[3] = {0, 0, 0};
//     
//     if ((x1 == x0 && x2 == x0) || (y1 == y0 && y2 == y0))
//     {
//       t->set_square_center (INFINITY, INFINITY, INFINITY, INFINITY);
// //       printf("--> infinite circumsquare\n");
//       return;
//     }
//     
//     xmin = x0; ixmin[0] = 0; ixmin[1] = -1;
//     if (x1 == xmin) {ixmin[1] = 1;}
//     else if (x1 < xmin) { xmin = x1; ixmin[0] = 1; ixmin[1] = -1; }
//     if (x2 == xmin) {ixmin[1] = 2;}
//     else if (x2 < xmin) { xmin = x2; ixmin[0] = 2; ixmin[1] = -1;}
//         
//     xmax = x0; ixmax[0] = 0; ixmax[1] = -1;
//     if (x1 == xmax) { ixmax[1] = 1;}
//     else if (x1 > xmax) { xmax = x1; ixmax[0] = 1; ixmax[1] = -1; }
//     if (x2 == xmax) { ixmax[1] = 2;}
//     else if (x2 > xmax) { xmax = x2; ixmax[0] = 2; ixmax[1] = -1; }
//         
//     ymin = y0; iymin[0] = 0; iymin[1] = -1;
//     if (y1 == ymin) { iymin[1] = 1;}
//     else if (y1 < ymin) { ymin = y1; iymin[0] = 1; iymin[1] = -1; }
//     if (y2 == ymin) { iymin[1] = 2;}
//     else if (y2 < ymin) { ymin = y2; iymin[0] = 2; iymin[1] = -1; }
//         
//     ymax = y0; iymax[0] = 0; iymax[1] = -1;
//     if (y1 == ymax) { iymax[1] = 1;}
//     else if (y1 > ymax) { ymax = y1; iymax[0] = 1; iymax[1] = -1; }    
//     if (y2 == ymax) { iymax[1] = 2;}
//     else if (y2 > ymax) { ymax = y2; iymax[0] = 2; iymax[1] = -1; }
//     
//     c[ixmin[0]]++;
//     c[ixmax[0]]++;
//     c[iymin[0]]++;
//     c[iymax[0]]++;
//     
//     if (ixmin[1] != -1) c[ixmin[1]]++;
//     if (ixmax[1] != -1) c[ixmax[1]]++;
//     if (iymin[1] != -1) c[iymin[1]]++;
//     if (iymax[1] != -1) c[iymax[1]]++;
//         
//     if (c[0] == 0 || c[1] == 0 || c[2] == 0 )
//     {
//       t->set_square_center (INFINITY, INFINITY, INFINITY, INFINITY);
// //       printf("--> infinite circumsquare\n");
//       return;
//     }
//     
// #ifdef DEBUG 
//     printf("tri : v0(%lf, %lf) v1(%lf, %lf) v2(%lf, %lf)\n", x0, y0, x1, y1, x2, y2);
//     printf("bbox --> xmin = %lf, xmax = %lf, ymin = %lf, ymax = %lf\n", xmin, xmax, ymin, ymax);
//     printf("ixmin --> %d %d, ixmax --> %d %d, iymin --> %d %d, iymax --> %d %d\n", ixmin[0], ixmin[1], ixmax[0], ixmax[1], iymin[0], iymin[1], iymax[0], iymax[1]);
//     printf("summary : v0 --> %d, v1 --> %d, v2 --> %d\n", c[0], c[1], c[2]);
// #endif    
//     
//     int cxmin = 1, cxmax = 1, cymin = 1, cymax = 1;
//     if ( c[ixmin[0]] == 1)
//         cxmin = 0;
//     if ( c[ixmax[0]] == 1)
//         cxmax = 0;
//     if ( c[iymin[0]] == 1)
//         cymin = 0;
//     if ( c[iymax[0]] == 1)
//         cymax = 0;
//         
//     if ( ixmin[1] != -1 && c[ixmin[1]] == 1)
//         cxmin = 0;
//     if ( ixmax[1] != -1 && c[ixmax[1]] == 1)
//         cxmax = 0;
//     if ( iymin[1] != -1 && c[iymin[1]] == 1)
//         cymin = 0;
//     if ( iymax[1] != -1 && c[iymax[1]] == 1)
//         cymax = 0;
// 
//   
// #ifdef DEBUG    
//      printf("cxmin --> %d, cxmax --> %d, cymin --> %d, cymax --> %d\n", cxmin, cxmax, cymin, cymax);
// #endif 
//      
//      double rx, ry;     
//      double factor = 0.5;
// 
//     if (xmax-xmin > ymax-ymin)
//     {      
// 	double enlarge = xmax-xmin-(ymax-ymin);
//         if (cymin && cymax)
//         {
//             ymin -= enlarge*factor;
//             ymax += enlarge*factor;
//         }        
// 	else if (cymin)	  
// 	  ymin -= enlarge;
// 	else if (cymax)
// 	  ymax += enlarge;
//     }
//     else if (xmax-xmin < ymax-ymin)
//     {
// //         sql = ymax-ymin;
// 	double enlarge = ymax-ymin-(xmax-xmin);
// 	if (cxmin && cxmax)
// 	{
// 	  xmin -= enlarge*factor;
// 	  xmax += enlarge*factor;
// 	}
// 	else if (cxmin)
// 	  xmin -= enlarge;
// 	else if (cxmax)
// 	  xmax += enlarge;
//     }      
//     
// #ifdef DEBUG     
//     printf("square --> xmin = %lf, xmax = %lf (sum %lf), ymin = %lf, ymax = %lf (sum %lf)\n", xmin, xmax, xmax-xmin, ymin, ymax, ymax-ymin);
// #endif      
//     
//     rx = (xmax-xmin)/2.0;
//     ry = (ymax-ymin)/2.0;
//     
//     assert(rx > 0);
//     assert(ry > 0);
//     
//     double x = (xmin+xmax)/2.0;
//     double y = (ymin+ymax)/2.0;
//     t->set_square_center (x, y, rx, ry);
// //     t->set_square_center (x, y, std::min(rx, ry), std::min(rx, ry));
// 
// // #ifdef DEBUG
// // 
// //     double u0, u1, u2;
// //     u0 = sqrsum (x - x0, y - y0);
// //     u1 = sqrsum (x - x1, y - y1);
// //     u2 = sqrsum (x - x2, y - y2);
// //     std::cerr << u0 << " = " << u1 << " = " << u2 << "\n";
// // #endif
// }

// #if defined(Linf)
// int Triangulation::bissector(Vertex *tv[3])
// {    
//     
//     double halfsqrt2 = sqrt(2.0)/2.0;
//     
//     // pp[i][0] * x + pp[i][1] * y + pp[i][2] = 0
//     double pp[6][3];
// 
//     for (int i = 0; i < 2; i++)
//     {
//         Vertex vv[2];
//         double dx = tv[i+1]->p[0]-tv[i]->p[0];
//         double dy = tv[i+1]->p[1]-tv[i]->p[1];
// 	
// 	int ii = 3*i;
// 	
//         if (fabs(dx) == fabs(dy))
//         {
//             // cas particulier
//         }
//         else
//         {
//             double rsquare1, rsquare2, dd1, dd2;
//             double t = 0.0;
//             if (fabs(dx) >= fabs(dy))
//             {
//                 t = (tv[i]->w-tv[i+1]->w)/(2*dx) + 0.5;
//             }
//             else
//             {
//                 t = (tv[i]->w-tv[i+1]->w)/(2*dy) + 0.5;
//             }
// 
//             vv[0]->p[0] = vv[1]->p[0] = tv[i]->p[0]+t*(tv[i+1]->p[0]-tv[i]->p[0]);
//             vv[0]->p[1] = vv[1]->p[1] = tv[i]->p[1]+t*(tv[i+1]->p[1]-tv[i]->p[1]);
// 
// // 	printf("t --> %lf\n", t);
// 
//             if (fabs(dx) >= fabs(dy))
//             {
//                 rsquare1 = t*fabs(dx);
//                 rsquare2 = (1.0-t)*fabs(dx);
// 
//                 double ddp = std::min(tv[i]->p[1]+rsquare1, tv[i+1]->p[1]+rsquare2);
//                 double ddm = std::max(tv[i]->p[1]-rsquare1, tv[i+1]->p[1]-rsquare2);
// 
//                 if (tv[i]->p[0] < tv[i+1]->p[0])
//                 {
//                     vv[0]->p[1] = ddp;
//                     vv[1]->p[1] = ddm;
//                 }
//                 else
//                 {
//                     vv[0]->p[1] = ddm;
//                     vv[1]->p[1] = ddp;
//                 }
// 
//                 if (tv[i]->p[0] < tv[i+1]->p[0])
//                 {
//                     //
//                     if (tv[i]->p[1]+rsquare1 < tv[i+1]->p[1]+rsquare2)
//                         pp[ii][0] = halfsqrt2;
//                     else
//                         pp[ii][0] = -halfsqrt2;
// 
//                     pp[ii][1] = halfsqrt2;
//                     pp[ii][2] = -pp[ii][0]*vv[0]->p[0] - pp[ii][1]*vv[0]->p[1];
// 
//                     //
//                     if (tv[i]->p[1]-rsquare1 < tv[i+1]->p[1]-rsquare2)
//                         pp[ii+2][0] = halfsqrt2;
//                     else
//                         pp[ii+2][0] = -halfsqrt2;
// 
//                     pp[ii+2][1] = halfsqrt2;
//                     pp[ii+2][2] = -pp[ii+2][0]*vv[1]->p[0] - pp[ii+2][1]*vv[1]->p[1];
//                 }
//                 else
//                 {
//                     if (tv[i]->p[1]-rsquare1 < tv[i+1]->p[1]-rsquare2)
//                         pp[ii][0] = -halfsqrt2;
//                     else
//                         pp[ii][0] = halfsqrt2;
// 
//                     pp[ii][1] = halfsqrt2;
//                     pp[ii][2] = -pp[ii][0]*vv[0]->p[0] - pp[ii][1]*vv[0]->p[1];
// 
//                     //
//                     if (tv[i]->p[1]+rsquare1 < tv[i+1]->p[1]+rsquare2)
//                         pp[ii+2][0] = -halfsqrt2;
//                     else
//                         pp[ii+2][0] = halfsqrt2;
// 
//                     pp[ii+2][1] = halfsqrt2;
//                     pp[ii+2][2] = -pp[ii+2][0]*vv[1]->p[0] - pp[ii+2][1]*vv[1]->p[1];
//                 }
// 
//                 pp[ii+1][0] = vv[0]->p[0];
//                 pp[ii+1][1] = 0.0;
//                 pp[ii+1][2] = 0.0;
//             }
//             else
//             {
//                 rsquare1 = t*fabs(dy);
//                 rsquare2 = (1.0-t)*fabs(dy);
// 
//                 double ddp = std::min(tv[i]->p[0]+rsquare1, tv[i+1]->p[0]+rsquare2);
//                 double ddm = std::max(tv[i]->p[0]-rsquare1, tv[i+1]->p[0]-rsquare2);
// 
//                 if (tv[i]->p[1] < tv[i+1]->p[1])
//                 {
//                     vv[0]->p[0] = ddm;
//                     vv[1]->p[0] = ddp;
//                 }
//                 else
//                 {
//                     vv[0]->p[0] = ddp;
//                     vv[1]->p[0] = ddm;
//                 }
// 
//                 if (tv[i]->p[1] < tv[i+1]->p[1])
//                 {
//                     //
//                     if (tv[i]->p[0]-rsquare1 < tv[i+1]->p[0]-rsquare2)
//                         pp[ii][0] = halfsqrt2;
//                     else
//                         pp[ii][0] = -halfsqrt2;
// 
//                     pp[ii][1] = halfsqrt2;
//                     pp[ii][2] = -pp[ii][0]*vv[0]->p[0] - pp[ii][1]*vv[0]->p[1];
// 
//                     //
//                     if (tv[i]->p[0]+rsquare1 < tv[i+1]->p[0]+rsquare2)
//                         pp[ii+2][0] = -halfsqrt2;
//                     else
//                         pp[ii+2][0] = halfsqrt2;
// 
//                     pp[ii+2][1] = halfsqrt2;
//                     pp[ii+2][2] = -pp[ii+2][0]*vv[1]->p[0] - pp[ii+2][1]*vv[1]->p[1];
//                 }
//                 else
//                 {
//                     if (tv[i]->p[0]+rsquare1 < tv[i+1]->p[0]+rsquare2)
//                         pp[ii][0] = -halfsqrt2;
//                     else
//                         pp[ii][0] = halfsqrt2;
// 
//                     pp[ii][1] = halfsqrt2;
//                     pp[ii][2] = -pp[ii][0]*vv[0]->p[0] - pp[ii][1]*vv[0]->p[1];
// 
//                     //
//                     if (tv[i]->p[0]-rsquare1 < tv[i+1]->p[0]-rsquare2)
//                         pp[ii+2][0] = -halfsqrt2;
//                     else
//                         pp[ii+2][0] = halfsqrt2;
// 
//                     pp[ii+2][1] = halfsqrt2;
//                     pp[ii+2][2] = -pp[ii+2][0]*vv[1]->p[0] - pp[ii+2][1]*vv[1]->p[1];
//                 }
// 
//                 pp[ii+1][0] = 0.0;
//                 pp[ii+1][1] = vv[0]->p[1];
//                 pp[ii+1][2] = 0.0;
//             }
//         }
//     }
// }
// #endif