#include <cassert>
#include <iostream>
#include <vector>
#include "geom.h"
#include "background.h"
#include "ANN/ANN.h"

#ifndef _tri_include
#define _tri_include

// #define DEBUG
// #define NDEBUG

// #define CHECK
// #define STEP_DEBUG
// #define ALIGN

// #define LLOYD
#define CLASSIC

#define Linf
// #define L2

#define ANGLE
#define DENSITY_CENTROID
// #define POWER
#define CIRCLE_DUMP

#define QUAD

// #define LAWSON_ALGORITHM
#define BOWYERWATSON_ALGORITHM

#define BUCKET_SIZE 100
#define SEARCH_LEVEL 10

#define CAVITY_CAPACITY 10000
#define FREE_CAPACITY 10000

#define ACTIVE_FACET 0
#define INACTIVE_FACET 1

#define DEFAULT_TRIANGLE_FLAG 0
#define LOCKED_TRIANGLE_FLAG  -1
#define FREED_TRIANGLE_FLAG  -2 


#if defined(DEBUG) 
#define dprintf(...) printf(__VA_ARGS__)
#else
#define dprintf(...) 
#endif 

class Vertex;
class Triangle;
class Cavity;
class Triangulation;
class Facet;

static long onvertex, onedge;

struct Vertex
{
  double p[2];
  double h;
  double w;
  double a;
  Vertex*n;

#if 1
  Triangle*t;
#endif

  

  int flag;
};

class Triangle
{
  static const int p[3];
  static const int n[3];
  int v[3];

  Triangle* adj[3];
  int ind[3];
  
  unsigned int flag;

  int pos;

    double center[4];   
    double a;
    double ca;
    double sa;

public:
  int pos2;
  int vertex(int i)const
  {
    return v[i];
  }
  static int prev(int i)
  {
    return p[i];
  }
  static int next(int i)
  {
    return n[i];
  }

  int org(int i) const
  {
    return v[next(i)];
  }

  int dest(int i)const
  {
    return v[prev(i)];
  }

  void set_vertex(int j,int v0)
  {
    v[j]= v0;
  }

  void set_adjacent(int j,Triangle*a)
  {
    adj[j]= a;
  }

  void set_opp(int j,int i)
  {
    ind[j]= i;
  }

  void set_members(int v0, int v1, int v2, Triangle *j0, Triangle *j1, Triangle *j2, int i0, int i1, int i2)
  {
    v[0]= v0;
    v[1]= v1;
    v[2]= v2;
    adj[0]= j0;
    adj[1]= j1;
    adj[2]= j2;
    ind[0]= i0;
    ind[1]= i1;
    ind[2]= i2;
  }

  int get_vertex(int j) const
  {
    return v[j];
  }

  Triangle* get_adjacent(int j) const
  {
    return adj[j];
  }

  int get_opp(int j)const
  {
    return ind[j];
  }

  bool is_modified()
  {
    return (flag & 1U) != 0;
  }

  void set_modified()
  {
    flag |= 1U;
  }

  void xor_modified()
  {
    flag ^= 1U;
  }

  const double*get_center()const
  {
    return&center[0];
  }
  
  void set_circle_center(double x,double y,double p)
  {
    center[0]= x;
    center[1]= y;
    center[2]= p;
  }
  void set_square_center(double x,double y, double rx, double ry)
  {
//     printf("center --> %lf %lf, radius = %lf\n", x, y, rx);
    center[0]= x;
    center[1]= y;
    center[2]= rx;
    center[3]= ry;
  }
  
  void set_angle(double angle, double cosa, double sina)
  {
    a = angle;    
    ca = cosa;
    sa = sina;
  }
  
  double get_angle()
  {
    return a;
  }

  void set_flag(int flag)
  {
  	 pos = flag;
  }

  int get_flag()
  {
  	 return pos;
  }  

};


class Triangulation
{
public:
  Triangulation();
  ~Triangulation();

  void init(double*s,int n,double epsilon,bool regular= false);
  void check_power(double *vx, int size);
  void export_vtk();
  void debug_step(Triangle *t, Triangle *s);
  
  void set_background(Background *bg)
  {
    _bg = bg;
  }
  
int check_delaunay();

  int ptind(int j)const
  {
    return j<<1;
  }

  bool insert_lawson(Vertex*pv);  
  bool check_previous(Triangle * t, int i, int &stack);
  bool check_next(Triangle * t, int i, int &stack);
  bool check_neigh (Triangle * t, int i/*, int p*/, int &stack);
  
  bool insert_bowyer_watson(Vertex*pv);
    bool build_initial_cavity(Cavity *cav, int p);
  bool build_final_cavity(Cavity *cav, int p);
  void ball(Vertex *v, int &nb, Triangle **tlist);
  bool correct_cavity(Cavity *cav, int p);
  void finalize(Cavity *cav);
  bool star_cavity(Cavity *cav, int p);
  int check_adjacencies();
  
  void align (double *vx, int size, double &xmin, double &xmax, double &ymin, double &ymax);
  void disalign();
  void disalign(double *p);

  Triangle* copy(Triangle *t)
  {
    Triangle *s= new_triangle();
    *s= *t;
    return s;
  }

  void change_id(Triangle *t, int i)
  {
    Triangle *s;
    int j;
    get_nb(t,i,s,j); // recupere les infos du ieme triangle voisin
    if(s==NULL)
      return;
    set_adjacent(s,j,t); // si s existe, marque t comme le jieme triangle adjacent de s
  }

  int getVertexCount() const
  {
    return nv;
  }

  int getTriangleCount() const
  {
    return ntri;
  }

  void circumcircle (Triangle * t, double *circle) const;
  void power_circumcircle (Triangle * t, double *circle) const;
  void update_circumcircle(Triangle*pt) const;  
#if defined(ANGLE)
    void circumsquare (Triangle * t, double *square, double angle) const;
#else
  void circumsquare (Triangle * t, double *square) const;
#endif
  void angled_circumsquare (Triangle * t, double *square, double angle) const;
  void update_circumsquare(Triangle*t) const;

  
  bool mutual_insquare(Triangle *t, int i, int p);
  bool mutual_insquare(Triangle *t, int i, Triangle*s, int j);
  bool flip(Triangle*t, int i, Triangle*s, int j, int&stack);
  bool checkflip2to2(Triangle*t, int i, Triangle*s, int j);
#if defined(ANGLE)
  bool checkflip2to2 (Triangle * t, int i, int p, double angle);
#else
  bool checkflip2to2(Triangle*t, int i, int p);
#endif
  void flip2to2(Triangle*t, int i, Triangle*s, int j);
  
#if defined(ANGLE)
  const double*get_center(Triangle*t)const
  {
    if(t->is_modified())
    {
      t->xor_modified();

#if defined(L2)
       update_circumcircle(t);
#endif
#if defined(Linf)
      update_circumsquare(t);
#endif      
    }
    return t->get_center();
  }
  
#else

  const double*get_center(Triangle*t)const
  {
    if(t->is_modified())
    {
      t->xor_modified();

#if defined(L2)
       update_circumcircle(t);
#endif
       
#if defined(Linf)
      update_circumsquare(t);
#endif 
      
    }
    return t->get_center();
  }
#endif

    double incircle(Triangle*t,int p)const
    {
      
#if defined(POWER)
	double d = -power_distance(get_center(t),v[p].p) + v[p].w;
#else
	double d = -power_distance(get_center(t),v[p].p);
#endif
// 	if (fabs(d) < 1e-12)
// 	  d = 0.0;

#ifdef DEBUG
        std::cerr<<"incircle test: triangle (";
        print_triangle(std::cerr,t)<<") with vertex "<<p<<" ("<<v[p].p[0]<<" "<<v[p].p[1]<<") ==> "<<d<<"\n";
#endif

        return d;
    }
    
  double insquare(Triangle*t,int p, double angle) const
  {
      double square[4];
      double d;
            
      circumsquare(t, square, angle);
      
// 	printf("--> %lf %lf\n", t->get_angle(), _angle);	
      
      if (square[0] == INFINITY)
	  d = INFINITY;
      else
      {
	  double pp[2];
	  Vertex *pv = v+p;
	  pp[0] = pv->p[0]*cos(angle) + pv->p[1]*sin(angle);
	  pp[1] = -pv->p[0]*sin(angle) + pv->p[1]*cos(angle);

#if defined(POWER)
	  d = -(max_distance(square,pp)) + v[p].w;
#else
	  d = -(max_distance(square,pp));
#endif
      }
// #ifdef DEBUG
//         std::cerr<<printf("insquare test: triangle (";
//         print_triangle(std::cerr,t)<<") with vertex "<<p<<" ("<<p[0]<<" "<<p[1]<<") ==> "<<d<<"\n";
// //         printf("square : %lf %lf, %lf %lf\n", get_center(t)[0]-get_center(t)[2], get_center(t)[1]-get_center(t)[2], get_center(t)[0]+get_center(t)[2], get_center(t)[1]+get_center(t)[2]);
// // 	printf("dist from center %lf, radius = %lf\n", max_distance(get_center(t),v[p].p), get_center(t)[2]);
// #endif
      return d;
  } 
    
  double insquare(Triangle*t,int p)const
  {
      double d;
      if (get_center(t)[0] == INFINITY)
	  d = INFINITY;
      else
      {
#if defined(POWER)	  
	  d = -(max_distance(get_center(t),v[p].p)) + v[p].w;
#else
	  d = -(max_distance(get_center(t),v[p].p));
#endif    
      }
#ifdef DEBUG
      std::cerr<<printf("insquare test: triangle (";
      print_triangle(std::cerr,t)<<") with vertex "<<p<<" ("<<v[p].p[0]<<" "<<v[p].p[1]<<") ==> "<<d<<"\n";
//         printf("square : %lf %lf, %lf %lf\n", get_center(t)[0]-get_center(t)[2], get_center(t)[1]-get_center(t)[2], get_center(t)[0]+get_center(t)[2], get_center(t)[1]+get_center(t)[2]);
// 	printf("dist from center %lf, radius = %lf\n", max_distance(get_center(t),v[p].p), get_center(t)[2]);
#endif
      return d;
  } 

  double power_distance(const double*a,const double*b)const
  {
    return sqrsum(b[0]-a[0],b[1]-a[1])-a[2];
  }

  double power_distance(const double*a,int i)const
  {
    return power_distance(a,v[i].p);
  }

  double power_distance(Triangle*t,int j)const
  {
    return power_distance(get_center(t),j);
  }
  
  double max_distance(const double*a,const double*b)const
  {
    if (fabs(b[0]-a[0]) >= fabs(b[1]-a[1]))
      return fabs(b[0]-a[0])-a[2];
    else
      return fabs(b[1]-a[1])-a[3];
  }

  double max_distance(const double*a,int i)const
  {
    return max_distance(a,v[i].p);
  }

  double max_distance(Triangle*t,int j)const
  {
    return max_distance(get_center(t),j);
  }
  
  double volume(Triangle *t);

  int locate(double*q,Triangle*&t,int&i);

  double lsize (double *q);

  long leftoncount;
  double lefton_test(double*q,Triangle*t,int i)
  {
    ++leftoncount;
    return ccw(q,v[org(t,i)].p,v[dest(t,i)].p);
  }

  double lefton_test(int i,int j,int k)
  {
    ++leftoncount;
    return ccw(v[i].p,v[j].p,v[k].p);
  }

  void print(const char*file)const;
  std::ostream&print_vertex(std::ostream&os,int i)const;
  std::ostream&print_triangle(std::ostream&os,Triangle*t,bool shift= false)const;

  bool is_valid() const;

  bool check_nb() const;

public :

  void getVertexPosition(int index, float & x, float &y)
  {
    x = v[index].p[0];
    y = v[index].p[1];
  }

  void getTriangleIndexes(int index, int &i, int &j, int &k)
  {
    i = tri[index].get_opp(0);
    j = tri[index].get_opp(1);
    k = tri[index].get_opp(2);
  }
public :

    Triangle * getTri()
    {
        return tri;
    }
    
    int get_triangle_index(Triangle *t)
    {
      return t-tri;
    }

    Vertex * getVertices()
    {
        return v;
    }
    
    int get_vertex_index(Vertex *pv)
    {
      return pv-v;
    }

    void set_angle(double angle)
    {
        _angle = angle;
	_cosa = cos(angle);
	_sina = sin(angle);
    }
    double get_angle()
    {
        return _angle;
    }
    
    void get_vertex_coord(Vertex *pv, double *p)
    {
        p[0] = pv->p[0]*_cosa + pv->p[1]*_sina;
        p[1] = -pv->p[0]*_sina + pv->p[1]*_cosa;
    }    

private :
  Vertex *v;
  Triangle *tri;
  int nv;
  int ntri;
  
  Triangle **free_t;
  int ifree_t;
  
  Cavity *cav;
  Background *_bg;
  
  double _factor;
  double _xmin;
  double _ymin;
  double _angle;
  double _cosa;
  double _sina;
  
  // If (t,i) and (t0,i0) refer to the same edge from two neighbor triangles, each one
  // is called a directed edge and each one is called symmetric to the other.
  void get_nb(Triangle *t, int i, Triangle *&t0, int &i0) const
  {
    t0 = t->get_adjacent(i);
    i0 = t->get_opp(i);
  }

  void set_nb(Triangle*t1,int i1,Triangle*t2,int i2)
  {
    if(t1!=NULL)
    {
      t1->set_adjacent(i1,t2);
      t1->set_opp(i1,i2);
    }

    if(t2!=NULL)
    {
      t2->set_adjacent(i2,t1);
      t2->set_opp(i2,i1);
    }
  }
  
  void turn_ccw(Triangle *&t, int &piv_t, Triangle *&s, int &piv_s)
  {
	get_nb (t, Triangle::prev (piv_t), s, piv_s);
	piv_s = Triangle::prev (piv_s);
  }
  void turn_cw(Triangle *&t, int &piv_t, Triangle *&s, int &piv_s)
  {
	get_nb (t, Triangle::next (piv_t), s, piv_s);
	piv_s = Triangle::next (piv_s);
  }
        
    
public :
  int get_vertex(Triangle*t,int i)const
  {
    return t->get_vertex(i);
  }
  
  Vertex* add_vertex(double*s)
  {
  	Vertex *pv = v+nv;
  	new_vertex(s);
  	return pv;
  }
  
  Vertex* get_vertex(int i)
  {
    return v+i;
  }
  
    void set_flag(Vertex *pv, int flag)
  {
  	 pv->flag = flag;
  }

  void inc_flag(Vertex *pv)
  {
  	 pv->flag++;
  }
  
  void dec_flag(Vertex *pv)
  {
  	 pv->flag--;
  }

  int get_index(Vertex *pv)
  {
	return pv-v;
  }

  int get_flag(Vertex *pv)
  {
  	 return pv->flag;
  }

private :
  void set_vertex(Triangle*t,int i,int v)
  {
    get_vertex(v)->t = t;
    return t->set_vertex(i,v);
  }

  Triangle* get_adjacent(Triangle*t,int i)const
  {
    return t->get_adjacent(i);
  }

  void set_adjacent(Triangle*t,int i,Triangle*s)
  {
    t->set_adjacent(i,s);
  }

  int get_index(Triangle*t,int i)const
  {
    return t->get_opp(i);
  }

  int org(Triangle*t,int i) const // premier sommet de l'arete i
  {
    return t->org(i);
  }

  int dest(Triangle*t,int i) const // deuxieme sommet de l'arete i
  {
    return t->dest(i);
  }

  int v_capacity;
  int tri_capacity;
  
#if defined (POWER)
  void new_vertex(double x,double y, double w)
  {
    assert(nv<v_capacity);
    Vertex*pv= v+nv++;
    pv->p[0]= x;
    pv->p[1]= y;
    pv->flag = 0;
#if defined (L2)
    pv->w = w*w;
#else
    pv->w = w;
#endif
//     pv->w = 0.0;
  }
#else
  void new_vertex(double x,double y)
  {
    assert(nv<v_capacity);
    Vertex*pv= v+nv++;
    pv->p[0]= x;
    pv->p[1]= y;
    pv->flag = 0;
  }
#endif


  void new_vertex(double*s)
  {
    assert(nv<v_capacity);
    Vertex*pv= v+nv++;
    pv->p[0]= s[0];
    pv->p[1]= s[1];
#if defined (POWER)
    pv->w = s[2]*s[2];
//     pv->w = 0.0;
//     printf("w --> %lf\n", pv->w);
#endif
    pv->flag = 0; 
    
#if defined (ANGLE)
    double angle = _bg->get_angle(pv->p);
    pv->a = angle;
#endif

    int ix= static_cast<int> ((s[0]-llx)*wr);
    int iy= static_cast<int> ((s[1]-lly)*hr);
    Vertex**bucket= bucket_array+ix+BUCKET_SIZE*iy;
    pv->n= (*bucket);
    (*bucket)= pv;
  }
  
    Triangle* new_triangle()
    {
      Triangle *t = NULL;
//         dprintf("ifree_t = %d\n", ifree_t);
        while (ifree_t > 0 && t == NULL)
        {
            t = free_t[--ifree_t];
            if (t->get_flag() == DEFAULT_TRIANGLE_FLAG)
                t = NULL;
        }
        if (t == NULL)
        {
            assert(ntri<tri_capacity);
            t = tri+ntri++;
        }

        t->set_members(-1, -1, -1, NULL, NULL, NULL, -1, -1, -1 );
	t->set_flag(DEFAULT_TRIANGLE_FLAG);
        t->pos2 = 0;
        return t;
    }
  
    void free_triangle(Triangle *t)
    {
        dprintf("--> erasing triangle %d --> %d %d %d\n", (int)(t-tri), t->get_vertex(0), t->get_vertex(1), t->get_vertex(2));
        free_t[ifree_t++] = t;
    }

    void clear_free_triangle()
    {
        ifree_t = 0;
    }

  Triangle *last_tri;
  int last_i;

/*   stack cstack; // stack for closers */

  double llx,lly,wr,hr;
  Vertex* find(Vertex**b,int x,int y)const;

  Vertex* bucket_array[BUCKET_SIZE*BUCKET_SIZE];  
};

class Facet
{
private:
    int _v[2];
    Facet *_adj[2];
    int _opp[2];
    
    Triangle *_t[2];
    int _f[2];
    
    Facet *_c[2];
    Facet *_p;
    
    int flag;
    
public:
    Facet()
    {
      _v[0] = _v[1] = -1;
      _adj[0] = _adj[1] = NULL;
      _opp[0] = _opp[1] = -1;
      _t[0] = _t[1] = NULL;
      _f[0] = _f[1] = -1;
      _c[0] = _c[1] = NULL;
      _p = NULL;
      flag = 0;
    }
    void set_members(int v0, int v1, Facet *f0, Facet *f1, int i0, int i1)
    {
        _v[0]= v0;
        _v[1]= v1;
        _adj[0]= f0;
        _adj[1]= f1;
        _opp[0]= i0;
        _opp[1]= i1;
    }
    
    void set_nb(int i, Facet *f, int opp)
    {
      _adj[i] = f;
      _opp[i] = opp;
    }
    
    void set_int_triangle(Triangle *t, int face)
    {
      _t[0] = t;
      _f[0] = face;
    }
    void set_ext_triangle(Triangle *t, int face)
    {
      _t[1] = t;
      _f[1] = face;
    }
    void set_parent(Facet *p)
    {
      _p = p;
    }
    void set_children(Facet *c0, Facet *c1)
    {
      _c[0] = c0;
      _c[1] = c1;
    }
    void reverse()
    {
        int itmp;
        itmp = _v[1];
        _v[1]= _v[2];
        _v[2] = itmp;

        itmp = _opp[1];
        _opp[1]= _opp[2];
        _opp[2] = itmp;
	
	itmp = _f[0];
        _f[0]= _f[1];
        _f[1] = itmp;

        Facet *ftmp;
        ftmp = _adj[1];
        _adj[1]= _adj[2];
        _adj[2] = ftmp;

        ftmp = _c[1];
        _c[1]= _c[2];
        _c[2] = ftmp;

        Triangle *ttmp;
        ttmp = _t[0];
        _t[0]= _t[1];
        _t[1] = ttmp;
    }
    void get_vertices(int &v0, int &v1)
    {
      v0 = _v[0];
      v1 = _v[1];
    }
    void get_neigh(Facet *&n0, Facet *&n1, int &opp0, int &opp1)
    {
      n0 = _adj[0];
      n1 = _adj[1];
      opp0 = _opp[0];
      opp1 = _opp[1];
    }
    Facet *get_adj(int i)
    {
      return _adj[i];
    }
    int get_opp(int i)
    {
      return _opp[i];
    }
    int get_vertex(int i)
    {
      return _v[i];
    }
    Triangle *get_int_triangle()
    {
      return _t[0];
    }
    int get_int_face()
    {
      return _f[0];
    }
    Triangle *get_ext_triangle()
    {
      return _t[1];
    }
    int get_ext_face()
    {
      return _f[1];
    }
    Facet *get_parent()
    {
      return _p;
    }
    Facet *get_child(int i)
    {
      return _c[i];
    }
    void set_flag(int i)
    {
      flag = i;
    }
    int get_flag()
    {
      return flag;
    }
};

class Cavity
{
    Triangulation *tri;
    
   Triangle *tbase;
   int fbase;
   int center;   
   
    Triangle *triangle_cavity[CAVITY_CAPACITY];
    int itcr;
    int itcw;    
     
    Vertex *vertex_cavity[CAVITY_CAPACITY];
    int ivcr;
    int ivcw; 
    
    Triangle *unmodified[CAVITY_CAPACITY];
    int iunmodr;
    int iunmodw;
    
    Facet cav_f[CAVITY_CAPACITY];
    int icavfw;
    int icavfr;
    
    Facet *free_f[CAVITY_CAPACITY];
    int ifree_f;
    
//     Facet *stack[CAVITY_CAPACITY];
//     int stackr;
//     int stackw;
    
    int old_triangle_flag;
    
public:    
    Cavity(Triangulation *tr)
    {
      tri = tr;
      center = -1;
      tbase = NULL;
      fbase = -1;
      itcr = 0;
      itcw = 0;
      ivcr = 0;
      ivcw = 0;
      icavfw = 0;
      iunmodr = 0;
      iunmodw = 0;
      ifree_f = 0;    
    }
    
//     void set_base(Triangle *t, int i, int p)
//     {
//       center = p;
//       tbase = t;
//       fbase = -1;
//       t->set_flag(get_new_flag());
//       add_triangle_cavity(t);
//     }
    void set_base(Triangle *t, int i, int p)
    {
        center = p;
        tbase = t;
        fbase = i;
        t->set_flag(get_new_flag());
        add_triangle_cavity(t);
        if (i != -1)
        {
            Triangle *s = t->get_adjacent(i);
            s->set_flag(get_new_flag());
            add_triangle_cavity(s);
        }
    }
//     Triangle *get_base()
//     {
//       return tbase;
//     }
    int is_base(Triangle *t)
    {
      if (t == tbase)
	return true;
      if (fbase != -1 && t == tbase->get_adjacent(fbase))
	return true;
      return false;
    }
    
    void init();
    
    void star_facet(Facet *fti);
    void recover_parent_facet(Facet *fti);
    
    void set_old_flag(int i)
    {
      old_triangle_flag = i;
    }  
    void inc_old_flag()
    {
      old_triangle_flag++;
      dprintf("new flag --> %d\n", old_triangle_flag);
    }    
    int get_old_flag()
    {
      return old_triangle_flag;
    }
    int get_new_flag()
    {
      return old_triangle_flag+1;
    }    
    void remove_facet(Facet *fti);
    Facet *restore_parent_facet(Facet *fti);
    
    void init_read_triangle_cavity()
    {
      itcr = 0;
    }
    void add_triangle_cavity(Triangle *t)
    {
      assert(t->get_flag() == get_new_flag());
      dprintf("adding tetra %d in cavity\n", tri->get_triangle_index(t));
      triangle_cavity[itcw++] = t;
    } 
    Triangle *next_triangle_cavity()
    {
        if (itcr < itcw)
            return triangle_cavity[itcr++];
        else
            return NULL;
    }
    int triangle_cavity_size()
    {
      return itcw;
    }
    void clear_triangle_cavity()
    { 
        itcr = 0;
        itcw = 0;
    }
    
    void init_read_vertex_cavity()
    {
      ivcr = 0;
    }
    void add_vertex_cavity(int p)
    {       
        Vertex *v = tri->get_vertex(p);
        if (v->flag == 0)
        {
            v->flag = 2;
            vertex_cavity[ivcw++] = v;
        }
//         else
//             v->flag += 3;
    }
    void inc_vertex_cavity_flag(int p)
    {
      Vertex *v = tri->get_vertex(p);
      v->flag++;
//       printf("inc vertex %d (%p) --> flag = %d\n", p, v, v->flag);
    }
    void dec_vertex_cavity_flag(int p)
    {
      Vertex *v = tri->get_vertex(p);
      v->flag--;
    }
    Vertex *next_vertex_cavity()
    {
      if (ivcr < ivcw)
	return vertex_cavity[ivcr++];
	else
	return NULL;
    }
    int vertex_cavity_size()
    {
      return ivcw;
    }    
    void clear_vertex_cavity()
    { 
      for (int i = 0; i < ivcw; i++)
	vertex_cavity[i]->flag = 0;
        ivcr = 0;
        ivcw = 0;
    }
    
    void init_unmodified()
    {
        iunmodr = 0;
        iunmodw = 0;
    }
    void add_unmodified_triangle(Triangle *t)
    {
      assert(t->get_flag() == -1);
      unmodified[iunmodw++] = t;
    } 
    Triangle *next_unmodified()
    {
      if (iunmodr < iunmodw)
	return unmodified[iunmodr++];
            else
	return NULL;
    }
    int unmodified_size()
    {
      return iunmodw;
    }
    Facet *get_facet(int i)
    {
	return &cav_f[i];
    }
    int get_nb_facet()
    {
        return icavfw;
    }    
    void init_read_facet()
    {
      icavfr = 0;
    }
    Facet *get_next_facet()
    {
      if (icavfr < icavfw)
	return &cav_f[icavfr++];
      else
	return NULL;
    }
    int get_index(Facet *ft)
    {
      return (int)(ft-cav_f);
    }
    void free_facet(Facet *ft)
    {
//       free_f[ifree_f++] = ft;
	assert(ft->get_flag() == 1);
    }
    Facet *new_facet()
    {
//       if (ifree_f > 0)
// 	return free_f[--ifree_f];
//       else
	return &cav_f[icavfw++];
    }
    void clear_cavity_facet()
    {
        icavfw = 0;
	icavfr = 0;
        ifree_f = 0;
    }    
//     void init_stack()
//     {
//         stackr = 0;
//         stackw = 0;
//     }
//     Facet *stack_next()
//     {
//       printf("read stack --> %d\n", stackr);
//       if (stackr < stackw)
// 	return stack[stackr++];
//       else
// 	return NULL;
//     }    
//     void add_to_stack(Facet *ft)
//     {
//       stack[stackw++] = ft;
//     }
//     int stack_size()
//     {
//       return stackw;
//     }
    
    void dump();
    void check();
    void export_vtk();
    
};

class progress_bar
{

protected:
    int _max;
    int _inc;
    int _last;

public:

    progress_bar()
    {
        _max = 100;
        _inc = 1;
        _last = 0;
    }
    progress_bar(int max)
    {
        _max = max;
        _inc = 1;
        _last = 0;
        if (_max > 100)
            _inc = floor((double)max/100.0);
    }

    inline void progress(int i)
    {
        if (i >=_last)
        {
            _last += _inc;
            printf("\033[M\033[A\033[M");
            int percent = (int)((double)i/(double)_max*100.0);
            printf("%d %% : ", percent);
            for (int j = 0; j < percent; j++)
                printf("#");
            printf("\n");
        }
    }
    inline void start()
    {
        printf("\n");
    }
    inline void start(int max)
    {
        _max = max;
        _inc = 1;
        _last = 0;
        if (_max > 100)
            _inc = floor((double)max/100.0);
        printf("\n");
    }
    inline double end()
    {
        printf("\033[M\033[A\033[M");
    }
};

#endif