#include "dcel.h"
#include "nutil.h"
#include <map>
#include <set>

#if __cplusplus >= 201103L
#define USEHASH
#include <unordered_map>
#endif

int dcelcallbacks::key_callback(const char *txt)
{
  if (txt[0]=='r')
  {
//  loop_subdivision(graph,graph2);
//  graph2.display(data,true,true,true);
    if (select>=0)
    {
      std::cout << select << std::endl;
      dcel* gg2= new dcel;
      dcel* gg=objects[select];
      doo_sabin_subdivision(*gg,*gg2); 
      objects[select]=gg2;
    }
    draw();
    return 1;
  }
  if (txt[0]=='c')
  {
    if (select>=0)
    {
      std::cout << select << std::endl;
      dcel* gg2= new dcel(*objects[select]);
      add_entity(gg2);
      draw();
    }
    return 1;
  }

  if (txt[0]=='p')
  {
    int sz=objects.size();
    for (int i=0;i<sz;++i)
    {
      objects[i]->print(std::cout);
    }
    return 1;
  }

  if (txt[0]=='d')
  {
    if (select>=0)
    {
      objects.erase(objects.begin()+select);
      draw();
    }
    return 1;
  }
  return 0;
}

int dcelcallbacks::pick_callback(int nb, pickinfo picks[])
{
  std::cout << "num=" << nb << std::endl;
  drag=false;
  select=-1;
  if (nb)
  {
    if (picks[0].nb_uid)
    {
      int objid=picks[0].uids[0];
      std::cout << "pick objid=" << objid << std::endl;
      picked=objects[objid];
      listp.clear();
      listp.reserve(picked->nb_vertices());
      for (vertex_handle it=picked->begin_vertices();it!=picked->end_vertices();++it)
        listp.push_back(it->point());
      return 0; // first element in the stack (closest to the eye) is chosen (-1 if no pick want
    }
  }
  return -1;
}


int dcelcallbacks::drag_callback(npoint3 orig,npoint3 p, pickinfo pick)
{
  int ret=0;
  select=-1;
  drag=true;
  npoint delta=npoint(p)-npoint(orig);
  std::cout << picked <<  "drag"  << pick.nb_uid <<  std::endl;
  if (picked)
  {
    if (pick.nb_uid>1) // second uid = point number
    {
      // move point
      int objid=pick.uids[0];
      int ptid=pick.uids[1];
      std::cout << "drag objid=" << objid << std::endl;
      vertex_handle it=picked->begin_vertices();
      for (int i=0;i<listp.size();++i,++it)
      {
        if (i==ptid)
        {
          npoint3 P=listp[i];
          it->point()=P+delta;
        }
//        objects[objid]->set_CP(i,picked->CP(i)+delta*P[3]);
      }
      draw(objid);
      ret=1;
    }
    else if (pick.nb_uid<=1)
    {
//      std::cout << "move part" << std::endl;
      // move whole part
      int objid=pick.uids[0];
      std::cout << "drag objid=" << objid << std::endl;
      vertex_handle it=picked->begin_vertices();
      for (int i=0;i<listp.size();++i,++it)
      {
        npoint3 P=listp[i];
        it->point()=P+delta;
//        objects[objid]->set_CP(i,picked->CP(i)+delta*P[3]);
      }
      draw(objid);
      ret=1;
    }
  }
  return ret;
}



int dcelcallbacks::release_callback(npoint3 orig,npoint3 p, pickinfo pick)
{
  int ret=0;
  if (picked)
  {
    if (drag==false)
    {
      if (pick.nb_uid>=1)
        select=pick.uids[0];
      else
        select=-1;
    }
    else
    {
      ret=drag_callback(orig,p,pick);
      picked=NULL;
      select=-1;
    }
    std::cout << "release : pick=" << select <<  std::endl;
  }
  return ret;
}

void dcelcallbacks::draw(void)
{
  for (int i=0;i<objects.size();++i) draw(i);
}

void dcelcallbacks::draw(int i)
{
  datas[i].clear();
  properties prop=prop_objects[i];
  prop.nb_uid=1;
  prop.uids[0]=i;
  datas[i].setpropall(prop);
  objects[i]->Display(datas[i]);
}




struct intpair
{
  intpair(int aa,int bb) {if (aa<bb) {a=aa;b=bb;} else {a=bb;b=aa;} }
  bool operator < (const intpair& other) const
  {
    if (a<other.a) return true;
    if (a>other.a) return false;
    if (b<other.b) return true;
    return false;
  }
  bool operator == (const intpair& other) const
  {
    if (a!=other.a) return false;
    if (b!=other.b) return false;
    return true;
  }
  int a,b;
};


#ifdef USEHASH
struct intpairhash
{
  size_t operator()(const intpair &p) const
  {
//    return (113 + std::hash<int>()(p.a)) * 113 + std::hash<int>()(p.b);
      return (p.a);
  }
};

  typedef std::unordered_map<intpair,std::vector<halfedge_handle>,intpairhash > edgevectype;
  typedef std::unordered_map<intpair,std::vector<halfedge_handle>,intpairhash >::iterator edgevectype_iterator;
  typedef std::unordered_map<intpair,int,intpairhash> edgescaltype;
  typedef std::unordered_map<intpair,int,intpairhash>::iterator edgescaltype_iterator;
#else  
  typedef std::map<intpair,std::vector<halfedge_handle> > edgevectype;
  typedef std::map<intpair,std::vector<halfedge_handle> >::iterator edgevectype_iterator;
  typedef std::map<intpair,int> edgescaltype;
  typedef std::map<intpair,int>::iterator edgescaltype_iterator;
#endif


dcel& dcel::operator=( const dcel& other) // deep copy
{
  clear();
  std::map<dcel_vertex*,vertex_handle> p2p;
  std::map<dcel_halfedge*,halfedge_handle> h2h;
  std::map<dcel_face*,face_handle> f2f;
  //create entities and store correspondence tables
  for (vertex_handle it=other.begin_vertices();it!=other.end_vertices();++it)
  {
    p2p[&*it]=add_vertex(*it);
  }

  for (halfedge_handle it=other.begin_halfedges();it!=other.end_halfedges();++it)
  {
    h2h[&*it]=add_halfedge(*it);
  }

  for (face_handle it=other.begin_faces();it!=other.end_faces();++it)
  {
    f2f[&*it]=add_face(*it);
  }
  std::cout << other.nv << " " << other.nh << " "<<  other.nf << std::endl;
  std::cout <<p2p.size() << " " << h2h.size() << " "<<  f2f.size() << std::endl;
  std::cout << nv << " " << nh << " "<<  nf << std::endl;
  
  // populate links

//   void set_incident_halfedge(halfedge_handle inc_)

  for (vertex_handle it=begin_vertices();it!=end_vertices();++it)
  {
    halfedge_handle incident_=it->incident_halfedge();
    if (incident_!=halfedge_handle()) it->set_incident_halfedge(h2h[&*incident_]); else it->set_incident_halfedge(halfedge_handle());
  }


//   void set_origin(vertex_handle ori_)
//   void set_twin(halfedge_handle tw_) 
//   void set_incident_face(face_handle inc_)
//   void set_next_halfedge(halfedge_handle next_)
//   void set_prev_halfedge(halfedge_handle prev_)

  for (halfedge_handle it=begin_halfedges();it!=end_halfedges();++it)
  {
    vertex_handle origin_=it->origin();
    halfedge_handle twin_=it->twin();
    face_handle face_=it->incident_face();
    halfedge_handle prev_=it->prev_halfedge();
    halfedge_handle next_=it->next_halfedge();
    if (origin_!=vertex_handle()) it->set_origin(p2p[&*origin_]); else it->set_origin(vertex_handle());
    if (twin_!=halfedge_handle()) it->set_twin(h2h[&*twin_]); else it->set_twin(halfedge_handle());
    if (face_!=face_handle()) it->set_incident_face(f2f[&*face_]); else  it->set_incident_face(face_handle());
    if (next_!=halfedge_handle()) it->set_next_halfedge(h2h[&*next_]); else it->set_next_halfedge(halfedge_handle()); 
    if (prev_!=halfedge_handle()) it->set_prev_halfedge(h2h[&*prev_]); else it->set_prev_halfedge(halfedge_handle());
  }


//   halfedge_handle outer_loop(void) const { return out ;}
//   halfedge_handle inner_loop(int i) const { return in[i] ;}
//   int nb_inner_loop(void) const { return in.size() ;}
//   void set_outer_loop(halfedge_handle out_) { out=out_ ;}
//   void add_inner_loop(halfedge_handle in_) { in.push_back(in_);}
//   void clear_inner_loops(void) { in.clear();}

  for (face_handle it=begin_faces();it!=end_faces();++it)
  {
    halfedge_handle outer_=it->outer_loop();
    if (outer_!=halfedge_handle()) it->set_outer_loop(h2h[&*outer_]); else it->set_outer_loop(halfedge_handle());
    for (int i=0;i<it->nb_inner_loop();++i)
    {
      halfedge_handle inner_=it->inner_loop(i);
      if (inner_!=halfedge_handle()) it->set_inner_loop(i,h2h[&*inner_]); else it->set_inner_loop(i,halfedge_handle());
    }
  }
  return *this;
}

dcel* dcel::clone(void) const
{
  return new dcel(*this);
}

dcel::~dcel(void)
{ 
  clear(); 
}

void dcel::clear(void)
{ 
  for (std::list<dcel_halfedge *>::iterator it=halfedges.begin();it!=halfedges.end();++it) delete *it;
  for (std::list<dcel_face *>::iterator it=faces.begin();it!=faces.end();++it) delete *it;
  for (std::list<dcel_vertex *>::iterator it=vertices.begin();it!=vertices.end();++it) delete *it;
  halfedges.clear();
  faces.clear();
  vertices.clear();
  nv=0;nf=0;nh=0;
}


void dcel::_display(data_container &data, bool vertices,bool lines, bool faces) const // for display purposes
{
  if (lines)
  for (halfedge_handle heh=begin_halfedges();heh!=end_halfedges();++heh)
  {
    vertex_handle v1h=heh->origin(),v2h;
    if (heh->twin()!=halfedge_handle()) 
      v2h=heh->twin()->origin();
    else
      if (heh->next_halfedge()!=halfedge_handle())
        v2h=heh->next_halfedge()->origin();
      else
        break;
    line l;
    l.pts[0]=v1h->point();
    l.pts[1]=v2h->point();
    data.add_line(l);
  }
  if (vertices)
  {
    properties p,p1;
    p=data.getproppoints();
    p.pointsize=5;
    p1=p;
    if (p.nb_uid<2) p.nb_uid=2;
    int i=0;
    for (vertex_handle vh=begin_vertices();vh!=end_vertices();++vh,++i)
    {
      p.uids[1]=i;
      data.setproppoints(p);
      data.add_point(vh->point());
    }
    data.setproppoints(p1);
  }
  if (faces)
  for (face_handle fh=begin_faces();fh!=end_faces();++fh)
  {
    npoint3 G(0.,0.,0.);
    int nb=0;
    halfedge_handle heh=fh->outer_loop();
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        npoint3 pt=heh->origin()->point();
        G+=pt;
        nb++;
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
    }
    G/=(double)nb;
    heh=fh->outer_loop();
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        npoint3 pt1=heh->origin()->point();
        heh=heh->next_halfedge();
        npoint3 pt2=heh->origin()->point();
        triangle t;
        t.pts[0]=G;
        t.pts[1]=pt1;
        t.pts[2]=pt2;
        data.add_triangle(t);
      }
      while(heh!=hbegin);
    }
  }
}

void dcel::print(std::ostream &os,int prec) const 
{
  os << "nv=" << nv << " nh=" << nh << " nf=" << nf << std::endl;
  os << "all vertices : " << std::endl;
  for (vertex_handle vh=begin_vertices();vh!=end_vertices();++vh)
  {
    os << "vertex " << vh->num << " " ;
    vh->point().print(os);
  }
  os << "all faces : " << std::endl;
  for (face_handle fh=begin_faces();fh!=end_faces();++fh)
  {
    os << "face " << fh->num << std::endl;
    halfedge_handle heh=fh->outer_loop();
    if (heh==halfedge_handle())
    {
      os << " no exterior loop : exterior face " << std::endl;
    }
    else
    {
      halfedge_handle hbegin=heh;
      os << " exterior loop : ";
      do
      {
        os << heh->origin()->num << " ";
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
      os << std::endl;
    }
    for (int i=0;i<fh->nb_inner_loop();++i)
    {
      heh=fh->inner_loop(i);
      halfedge_handle hbegin=heh;
      os << " interior loop : ";
      do
      {
        os << heh->origin()->num << " ";
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
      os << std::endl;
    }
  }
  os << std::endl;
}






// build dcel from list of points and faces (no internal loops). faces described by ordered point id as list (nbpt / pt1... ptX / etc...)
// point list begins with id 0 and is dense
void build_dcel_from_mesh(std::vector<npoint3> vp,std::vector<int> vf,dcel& graph)
{
  std::vector<vertex_handle> vhp;
  vhp.reserve(vp.size());
  int cnt=0;
  for (std::vector<npoint3>::iterator ip=vp.begin();ip!=vp.end();++ip)
  {
    vertex_handle v=graph.add_vertex(dcel_vertex());
    v->point()=*ip;
    v->num=cnt++;
    vhp.push_back(v);
  }
  cnt=0;
  int cnt2=0;

  edgevectype em;

  for (std::vector<int>::iterator iff=vf.begin();iff!=vf.end();)
  {
    int nb=*iff++;
    halfedge_handle hh2;
    halfedge_handle init;
    face_handle fh=graph.add_face(dcel_face());
    fh->num=cnt++;
    if (nb)
    {
      int orig=*iff;
      for (int p=0;p<nb;++p)
      {
        halfedge_handle hh=graph.add_halfedge(dcel_halfedge());
        hh->num=cnt2++;
        int pt=*iff++;
        int nextpt=(p+1==nb?orig:*iff);
        hh->set_origin(vhp[pt]);
        vhp[pt]->set_incident_halfedge(hh);
        if (p)
        {
          hh->set_prev_halfedge(hh2);
          hh2->set_next_halfedge(hh);
        }
        else 
        { 
          init=hh;
          fh->set_outer_loop(init);
        }
          
        hh->set_incident_face(fh);
        std::pair<edgevectype_iterator,bool> item=em.insert(make_pair(intpair(pt,nextpt),std::vector<halfedge_handle>())); 
        item.first->second.push_back(hh);
        if (item.first->second.size()==2)
        {
          item.first->second[0]->set_twin(item.first->second[1]);
          item.first->second[1]->set_twin(item.first->second[0]);
          em.erase(item.first);
        }
        hh2=hh; 
      }
      init->set_prev_halfedge(hh2);
      hh2->set_next_halfedge(init);
    }
  }
}




void doo_sabin_subdivision(dcel& graph, dcel& subd_graph)
{
  int nf=graph.nb_faces();
  int nv=graph.nb_vertices();
  int nh=graph.nb_halfedges();
  
  std::vector<int> hedgv(nh);
  std::vector<npoint3> vec;
  vec.reserve(nh);
  std::vector<int> ele;
  ele.reserve(nh+nv+nf+nh+nh/2);
  int nm=0;
  for (face_handle fh=graph.begin_faces();fh!=graph.end_faces();++fh)
  {
    halfedge_handle heh=fh->outer_loop();
    npoint3 G(0.,0.,0.);
    std::vector<npoint3> Ge;
    int nb=0;
  
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        G+=heh->origin()->point();
        npoint3 gh(0.,0.,0);
        gh+=heh->origin()->point();
        nb++;
        heh=heh->next_halfedge();
        gh+=heh->origin()->point();
        gh/=2.0;
        Ge.push_back(gh);
      }
      while(heh!=hbegin);
      G/=nb;
      
      heh=fh->outer_loop();
      double cnt=0;
      do
      {
        npoint3 gg(0.,0.,0);
        int nhv=heh->num;
        gg+=heh->origin()->point();
        gg+=G;
        gg+=Ge[cnt];
        gg+=Ge[(cnt==0?nb-1:cnt-1)];
        heh=heh->next_halfedge();
        gg/=4.0;
        hedgv[nhv]=nm++;
        vec.push_back(gg);
        cnt++;
      }
      while(heh!=hbegin);
    }
  }

  for (face_handle fh=graph.begin_faces();fh!=graph.end_faces();++fh)
  {
    halfedge_handle heh=fh->outer_loop();
    
    if (heh!=halfedge_handle())
    {
      int id=ele.size();
      ele.push_back(0);
      halfedge_handle hbegin=heh;
      do
      {
        ele[id]++;
        ele.push_back(hedgv[heh->num]);
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
    }
  }
  
  for (vertex_handle vh=graph.begin_vertices();vh!=graph.end_vertices();++vh)
  {
    halfedge_handle heh=vh->incident_halfedge();
    bool success=true;
    int id=ele.size();
    if (heh!=halfedge_handle())
    {

      ele.push_back(0);
      halfedge_handle hbegin=heh;
      do
      {
        ele[id]++;
        ele.push_back(hedgv[heh->num]);
        heh=heh->prev_halfedge();
        heh=heh->twin();
        if (heh==halfedge_handle()) {success=false;break;}
      }
      while(heh!=hbegin);
    }
    else
      success=false;
    if (!success) ele.resize(id);
  }
  std::set<intpair> ss; 
  for (halfedge_handle hh=graph.begin_halfedges();hh!=graph.end_halfedges();++hh)
  {
    if(hh->twin() !=halfedge_handle())
    {
      int num1=hh->num;
      int num2=hh->twin()->num;
      std::pair<std::set<intpair>::iterator,bool> p=ss.insert(intpair(num1,num2));
      if (!p.second)
      {
        ele.push_back(4);
        ele.push_back(hedgv[hh->num]);
        ele.push_back(hedgv[hh->twin()->next_halfedge()->num]);
        ele.push_back(hedgv[hh->twin()->num]);
        ele.push_back(hedgv[hh->next_halfedge()->num]);
      }
    }
  }
  build_dcel_from_mesh(vec,ele,subd_graph);
}


void catmull_clark_subdivision(dcel& graph, dcel& subd_graph)
{
  int nf=graph.nb_faces();
  int nv=graph.nb_vertices();
  int nh=graph.nb_halfedges();

  edgescaltype edgv;
  
  std::vector<int> facv(nf);
  std::vector<int> vertv(nv);
  std::vector<npoint3> vec;
  vec.reserve(nf+nh/2+nv);
  std::vector<int> ele;
  ele.reserve(nh+nv+nf+nh+nh/2);
  int nm=0;

  for (face_handle fh=graph.begin_faces();fh!=graph.end_faces();++fh)
  {
    halfedge_handle heh=fh->outer_loop();
    npoint3 G(0.,0.,0.);
    int nb=0;
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        G+=heh->origin()->point();
        nb++;
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
      G/=nb;
      vec.push_back(G);
      facv[fh->num]=nm++;
    }
  }
  for (halfedge_handle hh=graph.begin_halfedges();hh!=graph.end_halfedges();++hh)
  {
    npoint3 gh(0.,0.,0.);
    halfedge_handle heh=hh;
    int orig=heh->origin()->num;
    gh+=heh->origin()->point();
    gh+=vec[facv[heh->incident_face()->num]];
    heh=heh->twin();
    if (heh!=halfedge_handle())
    {
      gh+=heh->origin()->point();
      gh+=vec[facv[heh->incident_face()->num]];
      int dest=heh->origin()->num;
      gh/=4.0;
      std::pair<edgescaltype_iterator,bool> p= edgv.insert(std::make_pair(intpair(orig,dest),0));
      if (p.second)
      {
        vec.push_back(gh);
        p.first->second=nm++;
      }
    }
  }
  for (vertex_handle vh=graph.begin_vertices();vh!=graph.end_vertices();++vh)
  {
    npoint3 F(0.,0.,0.);
    npoint3 R(0.,0.,0.);
    halfedge_handle heh=vh->incident_halfedge();
    bool success=true;
    int n=0;
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        npoint3 e(0.,0.,0.);
        int orig=heh->origin()->num;
        e+=heh->origin()->point();
        F+=vec[facv[heh->incident_face()->num]];
        n++;
        heh=heh->twin();
        if (heh!=halfedge_handle())
        {
          int dest=heh->origin()->num;
          e+=heh->origin()->point();
          e/=2;
          R+=e;
          heh=heh->next_halfedge();
        }
        else
        {
          success=false;
          break;
        }
      }
      while(heh!=hbegin);
    }
    else
      success=false;
    if (success)
    {
      F/=n;
      R/=n;
      npoint3 K=(F+2*R+(n-3)*vh->point())/n;
      vec.push_back(K);
      vertv[vh->num]=nm++;
    }
    else
      vertv[vh->num]=-1;
  }
  
  for (face_handle fh=graph.begin_faces();fh!=graph.end_faces();++fh)
  {
    halfedge_handle heh=fh->outer_loop();
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        bool success=true;
        int ori=ele.size();
        ele.push_back(4);
        ele.push_back(facv[fh->num]);
        int origin=heh->origin()->num;
        if (heh->twin()!=halfedge_handle())
        {
          int dest=heh->twin()->origin()->num;
          edgescaltype_iterator it=edgv.find(intpair(origin,dest));
          if (it!=edgv.end())
          {
            ele.push_back(it->second);
            halfedge_handle next=heh->next_halfedge();
            int val=vertv[next->origin()->num];
            if (val>=0)
            {
              ele.push_back(vertv[next->origin()->num]);
              origin=next->origin()->num;
              if (next->twin()!=halfedge_handle())
              {
                dest=next->twin()->origin()->num;
                it=edgv.find(intpair(origin,dest));
                if (it!=edgv.end())
                {
                  ele.push_back(it->second);
                }
                else success=false;
              }
              else success=false;
            }
            else success=false;
          }
          else success=false;
        } 
        else success=false;
        if (!success) 
          ele.resize(ori);
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
    }
  }
  build_dcel_from_mesh(vec,ele,subd_graph);
}



void loop_subdivision(dcel& graph, dcel& subd_graph)
{
  int nf=graph.nb_faces();
  int nv=graph.nb_vertices();
  int nh=graph.nb_halfedges();
  edgescaltype edgv;
  std::vector<int> facv(nf);
  std::vector<int> vertv(nv);
  std::vector<npoint3> vec;
  std::vector<npoint3> vef;
  vec.reserve(nh/2+nv);
  std::vector<int> ele;
  ele.reserve(4*nf);
  int nm=0;
  int nmf=0;
  for (face_handle fh=graph.begin_faces();fh!=graph.end_faces();++fh)
  {
    halfedge_handle heh=fh->outer_loop();
    npoint3 G(0.,0.,0.);
    int nb=0;
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        G+=heh->origin()->point();
        nb++;
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
      G/=nb;
      vef.push_back(G);
//      vec.push_back(G);
      facv[fh->num]=nmf++;
    }
  }
  for (halfedge_handle hh=graph.begin_halfedges();hh!=graph.end_halfedges();++hh)
  {
    npoint3 gh(0.,0.,0.);
    halfedge_handle heh=hh;
    int orig=heh->origin()->num;
    gh+=heh->origin()->point()*1./8.;
    gh+=vef[facv[heh->incident_face()->num]]*3./8.;
    heh=heh->twin();
    if (heh!=halfedge_handle())
    {
      gh+=heh->origin()->point()*1./8.;
      gh+=vef[facv[heh->incident_face()->num]]*3./8.;
      int dest=heh->origin()->num;
      std::pair<edgescaltype_iterator,bool> p= edgv.insert(std::make_pair(intpair(orig,dest),0));
      if (p.second)
      {
        vec.push_back(gh);
        p.first->second=nm++;
      }
    }
  }
  for (vertex_handle vh=graph.begin_vertices();vh!=graph.end_vertices();++vh)
  {
    npoint3 M(0.,0.,0.);
    halfedge_handle heh=vh->incident_halfedge();
    bool success=true;
    int n=0;
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      do
      {
        heh=heh->twin();
        if (heh!=halfedge_handle())
        {
          n++;
          int dest=heh->origin()->num;
          M+=heh->origin()->point();
          heh=heh->next_halfedge();
        }
        else
        {
          success=false;
          break;
        }
      }
      while(heh!=hbegin);
    }
    else
      success=false;
    if (success)
    {
      double v=3.+2.*cos(2*3.1415926536/n);
      double beta=beta=5./8. - v*v/64.;
//      std::cout << beta << " " << n << " " ;
      npoint3 K=vh->point()*(1.-beta)+M*beta/n;
      vec.push_back(K);
      vertv[vh->num]=nm++;
//      std::cout << "yé" << std::endl;
//      K.print(std::cout);
    }
    else
    {
      vertv[vh->num]=-1;
//      std::cout << "NON" << std::endl;
    }
  }
  

  for (face_handle fh=graph.begin_faces();fh!=graph.end_faces();++fh)
  {
    halfedge_handle heh=fh->outer_loop();
    if (heh!=halfedge_handle())
    {
      halfedge_handle hbegin=heh;
      bool success=true;
      int ori=ele.size();
      ele.push_back(0);
      int nel=0;
      do
      {
        int origin=heh->origin()->num;
        if (heh->twin()!=halfedge_handle())
        {
          int dest=heh->twin()->origin()->num;
          edgescaltype_iterator it=edgv.find(intpair(origin,dest));
          if (it!=edgv.end())
          {
            ele.push_back(it->second);
            nel++;
            halfedge_handle next=heh->next_halfedge();
          }
          else success=false;
        } 
        else success=false;
        if (!success)
        {
          ele.resize(ori);
          break;
        }
        heh=heh->next_halfedge();
      }
      while(heh!=hbegin);
      if (success)
        ele[ori]=nel;
    }
  }
  for (vertex_handle vh=graph.begin_vertices();vh!=graph.end_vertices();++vh)
  {
    int vor=vertv[vh->num];
//    std::cout << vh->num << "v" << vor << std::endl;
    int n=0;
    if (vor>=0)
    {
      halfedge_handle hbegin=vh->incident_halfedge();
      if (hbegin!=halfedge_handle())
      {
        halfedge_handle hh=hbegin;

        int vor1=-1;
        int vor2=-1;
        int vord=-1;
        do
        {
          int origin=hh->origin()->num;
//          std::cout << "o" << hh->origin()->num  << " " ;
          hh=hh->twin();
          if (hh!=halfedge_handle())
          {
//            std::cout << "d" << hh->origin()->num  << " ";
            int dest=hh->origin()->num;
            #ifdef USEHASH
            std::unordered_map<intpair,int,intpairhash>::iterator it=edgv.find(intpair(origin,dest));
            #else  
            std::map<intpair,int>::iterator it=edgv.find(intpair(origin,dest));
            #endif            
            if (it!= edgv.end())
            {
              vor1=vor2;
              vor2=it->second;
//              std::cout << "e" << vor2 << " ";
              if (vor1!=-1)
              {
                ele.push_back(3);
                ele.push_back(vor);
                ele.push_back(vor2);
                ele.push_back(vor1);
//                std::cout << " han " << vor << " " << vor1 << " " << vor2 << std::endl;
              }
              else
              {
                vord=vor2;
              }
            }
            else
              break;
            hh=hh->next_halfedge();
            if (hh!=halfedge_handle())
            {
//              std::cout << hh->origin()->num  << " ";
              n++;
            }
            else break;
          }
          else break;
        }
        while (hh!=hbegin);
        ele.push_back(3);
        ele.push_back(vor);
        ele.push_back(vord);
        ele.push_back(vor2);
      }
    }
//    std::cout << n << std::endl;
//    break;
  }
  build_dcel_from_mesh(vec,ele,subd_graph);
}