Mesh_object.cpp

Go to the documentation of this file.

#include <Mesh_object.h>


Mesh_object::Mesh_object()
{name="NO_NAME";}

Mesh_object::~Mesh_object()
{destroy();}

Mesh_object::Mesh_object(const Mesh_object& m)
{
  connectivity = m.connectivity;
  point_set = m.point_set;
  parameters = m.parameters;
  name = m.name;
}

int Mesh_object::destroy()
{
  connectivity.destroy();
  point_set.destroy();
  parameters.resize(0);
  name="NO_NAME";
  return 0;
}


int Mesh_object::get_vertex_number() const 
{return point_set.size();}
double& Mesh_object::get_vertex(int k_vertex,int k_dim) 
{return point_set(k_vertex,k_dim);}
double& Mesh_object::get_vertex(int k_polygon,int k_vertex,int k_dim) 
{return point_set(connectivity(k_polygon,k_vertex),k_dim);}
double  Mesh_object::get_vertex(int k_vertex,int k_dim) const 
{return point_set(k_vertex,k_dim);}
double  Mesh_object::get_vertex(int k_polygon,int k_vertex,int k_dim) const 
{return point_set(connectivity(k_polygon,k_vertex),k_dim);}
const double* Mesh_object::get_vertex_p(int k_vertex)
{return point_set.get_vertex_p(k_vertex);}
const double* Mesh_object::get_vertex_p(int k_polygon,int k_vertex)
{return point_set.get_vertex_p(connectivity(k_polygon,k_vertex));}
V_3D Mesh_object::get_vertex_v3d(int k_vertex) const 
{return point_set(k_vertex);}
V_3D Mesh_object::get_vertex_v3d(int k_polygon,int k_vertex) const 
{return point_set(connectivity(k_polygon,k_vertex));}


int Mesh_object::set_vertex(int k_vertex,double x,double y,double z)
{return point_set.set_vertex(k_vertex,x,y,z);}
int Mesh_object::set_vertex(int k_vertex,int k_dim,double value)
{return point_set.set_vertex(k_vertex,k_dim,value);}
int Mesh_object::set_vertex(int k_vertex,V_3D v)
{return point_set.set_vertex(k_vertex,v);}
int Mesh_object::set_vertex(int k_vertex,double* v)
{return point_set.set_vertex(k_vertex,v);}
int Mesh_object::set_vertex(int k_vertex,int dim,double *v)
{return point_set.set_vertex(k_vertex,dim,v);}
int Mesh_object::set_vertex(const Point_set& p)
{return point_set.set_vertex(p);}



int Mesh_object::add_vertex(double x,double y,double z)
{return point_set.add_vertex(x,y,z);}
int Mesh_object::add_vertex(const V_3D& v)
{return point_set.add_vertex(v);}
int Mesh_object::add_vertex(const std::vector <double>& to_add)
{return point_set.add_vertex(to_add);}
int Mesh_object::add_vertex(const std::vector <V_3D>& to_add)
{return point_set.add_vertex(to_add);}
int Mesh_object::add_vertex(const double* v)
{return point_set.add_vertex(v);}
int Mesh_object::add_vertex(int dim,double *v)
{return point_set.add_vertex(dim,v);}
int Mesh_object::add_vertex(const Point_set& p)
{return point_set.add_vertex(p);}




int Mesh_object::add_unique(double x,double y,double z,double epsilon)
{return point_set.add_unique(x,y,z,epsilon);}
int Mesh_object::add_unique(const V_3D& v,double epsilon)
{return point_set.add_unique(v,epsilon);}
int Mesh_object::add_unique(const double* v,double epsilon)
{return point_set.add_unique(v,epsilon);}  

int Mesh_object::exist(const V_3D& v) const
{return point_set.exist(v);}
int Mesh_object::exist(const V_3D& v,double epsilon) const
{return point_set.exist(v,epsilon);}
int Mesh_object::exist(double x,double y,double z) const
{return point_set.exist(x,y,z);}
int Mesh_object::exist(double x,double y,double z,double epsilon) const
{return point_set.exist(x,y,z,epsilon);}


int Mesh_object::eliminate(int k_index)
{return point_set.eliminate(k_index);}
int Mesh_object::eliminate(const std::vector <int>& k_index)
{return point_set.eliminate(k_index);}
int Mesh_object::eliminate(const std::vector <double>& k_index)
{return point_set.eliminate(k_index);}
int Mesh_object::eliminate(double x,double y,double z)
{return point_set.eliminate(x,y,z);}
int Mesh_object::eliminate(double x,double y,double z,double epsilon)
{return point_set.eliminate(x,y,z,epsilon);}
int Mesh_object::eliminate(const V_3D& x,double epsilon)
{return point_set.eliminate(x,epsilon);}
int Mesh_object::eliminate(const V_3D& x)
{return point_set.eliminate(x);}
int Mesh_object::eliminate(const double* x,double epsilon)
{return point_set.eliminate(x,epsilon);}
int Mesh_object::eliminate(const double* x)
{return point_set.eliminate(x);}





int Mesh_object::get_index(int k_index) const
{return connectivity.get_index(k_index);}
int& Mesh_object::get_index(int k_index)
{return connectivity.get_index(k_index);}
int Mesh_object::get_index(int k_polygon,int k_index) const
{return connectivity.get_index(k_polygon,k_index);}
int& Mesh_object::get_index(int k_polygon,int k_index)
{return connectivity.get_index(k_polygon,k_index);}
std::vector <int> Mesh_object::get_index_of_polygon(int k_polygon) const
{return connectivity.get_index_of_polygon(k_polygon);}
const int* Mesh_object::get_index() const
{return connectivity.get_index();}


int Mesh_object::size_index() const
{return connectivity.size_index();}
int Mesh_object::polygon_number() const
{return connectivity.polygon_number();}
int Mesh_object::polygon_size(int k_polygon) const
{return connectivity.polygon_size(k_polygon);}
int Mesh_object::vertex_number() const
{return get_vertex_number();}


int Mesh_object::get_index_triangulated(int k_index) const
{return connectivity.get_index_triangulated(k_index);}
int& Mesh_object::get_index_triangulated(int k_index)
{return connectivity.get_index_triangulated(k_index);}
int Mesh_object::get_index_triangulated(int k_triangle,int k_index) const
{return connectivity.get_index_triangulated(k_triangle,k_index);}
int& Mesh_object::get_index_triangulated(int k_triangle,int k_index)
{return connectivity.get_index_triangulated(k_triangle,k_index);}
const int* Mesh_object::get_index_triangulated() const
{return connectivity.get_index_triangulated();}
int Mesh_object::get_triangle_number() const
{return connectivity.get_triangle_number();}
int Mesh_object::find_polygon(const std::vector <int>& polygon) const
{return connectivity.find_polygon(polygon);}



int Mesh_object::add_polygon(const std::vector <int>& polygon)
{return connectivity.add_polygon(polygon);}
int Mesh_object::add_polygon(int index_0,int index_1,int index_2,int index_3)
{return connectivity.add_polygon(index_0,index_1,index_2,index_3);}
int Mesh_object::add_polygon(int index_0,int index_1,int index_2)
{return connectivity.add_polygon(index_0,index_1,index_2);}
int Mesh_object::add_index(int index)
{return connectivity.add_index(index);}
int Mesh_object::set_polygon(int k_polygon,const std::vector <int>& polygon)
{return connectivity.set_polygon(k_polygon,polygon);}
int Mesh_object::set_index(int k_polygon,int k_index,int index)
{return connectivity.set_index(k_polygon,k_index,index);}
int Mesh_object::delete_polygon(int k_polygon)
{return connectivity.delete_polygon(k_polygon);}
int Mesh_object::delete_given_polygon(const std::vector <int> polygon)
{return connectivity.delete_given_polygon(polygon);}
int Mesh_object::flip_polygon(int k_polygon)
{return connectivity.flip_polygon(k_polygon);}
int Mesh_object::flip_polygon()
{return connectivity.flip_polygon();}




int Mesh_object::load_off_file(const char* filename)
{return load_off_file(std::string(filename));}
int Mesh_object::load_off_file(const std::string& filename)
{
  int ok=0;
  ok += destroy();
  ok += point_set.load_off_file(filename);
  ok += connectivity.load_off_file(filename);
  return ok;
}

int Mesh_object::load_g_file(const char* filename)
{return load_g_file(std::string(filename));}
int Mesh_object::load_g_file(const std::string& filename)
{
  int ok=0;
  ok += destroy();
  ok += point_set.load_g_file(filename);
  ok += connectivity.load_g_file(filename);
  return ok;
}

int Mesh_object::save_g_file(const char* filename) const
{return save_g_file(std::string(filename));}
int Mesh_object::save_g_file(const std::string& filename) const
{

  int ok=0;
  FILE *fid=NULL;
  
  fid=fopen(filename.data(),"w");
  if(fid==NULL)
    {printf("ERROR loading %s in save_g_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}

  //Just header
  fprintf(fid,"%d %d\n",point_set.size(),connectivity.polygon_number());

  if(fclose(fid)!=0)
    {printf("ERROR closing %s in write_off_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  //write geometry
  point_set.write_g_type(filename);
  
  //write connectivity
  connectivity.write_g_type(filename);



  //write end
  fid=NULL;
  fid=fopen(filename.data(),"a");
  if(fid==NULL)
    {printf("ERROR loading %s in save_off_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}

  fprintf(fid,"EOF Mesh Geometry %s [OK]\n",filename.data());

  if(fclose(fid)!=0)
    {printf("ERROR closing %s in write_g_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  return ok;

}



int Mesh_object::save_off_file(const char* filename) const
{return save_off_file(std::string(filename));}
int Mesh_object::save_off_file(const std::string& filename) const
{
  int ok=0;
  FILE *fid=NULL;
  
  fid=fopen(filename.data(),"w");
  if(fid==NULL)
    {printf("ERROR loading %s in save_off_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  //Just header

  fprintf(fid,"# Geometrical Mesh [%s] in [%s]\n",name.data(),filename.data());
  time_t current_time=time(NULL);
  fprintf(fid,"# Mesh Exported on %s",asctime(localtime(&current_time)));
  fprintf(fid,"#\n");
  fprintf(fid,"# N_vertices [%d]\n",point_set.size());
  fprintf(fid,"# N_polygons [%d]\n",connectivity.polygon_number());
  if(connectivity.edge_number()!=0)
    fprintf(fid,"# N_edges [%d]\n",connectivity.edge_number());
  fprintf(fid,"#\n");

  //OFF Header
  fprintf(fid,"OFF\n");
  fprintf(fid,"%d %d %d\n",point_set.size(),connectivity.polygon_number(),connectivity.edge_number());

  if(fclose(fid)!=0)
    {printf("ERROR closing %s in write_off_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  //write geometry
  point_set.write_off_type(filename);
  
  //write connectivity
  connectivity.write_off_type(filename);



  //write end
  fid=NULL;
  fid=fopen(filename.data(),"a");
  if(fid==NULL)
    {printf("ERROR loading %s in save_off_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}

  fprintf(fid,"EOF Mesh Geometry %s [OK]\n",filename.data());

  if(fclose(fid)!=0)
    {printf("ERROR closing %s in write_off_file in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  return ok;

}


int Mesh_object::load_collada_file(const char* filename)
{return load_collada_file(std::string(filename),std::string("_ALL_"));}
int Mesh_object::load_collada_file(const std::string& filename)
{return load_collada_file(filename,std::string("_ALL_"));}
int Mesh_object::load_collada_file(const char* filename,const char* geometry_name)
{return load_collada_file(std::string(filename),std::string(geometry_name));}
int Mesh_object::load_collada_file(const std::string& filename,const std::string& geometry_name)
{
  int ok=0;
  ok += destroy();
  ok += point_set.load_collada_file(filename,geometry_name);
  ok += connectivity.load_collada_file(filename,geometry_name);
  return ok;
}




int Mesh_object::build_normal_per_polygon()
{
  if(parameters.size()<1)
    parameters.resize(1);
  parameters[0].destroy();

  int N_polygons=polygon_number();

  // per polygon type
  parameters[0].resize_int(1);
  parameters[0].get_int(0).resize(1);
  parameters[0].get_int(0)[0]=1; //PER_POLYGON

  // fill polygon normals
  parameters[0].resize_double(1);
  parameters[0].get_double(0).resize(3*N_polygons);
  int k_polygon=0;
  V_3D x0,x1,x2;
  V_3D n;

  int k_dim=0;
  for(k_polygon=0;k_polygon<N_polygons;k_polygon++)
    {
      //the first triangle is enough
      x0=get_vertex_v3d(k_polygon,0);
      x1=get_vertex_v3d(k_polygon,1);
      x2=get_vertex_v3d(k_polygon,2);
      
      n=((x1-x0).vector_prod(x2-x0)).normalized();
      for(k_dim=0;k_dim<3;k_dim++)
        parameters[0].get_double(0)[3*k_polygon+k_dim]=n[k_dim];
    }
  
  return 0;

}


int Mesh_object::build_normal_intermediate(double angle_threshold)
{

  int ok=0;

  //convert angle into degree
  if(angle_threshold>3.14159)
    angle_threshold = angle_threshold*3.14159/180.0;

  //first build one normal per polygon
  ok += build_normal_per_polygon();

  // second int vector will give the position
  parameters[0].resize_int(2);
  parameters[0].get_int(1).resize(polygon_number());



  int k_polygon=0,N_polygon=polygon_number();
  int k_vertex=0,N_vertex=0;
  V_3D n0,n1,normal,current_normal;
  int k_ring=0,k_ring_2=0,ring_size=0;

  double_vector normal_vector;
  int_vector polygon_ring;
  int is_sharp_edge=0;
  int count=0;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      // integrated position in the vector for a direct access
      parameters[0].get_int(1)[k_polygon] = count;

      N_vertex=polygon_size(k_polygon);
      for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
        {

          //get polygon ring of this vertex
          polygon_ring=get_polygon_ring(k_polygon,k_vertex);
          ring_size=polygon_ring.size();

          //take the dot product between every normals of every polygon surrounding the current vertex
          is_sharp_edge=0;
          for(k_ring=0;(is_sharp_edge==0) && (k_ring<ring_size);k_ring++)
            {
              n0=get_polygon_normal(polygon_ring[k_ring]);
              for(k_ring_2=k_ring+1;(is_sharp_edge==0) && (k_ring_2<ring_size);k_ring_2++)
                {
                  n1=get_polygon_normal(polygon_ring[k_ring_2]);
                  
                  if( n0.angle(n1)>angle_threshold )
                    is_sharp_edge=1;
                }
            }


          // if is a sharp edge, then copy the polygon normal for the vertex
          normal = get_polygon_normal(k_polygon);
          if(is_sharp_edge==1)
            {normal_vector.add(normal[0]);normal_vector.add(normal[1]);normal_vector.add(normal[2]);}

          
          // if is a smooth edge, then copy the smooth normal for the vertex
          else
            {

              normal.set(0.0,0.0,0.0);
              for(k_ring=0;k_ring<ring_size;k_ring++) 
                {
                  current_normal = get_polygon_normal(polygon_ring[k_ring]);
                  normal += current_normal;
                }
              normal = normal.normalized();
              normal_vector.add(normal[0]);normal_vector.add(normal[1]);normal_vector.add(normal[2]);
            }


          count++;
        }
    }


  //replace the normal vector
  parameters[0].get_double(0) = normal_vector;
  parameters[0].get_int(0)[0]=3; //PER_VERTEX_AND_POLYGON

  return ok;
}


int Mesh_object::build_normal_per_vertex()
{
  
  int ok=0;

  //first build one normal per polygon
  ok += build_normal_per_polygon();

  int k_ring=0,ring_size=0;
  int_vector polygon_ring;
  V_3D normal,n;

  double_vector normal_vector;
  int k_vertex=0,N_vertex=vertex_number();
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    {
      polygon_ring = get_polygon_ring(k_vertex);
      
      ring_size=polygon_ring.size();
      normal.set(0,0,0);
      for(k_ring=0;k_ring<ring_size;k_ring++)
        {
          n = get_polygon_normal(polygon_ring[k_ring]);
          normal += n;
        }
      normal = normal.normalized();
      normal_vector.add(normal[0]);normal_vector.add(normal[1]);normal_vector.add(normal[2]);
    }


  //replace the normal vector
  parameters[0].get_double(0) = normal_vector;
  parameters[0].get_int(0)[0]=2; //PER_VERTEX



  return ok;
  
}


int Mesh_object::normal_type() const
{
  if(parameters.size()<1 || parameters[0].size_int()<1 || parameters[0].get_int(0).size()<1)
    return -1;
  else
    return parameters[0].get_int(0)[0];
}

V_3D Mesh_object::get_polygon_normal(int k_polygon) const
{
  if(normal_type()!=1 || parameters[0].get_double(0).size()!=3*polygon_number())
    {printf("Error in get_polygon_normal(%d), normal are not calculated for each polygons (type=%d)\n",k_polygon,normal_type());exit(-1);}
  
  V_3D normal;
  for(int k_dim=0;k_dim<3;k_dim++)
    normal[k_dim]=parameters[0].get_double(0)[3*k_polygon+k_dim];

  return normal;
}

V_3D Mesh_object::get_polygon_and_vertex_normal(int k_polygon,int k_vertex) const
{
  if(normal_type()!=3)
    {printf("Error in get_polygon_and_vertex_normal(%d,%d), normal are not calculated for each polygons_and_vertex (type=%d)\n",k_polygon,k_vertex,normal_type());exit(-1);}
  if(
     (parameters.size()<1) || 
     (parameters[0].size_int()<2) || 
     (parameters[0].get_int(1).size()<polygon_number())||
     (parameters[0].size_double()<1)
     )
    {printf("Error in get_polygon_and_vertex_normal(%d,%d), size are not correct\n",k_polygon,k_vertex);exit(-1);}
  

  int k_position=0;
  V_3D normal;
  k_position = parameters[0].get_int(1)[k_polygon]+k_vertex;
  for(int k_dim=0;k_dim<3;k_dim++)
    {
      if(3*k_position+2>=parameters[0].get_double(0).size())
        {printf("Error in get_polygon_and_vertex_normal(%d,%d), recored normal have size %d and not %d\n",k_polygon,k_vertex,parameters[0].get_double(0).size(),k_position);exit(-1);}
      normal[k_dim]=parameters[0].get_double(0)[3*k_position+k_dim];
    }
  return normal;
}
V_3D Mesh_object::get_vertex_normal(int k_vertex) const
{
  if(normal_type()!=2 || parameters[0].get_double(0).size()!=3*vertex_number())
    {printf("Error in get_vertex_normal(%d), normal are not calculated for each vertex (type=%d)\n",k_vertex,normal_type());exit(-1);}
 

  V_3D normal;
  for(int k_dim=0;k_dim<3;k_dim++)
    normal[k_dim]=parameters[0].get_double(0)[3*k_vertex+k_dim];

  return normal;
}
int Mesh_object::set_vertex_normal(int k_vertex,const V_3D& n)
{
  if(normal_type()!=2 || parameters[0].get_double(0).size()!=3*vertex_number())
    {printf("Error in get_polygon_normal(%d), normal are not calculated for each polygons (type=%d)\n",k_vertex,normal_type());exit(-1);}
 

  for(int k_dim=0;k_dim<3;k_dim++)
    parameters[0].get_double(0)[3*k_vertex+k_dim]=n[k_dim];

  return 0;
}

V_3D Mesh_object::get_vertex_normal(int k_polygon,int k_vertex) const
{int k_index=get_index(k_polygon,k_vertex);return get_vertex_normal(k_index);}
int Mesh_object::set_vertex_normal(int k_polygon,int k_vertex,const V_3D& n)
{int k_index=get_index(k_polygon,k_vertex);return set_vertex_normal(k_index,n);}



int_vector Mesh_object::get_one_ring(int k_polygon,int k_vertex) const
{return connectivity.get_one_ring(k_polygon,k_vertex);}
int_vector Mesh_object::get_one_ring(int k_vertex) const
{return connectivity.get_one_ring(k_vertex);}

int_vector Mesh_object::get_polygon_ring(int k_polygon,int k_vertex) const
{return connectivity.get_polygon_ring(k_polygon,k_vertex);}
int_vector Mesh_object::get_polygon_ring(int k_vertex) const
{return connectivity.get_polygon_ring(k_vertex);}

int Mesh_object::load_obj_file(const char* filename)
{std::string file=filename;return load_obj_file(file);}
int Mesh_object::load_obj_file(const std::string& filename)
{
  int ok=0;
  ok += destroy();
  ok += point_set.load_obj_file(filename);
  ok += connectivity.load_obj_file(filename);

  return ok;
}

int Mesh_object::save_obj_file(const char* filename)
{std::string file=filename;return save_obj_file(file);}
int Mesh_object::save_obj_file(const std::string& filename)
{
  int ok=0;
  FILE *fid=NULL;

  //Just the header
  fid=fopen(filename.data(),"w");
  if(fid==NULL)
    {printf("ERROR loading %s in save_obj_type in Mesh [%s]\n",filename.data(),name.data());exit(-1);}

  fprintf(fid,"# Geometrical Mesh [%s] in [%s]\n",name.data(),filename.data());
  time_t current_time=time(NULL);
  fprintf(fid,"# Mesh Exported on %s",asctime(localtime(&current_time)));
  fprintf(fid,"#\n");
  fprintf(fid,"# N_vertices [%d]\n",point_set.size());
  fprintf(fid,"# N_polygons [%d]\n",connectivity.polygon_number());
  if(connectivity.edge_number()!=0)
    fprintf(fid,"# N_edges [%d]\n",connectivity.edge_number());
  fprintf(fid,"#\n");
  
  if(fclose(fid)!=0)
    {printf("ERROR closing %s in write_off_type in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  //write geometry
  point_set.write_obj_type(filename);

  //write connectivity
  connectivity.write_obj_type(filename);
  
  //write end
  fid=NULL;
  fid=fopen(filename.data(),"a");
  if(fid==NULL)
    {printf("ERROR loading %s in save_obj_type in Mesh [%s]\n",filename.data(),name.data());exit(-1);}

  fprintf(fid,"EOF Mesh Geometry %s [OK]\n",filename.data());

  if(fclose(fid)!=0)
    {printf("ERROR closing %s in write_obj_type in Mesh [%s]\n",filename.data(),name.data());exit(-1);}


  return ok;
}




int Mesh_object::matrix_transform(const Matrix& M) 
{return point_set.matrix_transform(M);}

int Mesh_object::rotate(const V_3D& n,double angle)
{return point_set.rotate(n,angle);}
int Mesh_object::rotate(double n_x,double n_y,double n_z,double angle)
{return point_set.rotate(n_x,n_y,n_z,angle);}
int Mesh_object::rotate(const double *n,double angle)
{return point_set.rotate(n,angle);}


int Mesh_object::rotate_centered(const V_3D& n,double angle)
{return point_set.rotate_centered(n,angle);}
int Mesh_object::rotate_centered(double n_x,double n_y,double n_z,double angle)
{return point_set.rotate_centered(n_x,n_y,n_z,angle);}
int Mesh_object::rotate_centered(const double *n,double angle)
{return point_set.rotate_centered(n,angle);}

int Mesh_object::rotate(const V_3D& n,double angle,const V_3D& center)
{return point_set.rotate(n,angle,center);}
int Mesh_object::rotate(double n_x,double n_y,double n_z,double angle,const V_3D& center)
{return point_set.rotate(n_x,n_y,n_z,angle,center);}
int Mesh_object::rotate(const double *n,double angle,const V_3D& center)
{return point_set.rotate(n,angle,center);}

int Mesh_object::translate(double t_x,double t_y,double t_z)
{return point_set.translate(t_x,t_y,t_z);}
int Mesh_object::translate(const double *t)
{return point_set.translate(t);}
int Mesh_object::translate(const V_3D& t)
{return point_set.translate(t);}

int Mesh_object::scale(double s_x,double s_y,double s_z)
{return point_set.scale(s_x,s_y,s_z);}
int Mesh_object::scale(const double *s)
{return point_set.scale(s);}
int Mesh_object::scale(const V_3D& s)
{return point_set.scale(s);}
int Mesh_object::scale(double s)
{return point_set.scale(s);}

int Mesh_object::scale_centered(double s_x,double s_y,double s_z)
{return point_set.scale_centered(s_x,s_y,s_z);}
int Mesh_object::scale_centered(const double *s)
{return point_set.scale_centered(s);}
int Mesh_object::scale_centered(const V_3D& s)
{return point_set.scale_centered(s);}
int Mesh_object::scale_centered(double s)
{return point_set.scale_centered(s);}

double Mesh_object::get_closest_distance(double x,double y,double z,int *index) const
{return point_set.get_closest_distance(x,y,z,index);}
double Mesh_object::get_closest_distance(const double* x,int *index) const
{return point_set.get_closest_distance(x,index);}
double Mesh_object::get_closest_distance(const V_3D& x,int *index) const
{return point_set.get_closest_distance(x,index);}

double Mesh_object::get_further_distance(double x,double y,double z,int *index) const
{return point_set.get_further_distance(x,y,z,index);}
double Mesh_object::get_further_distance(const double* x,int *index) const
{return point_set.get_further_distance(x,index);}
double Mesh_object::get_further_distance(const V_3D& x,int *index) const
{return point_set.get_further_distance(x,index);}

double Mesh_object::distance_to_index(double x,double y,double z,int index) const
{return point_set.distance_to_index(x,y,z,index);}
double Mesh_object::distance_to_index(const double* x,int index) const
{return point_set.distance_to_index(x,index);}
double Mesh_object::distance_to_index(const V_3D& x,int index) const
{return point_set.distance_to_index(x,index);}

std::vector <int> Mesh_object::N_closest(const V_3D& center,int _N_closest) const
{return point_set.N_closest(center,_N_closest);}
std::vector <int> Mesh_object::N_closest(double x,double y,double z,int _N_closest) const
{return point_set.N_closest(x,y,z,_N_closest);}
std::vector <int> Mesh_object::N_closest(double *x,int _N_closest) const
{return point_set.N_closest(x,_N_closest);}


Polygon Mesh_object::get_polygon(int k_polygon) const
{
  if(k_polygon<0 || k_polygon>polygon_number())
    {printf("Error in get_polygon(%d) in Mesh_object, size=%d\n",k_polygon,polygon_number());exit(-1);}

  std::vector <int> poly_index=get_index_of_polygon(k_polygon);
  std::vector <V_3D> poly_geom;
  for(int k=0;k<int(poly_index.size());k++)
    poly_geom.push_back(get_vertex_v3d(poly_index[k]));
  
  Polygon P(poly_geom);
  return P;
}

V_3D Mesh_object::closest_mesh_point(const V_3D& x) const
{

  // check the distance to every polygon
  int k_polygon=0,N_polygon=polygon_number();
  double current_dist=0.0,min_dist=9999.9;
  V_3D closest_point;

  Polygon current_polygon;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      current_polygon=get_polygon(k_polygon);
      current_dist=current_polygon.shortest_distance_to_point(x);

      if(current_dist<min_dist)
        {
          min_dist=current_dist;
          closest_point=current_polygon.closest_point(x);
        }
    }
  return closest_point;
}



std::vector <Curve_3D> Mesh_object::plane_intersection(const V_3D& n,const V_3D& x0) const
{
  std::vector <Curve_3D> curves;
  int k_polygon=0,N_polygon=polygon_number();

  std::vector <Segment> s;int type=-1;
  std::vector <V_3D> temp;
  //check every intersection for every polygons
  //store the intersection in an unordered vector of Segment
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      //check the intersection
      temp = get_polygon(k_polygon).plane_intersection(n,x0,&type);
      if(type==1 || type==3)
        s.push_back(Segment(temp[0],temp[1]));
    }

  int segment_number=s.size();


  // now order the segments in curve
  std::vector <int> is_segment_added;
  is_segment_added.resize(segment_number);
  
  Curve_3D current_curve;
  current_curve.add(s[0][0]);
  current_curve.add(s[0][1]);
  is_segment_added[0]=1;
  
  int k_segment=-1,is_addition=1;
  while(is_addition==1)//as long as there is segment to add
    {
      //check if there is at least one added segment during the whole pass
      is_addition=0;

      //add forward
      for(k_segment=0;k_segment<segment_number;k_segment++)
        if(is_segment_added[k_segment]==0)
          {
            if(s[k_segment][0]==current_curve[current_curve.size()-1])
              {current_curve.add(s[k_segment][1]);is_segment_added[k_segment]=1;is_addition=1;}
            else if(s[k_segment][1]==current_curve[current_curve.size()-1])
              {current_curve.add(s[k_segment][0]);is_segment_added[k_segment]=1;is_addition=1;}
          }
      //add back
      for(k_segment=0;k_segment<segment_number;k_segment++)
        if(is_segment_added[k_segment]==0)
          {
            if(s[k_segment][0]==current_curve[0])
              {current_curve.add_first(s[k_segment][1]);is_segment_added[k_segment]=1;is_addition=1;}
            else if(s[k_segment][1]==current_curve[0])
              {current_curve.add_first(s[k_segment][0]);is_segment_added[k_segment]=1;is_addition=1;}
          }


      //if there is no added segment, it might be non-connected cases
      if(is_addition==0)
        {
          // add this curve to the vector
          curves.push_back(current_curve);
          current_curve.destroy();

          //check if every segment has been added
          for(k_segment=0;is_addition==0 && k_segment<segment_number;k_segment++)
            if(is_segment_added[k_segment]==0)//Then add it in the next curve
              {
                is_addition=1;
                current_curve.add(s[k_segment][0]);current_curve.add(s[k_segment][1]);
                is_segment_added[k_segment]=1;
              }
        }
    }
  

  return curves;
}


std::vector <int> Mesh_object::add_unique_polygon(const Polygon& polygon)
{

  int k=0;
  int N_polygon=polygon.size();
  int index=0;
  std::vector <int> index_polygon;
  
  for(k=0;k<N_polygon;k++)
    {
      index = exist(polygon[k],0.00001);
      if(index==-1) //new polygon
        {
          add_vertex(polygon[k]);
          index_polygon.push_back(get_vertex_number()-1);
        }
      else //already exists
        index_polygon.push_back(index);
    }

  add_polygon(index_polygon);
  return index_polygon;

}


Mesh_object Mesh_object::half_space_intersection(const V_3D& n,const V_3D& x0)
{int type=0;return half_space_intersection(n,x0,&type);}
Mesh_object Mesh_object::half_space_intersection(const V_3D& n,const V_3D& x0,int *type)
{
  Mesh_object new_mesh;
  int k_polygon=0;
  int N_polygon=polygon_number();

  *type = 0;
  Polygon p; int t=0;
  Polygon q;
  //cut every polygon
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      p = get_polygon(k_polygon);
      q = p.half_space_intersection(n,x0,&t);

      if(t!=2)//add the polygon and build the connectivity
        new_mesh.add_unique_polygon(q);
      if(t!=0)//at least one polygon is cutted
        *type=1;
    }


  if(new_mesh.vertex_number()==0)//no more mesh at all
    *type=2;
  
  new_mesh.build_connectivity_index();

  return new_mesh;
}

int Mesh_object::build_connectivity_index()
{return connectivity.build_all_index();}


int Mesh_object::color_type() const
{
  if(parameters.size()<2 || parameters[1].size_int()<1 || parameters[1].size_int(0)<1)
    return -1;
  else
    return parameters[1].get_int(0)[0];
}

int Mesh_object::set_per_vertex_color()
{
  if(parameters.size()<2)
    parameters.resize(2);

  // set double size
  if(parameters[1].size_double()<1)
    parameters[1].resize_double(1);
  if(parameters[1].size_double(0)!=4*vertex_number())
    parameters[1].resize_double(0,4*vertex_number());

  // set per polygon type
  if(parameters[1].size_int()<1)
    parameters[1].resize_int(1);
  if(parameters[1].size_int(0)<1)
    parameters[1].resize_int(0,1);

  parameters[1].get_int(0)[0]=1;


  return 0;
}
int Mesh_object::set_per_polygon_color()
{
  if(parameters.size()<2)
    parameters.resize(2);

  // set double size
  if(parameters[1].size_double()<1)
    parameters[1].resize_double(1);
  if(parameters[1].size_double(0)!=4*polygon_number())
    parameters[1].resize_double(0,4*polygon_number());

  // set per polygon type
  if(parameters[1].size_int()<1)
    parameters[1].resize_int(1);
  if(parameters[1].size_int(0)<1)
    parameters[1].resize_int(0,1);

  parameters[1].get_int(0)[0]=2;


  return 0;
}


int Mesh_object::fill_color(double r,double g,double b,double t)
{
  if(color_type()!=1 || parameters[1].size_double(0)!=4*vertex_number())
    set_per_vertex_color();

  int k_vertex=0,N_vertex=vertex_number();
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    {
      parameters[1].get_double(0)[4*k_vertex+0]=r;
      parameters[1].get_double(0)[4*k_vertex+1]=g;
      parameters[1].get_double(0)[4*k_vertex+2]=b;
      parameters[1].get_double(0)[4*k_vertex+3]=t;

    }
  return 0;
}
int Mesh_object::set_per_vertex_color(int k_vertex,double r,double g,double b,double t)
{
  if(color_type()!=1 || parameters[1].size_double(0)!=4*vertex_number())
    set_per_vertex_color();

  parameters[1].get_double(0)[4*k_vertex+0]=r;
  parameters[1].get_double(0)[4*k_vertex+1]=g;
  parameters[1].get_double(0)[4*k_vertex+2]=b;
  parameters[1].get_double(0)[4*k_vertex+3]=t;

  return 0;
}

int Mesh_object::get_per_polygon_color(int k_polygon,double *r,double *g,double *b,double *t) const
{
  if(color_type()!=2 || parameters[1].size_double(0)!=4*polygon_number())
    return -1;
  else
    {                                           
      *r = parameters[1].get_double(0)[4*k_polygon+0];
      *g = parameters[1].get_double(0)[4*k_polygon+1];
      *b = parameters[1].get_double(0)[4*k_polygon+2];
      *t = parameters[1].get_double(0)[4*k_polygon+3];
    }
  return 0;

}

int Mesh_object::set_per_polygon_color(int k_polygon,double r,double g,double b,double t)
{
  if(color_type()!=2 || parameters[1].size_double(0)!=4*polygon_number())
    set_per_polygon_color();

  parameters[1].get_double(0)[4*k_polygon+0]=r;
  parameters[1].get_double(0)[4*k_polygon+1]=g;
  parameters[1].get_double(0)[4*k_polygon+2]=b;
  parameters[1].get_double(0)[4*k_polygon+3]=t;

  return 0;
}

int Mesh_object::get_per_vertex_color(int k_vertex,double *r,double *g,double *b,double *t) const
{
  if(color_type()!=1 || parameters[1].size_double(0)!=4*vertex_number())
    return -1;
  else
    {                                           
      *r = parameters[1].get_double(0)[4*k_vertex+0];
      *g = parameters[1].get_double(0)[4*k_vertex+1];
      *b = parameters[1].get_double(0)[4*k_vertex+2];
      *t = parameters[1].get_double(0)[4*k_vertex+3];
    }
  return 0;
}

int Mesh_object::get_per_vertex_color(int k_polygon,int k_vertex,double *r,double *g,double *b,double *t) const
{int index = get_index(k_polygon,k_vertex); return get_per_vertex_color(index,r,g,b,t);}
int Mesh_object::set_per_vertex_color(int k_polygon,int k_vertex,double r,double g,double b,double t)
{int index = get_index(k_polygon,k_vertex); return set_per_vertex_color(index,r,g,b,t);}



int Mesh_object::set_per_polygon_and_vertex_color()
{
  if(parameters.size()<2)
    parameters.resize(2);

  // set double size
  if(parameters[1].size_double()<1)
    parameters[1].resize_double(1);
  if(parameters[1].size_double(0)!=4*polygon_number()*vertex_number())
    parameters[1].resize_double(0,4*polygon_number()*vertex_number());

  // set per polygon type
  if(parameters[1].size_int()<2)
    parameters[1].resize_int(2);
  if(parameters[1].size_int(0)<1)
    parameters[1].resize_int(0,1);

  parameters[1].get_int(0)[0]=3;


  if(parameters[1].size_int(1)<polygon_number())
    parameters[1].resize_int(1,polygon_number());

  //now fill the access index
  int count=0;int N_polygon=polygon_number();
  for(int k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      parameters[1].get_int(1)[k_polygon]=count;
      count+=polygon_size(k_polygon);
    }


  return 0;
}

int Mesh_object::set_per_polygon_and_vertex_color(int k_polygon,int k_vertex,double r,double g,double b,double t)
{
  if(color_type()!=3 || parameters[1].get_int(1).size()!=polygon_number() || parameters[1].get_double(0).size()!=4*vertex_number()*polygon_number())
    {set_per_polygon_and_vertex_color();}


  int k_position = parameters[1].get_int(1)[k_polygon]+k_vertex;
  if(4*k_position+3>=parameters[1].get_double(0).size())
    {printf("Error in set_per_polygon_and_vertex_color(%d,%d,%f,%f,%f,%f), size are not ok %d-%d\n",k_polygon,k_vertex,r,g,b,t,4*k_position+3,parameters[1].get_double(0).size());exit(-1);}

  parameters[1].get_double(0,4*k_position+0)=r;
  parameters[1].get_double(0,4*k_position+1)=g;
  parameters[1].get_double(0,4*k_position+2)=b;
  parameters[1].get_double(0,4*k_position+3)=t;



  return 0;
}

int Mesh_object::get_per_polygon_and_vertex_color(int k_polygon,int k_vertex,double *r,double *g,double *b,double *t) const
{
  if(color_type()!=3 || parameters[1].get_int(1).size()!=polygon_number() || parameters[1].get_double(0).size()!=4*vertex_number()*polygon_number())
    {printf("Error in set_per_polygon_and_vertex_color(%d,%d,...), size is not correct\n",k_polygon,k_vertex);exit(-1);}

  int k_position = parameters[1].get_int(1)[k_polygon]+k_vertex;
  if(4*k_position+3>=parameters[1].get_double(0).size())
    {printf("Error in set_per_polygon_and_vertex_color(%d,%d,...), size is not ok %d-%d\n",k_polygon,k_vertex,4*k_position+3,parameters[1].get_double(0).size());exit(-1);}

  *r = parameters[1].get_double(0,4*k_position+0);
  *g = parameters[1].get_double(0,4*k_position+1);
  *b = parameters[1].get_double(0,4*k_position+2);
  *t = parameters[1].get_double(0,4*k_position+3);

  return 0;
}

V_3D Mesh_object::get_middle_point() const {return point_set.get_middle_point();}
V_3D Mesh_object::get_min_point() const {return point_set.get_min_point();}
V_3D Mesh_object::get_max_point() const {return point_set.get_max_point();}
V_3D Mesh_object::get_size_mesh() const {return point_set.get_size_set();}
int Mesh_object::center_mesh() {return point_set.center_point_set();}
int Mesh_object::unitize_size() {return point_set.unitize_size();}

int Mesh_object::set_name(const char* _name){name=_name;return 0;}
int Mesh_object::set_name(const std::string& _name){name=_name;return 0;}
std::string  Mesh_object::get_name() const{return name;}
std::string& Mesh_object::get_name()      {return name;}

int Mesh_object::build_cube()
{
  V_3D x0(0,0,0),x1(1,0,0),x2(1,1,0),x3(0,1,0),x4(0,0,1),x5(1,0,1),x6(1,1,1),x7(0,1,1);
  destroy();

  //add the vertices
  add_vertex(x0);add_vertex(x1);add_vertex(x2);add_vertex(x3);
  add_vertex(x4);add_vertex(x5);add_vertex(x6);add_vertex(x7);

  // build the connectivity
  add_polygon(0,1,2,3);
  add_polygon(0,3,7,4);
  add_polygon(4,7,6,2);
  add_polygon(2,1,5,6);
  add_polygon(0,4,5,1);
  add_polygon(3,2,6,7);

  // build the index
  build_connectivity_index();


  return 0;

}

int Mesh_object::subdivide_mid_edge()
{
  // create a new temp mesh
  Mesh_object mesh2;

  // just to save the order of the vertices
  int N_vertex=vertex_number();
  for(int k=0;k<N_vertex;k++)
    mesh2.add_vertex(get_vertex_v3d(k));

  // now get every new polygons
  std::vector <Polygon> new_polygons;
  int N_polygon=polygon_number();
  int k=0,k2=0;int poly_size=0;
  for(k=0;k<N_polygon;k++)
    {
      new_polygons = get_polygon(k).subdivide_mid_edge();
      poly_size=new_polygons.size();

      for(k2=0;k2<poly_size;k2++)
        mesh2.add_unique_polygon(new_polygons[k2]);
    }

  //exchange the meshes
  (*this)=mesh2;

  return build_connectivity_index();
}

int Mesh_object::subdivide_barycenter_mid_edge()
{
  // create a new temp mesh
  Mesh_object mesh2;

  // just to save the order of the vertices
  int N_vertex=vertex_number();
  for(int k=0;k<N_vertex;k++)
    mesh2.add_vertex(get_vertex_v3d(k));

  // now get every new polygons
  std::vector <Polygon> new_polygons;
  int N_polygon=polygon_number();
  int k=0,k2=0;int poly_size=0;
  for(k=0;k<N_polygon;k++)
    {
      new_polygons = get_polygon(k).subdivide_barycenter_mid_edge();
      poly_size=new_polygons.size();

      for(k2=0;k2<poly_size;k2++)
        mesh2.add_unique_polygon(new_polygons[k2]);
    }

  //exchange the meshes
  (*this)=mesh2;

  return build_connectivity_index();
}

int Mesh_object::subdivide_mixed_mid_edge()
{
  // create a new temp mesh
  Mesh_object mesh2;

  // just to save the order of the vertices
  int N_vertex=vertex_number();
  for(int k=0;k<N_vertex;k++)
    mesh2.add_vertex(get_vertex_v3d(k));

  // now get every new polygons
  std::vector <Polygon> new_polygons;
  int N_polygon=polygon_number();
  int k=0,k2=0;int poly_size=0;
  Polygon current;
  for(k=0;k<N_polygon;k++)
    {
      current = get_polygon(k);
      if(current.size()==3) //if a triangle
        new_polygons = current.subdivide_mid_edge();
      else // not a triangle => barycenter
        new_polygons = current.subdivide_barycenter_mid_edge();
      poly_size=new_polygons.size();

      for(k2=0;k2<poly_size;k2++)
        mesh2.add_unique_polygon(new_polygons[k2]);
    }

  //exchange the meshes
  (*this)=mesh2;

  return build_connectivity_index();  
}


int Mesh_object::fast_subdivide_mid_edge()
{

  //the subdivided Mesh
  Mesh_object mesh2;
  
  std::vector <int_vector> index_mid_point; //the index of the mid point
  std::vector <int_vector> already_existing_mid_point; //the already existing edges
  index_mid_point.resize(vertex_number());
  already_existing_mid_point.resize(vertex_number());

  int k_polygon=0;
  int N_polygon=polygon_number();
  int size_polygon=0;
  int k_edge=0;

  int index_v0=0,index_v1=0;
  int k_position_index=0;


  //first add the old vertices
  int k_vertex=0,N_vertex=vertex_number();
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    mesh2.add_vertex(get_vertex_v3d(k_vertex));


  // record the mid_points first
  int current_vertex=vertex_number();
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      size_polygon=polygon_size(k_polygon);
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0=get_index(k_polygon,k_edge);
          index_v1=get_index(k_polygon,(k_edge+1)%size_polygon);

          k_position_index=index_mid_point[index_v0].exists(index_v1);
          if( k_position_index==-1 )//edge does not exists yet
            {
              // add the order of the positioning of the mid_points
              index_mid_point[index_v0].add(index_v1);
              index_mid_point[index_v1].add(index_v0);

              // add the index in the vector
              already_existing_mid_point[index_v0].add(current_vertex);
              already_existing_mid_point[index_v1].add(current_vertex);

              // the grometrical mid_point
              mesh2.add_vertex(0.5*(get_vertex_v3d(index_v0)+get_vertex_v3d(index_v1)));

              current_vertex++;
            }
        }
    }

  int k_position_index2=0,index_mid2=0;
  int index_vm1=0;
  int index_mid=0;
  std::vector <int> temp_poly;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      // the triangles/mid_edges-barycenter
      size_polygon=polygon_size(k_polygon);
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0  = get_index(k_polygon,k_edge);
          index_v1  = get_index(k_polygon,(k_edge+1)%size_polygon);
          index_vm1 = get_index(k_polygon,(k_edge-1>=0?k_edge-1:size_polygon-1));

          k_position_index  = index_mid_point[index_v0].exists(index_v1);
          k_position_index2 = index_mid_point[index_v0].exists(index_vm1);

          index_mid  = already_existing_mid_point[index_v0][k_position_index];
          index_mid2 = already_existing_mid_point[index_v0][k_position_index2];


          mesh2.add_polygon(index_v0,index_mid,index_mid2);
        }

      //add the non triangular part linking every barycenters
      temp_poly.resize(0);
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0  = get_index(k_polygon,k_edge);
          index_v1  = get_index(k_polygon,(k_edge+1)%size_polygon);

          k_position_index = index_mid_point[index_v0].exists(index_v1);
          index_mid=already_existing_mid_point[index_v0][k_position_index];
          temp_poly.push_back(index_mid);
        }
      mesh2.add_polygon(temp_poly);
    }
  


  mesh2.build_connectivity_index();
  (*this)=mesh2;
  return 0;
  
}


int Mesh_object::fast_subdivide_barycenter_mid_edge()
{


  //the subdivided Mesh
  Mesh_object mesh2;
  
  std::vector <int_vector> index_mid_point; //the index of the mid point
  std::vector <int_vector> already_existing_mid_point; //the already existing edges
  index_mid_point.resize(vertex_number());
  already_existing_mid_point.resize(vertex_number());

  //the index of the barycenter
  int_vector barycenter_index;
  barycenter_index.resize(polygon_number());

  int k_polygon=0;
  int N_polygon=polygon_number();
  int size_polygon=0;
  int k_edge=0;

  int index_v0=0,index_v1=0;
  int k_position_index=0;


  //first add the old vertices
  int k_vertex=0,N_vertex=vertex_number();
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    mesh2.add_vertex(get_vertex_v3d(k_vertex));


  // record the mid_points first
  int current_vertex=vertex_number();
  V_3D temp_barycenter;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      size_polygon=polygon_size(k_polygon);
      temp_barycenter.set(0,0,0);
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0=get_index(k_polygon,k_edge);
          index_v1=get_index(k_polygon,(k_edge+1)%size_polygon);

          k_position_index=index_mid_point[index_v0].exists(index_v1);
          if( k_position_index==-1 )//edge does not exists yet
            {
              // add the order of the positioning of the mid_points
              index_mid_point[index_v0].add(index_v1);
              index_mid_point[index_v1].add(index_v0);

              // add the index in the vector
              already_existing_mid_point[index_v0].add(current_vertex);
              already_existing_mid_point[index_v1].add(current_vertex);

              // the geometrical mid_point
              mesh2.add_vertex(0.5*(get_vertex_v3d(index_v0)+get_vertex_v3d(index_v1)));

              current_vertex++;
            }
          temp_barycenter+=get_vertex_v3d(index_v0);
        }

      temp_barycenter/=double(size_polygon);
      barycenter_index[k_polygon]=current_vertex;
      mesh2.add_vertex(temp_barycenter);
      current_vertex++;
    }






  //2. Fill the new mesh with the connectivity
  int k_position_index2=0,index_mid2=0;
  int index_vm1=0;
  int index_mid=0;
  int index_barycenter=0;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      // the triangles/mid_edges-barycenter
      size_polygon=polygon_size(k_polygon);
      index_barycenter = barycenter_index[k_polygon];
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0  = get_index(k_polygon,k_edge);
          index_v1  = get_index(k_polygon,(k_edge+1)%size_polygon);
          index_vm1 = get_index(k_polygon,(k_edge-1>=0?k_edge-1:size_polygon-1));

          k_position_index  = index_mid_point[index_v0].exists(index_v1);
          k_position_index2 = index_mid_point[index_v0].exists(index_vm1);

          index_mid  = already_existing_mid_point[index_v0][k_position_index];
          index_mid2 = already_existing_mid_point[index_v0][k_position_index2];

          mesh2.add_polygon(index_v0,index_mid,index_barycenter,index_mid2);
        }
    }



  mesh2.build_connectivity_index();
  (*this)=mesh2;
  return 0;



}




int Mesh_object::fast_subdivide_mixed_mid_edge()
{


  //the subdivided Mesh
  Mesh_object mesh2;
  
  std::vector <int_vector> index_mid_point; //the index of the mid point
  std::vector <int_vector> already_existing_mid_point; //the already existing edges
  index_mid_point.resize(vertex_number());
  already_existing_mid_point.resize(vertex_number());

  //the index of the barycenter
  int_vector barycenter_index;
  barycenter_index.resize(polygon_number());

  int k_polygon=0;
  int N_polygon=polygon_number();
  int size_polygon=0;
  int k_edge=0;

  int index_v0=0,index_v1=0;
  int k_position_index=0;


  //first add the old vertices
  int k_vertex=0,N_vertex=vertex_number();
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    mesh2.add_vertex(get_vertex_v3d(k_vertex));


  // record the mid_points first
  int current_vertex=vertex_number();
  V_3D temp_barycenter;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      size_polygon=polygon_size(k_polygon);
      temp_barycenter.set(0,0,0);
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0=get_index(k_polygon,k_edge);
          index_v1=get_index(k_polygon,(k_edge+1)%size_polygon);

          k_position_index=index_mid_point[index_v0].exists(index_v1);
          if( k_position_index==-1 )//edge does not exists yet
            {
              // add the order of the positioning of the mid_points
              index_mid_point[index_v0].add(index_v1);
              index_mid_point[index_v1].add(index_v0);

              // add the index in the vector
              already_existing_mid_point[index_v0].add(current_vertex);
              already_existing_mid_point[index_v1].add(current_vertex);

              // the geometrical mid_point
              mesh2.add_vertex(0.5*(get_vertex_v3d(index_v0)+get_vertex_v3d(index_v1)));

              current_vertex++;
            }

          temp_barycenter+=get_vertex_v3d(index_v0);
        }

      temp_barycenter/=double(size_polygon);
      barycenter_index[k_polygon]=current_vertex;
      if(size_polygon!=3)
        {
          mesh2.add_vertex(temp_barycenter);
          current_vertex++;
        }
    }






  //2. Fill the new mesh with the connectivity
  int k_position_index2=0,index_mid2=0;
  int index_vm1=0;
  int index_mid=0;
  int index_barycenter=0;
  std::vector <int> temp_poly;
  for(k_polygon=0;k_polygon<N_polygon;k_polygon++)
    {
      // the triangles/mid_edges-barycenter
      size_polygon=polygon_size(k_polygon);
      index_barycenter = barycenter_index[k_polygon];
      for(k_edge=0;k_edge<size_polygon;k_edge++)
        {
          index_v0  = get_index(k_polygon,k_edge);
          index_v1  = get_index(k_polygon,(k_edge+1)%size_polygon);
          index_vm1 = get_index(k_polygon,(k_edge-1>=0?k_edge-1:size_polygon-1));

          k_position_index  = index_mid_point[index_v0].exists(index_v1);
          k_position_index2 = index_mid_point[index_v0].exists(index_vm1);

          index_mid  = already_existing_mid_point[index_v0][k_position_index];
          index_mid2 = already_existing_mid_point[index_v0][k_position_index2];

          if(size_polygon!=3)
            mesh2.add_polygon(index_v0,index_mid,index_barycenter,index_mid2);
          else
            mesh2.add_polygon(index_v0,index_mid,index_mid2);
        }


      //add the non triangular part linking every barycenters
      if(size_polygon==3)
        {
          temp_poly.resize(0);
          for(k_edge=0;k_edge<size_polygon;k_edge++)
            {
              index_v0  = get_index(k_polygon,k_edge);
              index_v1  = get_index(k_polygon,(k_edge+1)%size_polygon);
              
              k_position_index = index_mid_point[index_v0].exists(index_v1);
              index_mid=already_existing_mid_point[index_v0][k_position_index];
              temp_poly.push_back(index_mid);
            }
          mesh2.add_polygon(temp_poly);
        }

    }



  mesh2.build_connectivity_index();
  (*this)=mesh2;
  return 0;



}


int Mesh_object::Laplacien_smooth(double lambda)
{
  int k_vertex=0;
  int N_vertex=vertex_number();

  int_vector one_ring;

  int k_one_ring=0;
  int N_one_ring=0;
  V_3D barycenter;
  V_3D x,u;

  std::vector <V_3D> temp;
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    {
      //get barycenter of the one ring
      barycenter.set(0,0,0);
      one_ring = get_one_ring(k_vertex);
      N_one_ring=one_ring.size();
      for(k_one_ring=0;k_one_ring<N_one_ring;k_one_ring++)
        barycenter+=get_vertex_v3d(one_ring[k_one_ring]);
      barycenter/=N_one_ring;

      // set x = x+lambda*Delta(x)
      x = get_vertex_v3d(k_vertex);
      u = barycenter-x;
      temp.push_back(x+lambda*u);
    }
  for(k_vertex=0;k_vertex<N_vertex;k_vertex++)
    set_vertex(k_vertex,temp[k_vertex]);



  return 0;
}


int Mesh_object::triangulate()
{return connectivity.triangulate();}




int Mesh_object::set_intermediate_texture()
{
  if(parameters.size()<3)
    parameters.resize(3);

  // set double size
  if(parameters[2].size_double()<1)
    parameters[2].resize_double(1);
  parameters[2].resize_double(0,0);

  // set per polygon type
  if(parameters[2].size_int()<1)
    parameters[2].resize_int(3);
  if(parameters[2].size_int(0)<1)
    parameters[2].resize_int(0,1);
  parameters[2].resize_int(1,1);
  parameters[2].resize_int(2,0);

  parameters[2].get_int(0)[0]=0;
  parameters[2].get_int(1)[0]=-1;



  return 0;
}


int Mesh_object::texture_type() const
{
  if(parameters.size()<3 || parameters[2].size_int()<1 || parameters[2].size_int(0)<1)
    return -1;
  else
    return parameters[2].get_int(0)[0];
}

int Mesh_object::set_intermediate_texture(int k_polygon,int k_vertex,double t0,double t1)
{
  if(k_vertex>=vertex_number() || k_polygon>=polygon_number())
    {printf("Error in set_texture, k_vertex(%d),k_polygon(%d) is too large(%d)(%d)\n",k_vertex,k_polygon,vertex_number(),polygon_number());exit(-1);}

  if(texture_type()!=0)
    set_intermediate_texture();
  
  int k_polygon_access=get_polygon_access(k_polygon);
  int k_index=k_polygon_access+k_vertex;
  if(k_index<0 || k_index>parameters[2].get_int(2).size())
    {printf("Error in set_intermediate_texture, k_index(%d) is wrong size(%d)\n",k_index,parameters[2].get_int(2).size());exit(-1);}

  int k_double_index = parameters[2].get_int(2)[k_index];

  parameters[2].get_double(0)[2*k_double_index+0]=t0;
  parameters[2].get_double(0)[2*k_double_index+1]=t1;
  return 0;
}


int Mesh_object::get_intermediate_texture(int k_polygon,int k_vertex,double* t0,double* t1) const
{
  if(texture_type()!=0)
    return -1;
  if(k_vertex>=vertex_number() || k_polygon>=polygon_number())
    return -1;

  int k_polygon_access=get_polygon_access(k_polygon);
  int k_index=k_polygon_access+k_vertex;
  if(k_index<0 || k_index>parameters[2].get_int(2).size())
    {printf("Error in get_intermediate_texture, k_index(%d) is wrong size(%d)\n",k_index,parameters[2].get_int(2).size());exit(-1);}
  int k_double_index = parameters[2].get_int(2)[k_index];


  if(2*k_double_index+0<0 || 2*k_double_index+1>parameters[2].get_double(0).size())
    {printf("Error in get_intermediate_texture in Mesh_object, texture_polygon %d in vertex %d has k_polyon_access %d, index %d, and array of texture has a size of %d\n",k_polygon,k_vertex,k_polygon_access,k_index,parameters[2].get_double(0).size()/2);exit(-1);}

  *t0=parameters[2].get_double(0)[2*k_double_index+0];
  *t1=parameters[2].get_double(0)[2*k_double_index+1];
  return 0;
  
}

int Mesh_object::add_intermediate_texture(double t0,double t1)
{
  if(texture_type()!=0)
    set_intermediate_texture();
    
  parameters[2].get_double(0).add(t0);
  parameters[2].get_double(0).add(t1);
  return 0;
}
int Mesh_object::add_polygon_texture(const std::vector <int>& index)
{
  if(texture_type()!=0)
    set_intermediate_texture();

  parameters[2].get_int(2).add(index);
  parameters[2].get_int(2).add(-1);
  return 0;
}

int Mesh_object::get_polygon_access(int k_polygon) const
{return connectivity.get_polygon_access(k_polygon);}

int Mesh_object::get_texture_polygon_number() const
{
  if(texture_type()!=0)
    return -1;

  int polygon_count=0;
  int k=0;
  for(k=0;k<parameters[2].get_int(2).size();k++)
    if(parameters[2].get_int(2)[k]==-1)
      polygon_count++;
  return polygon_count;
}

int Mesh_object::set_texture_number(int k_texture)
{
  if(texture_type()!=0)
    set_intermediate_texture();

  parameters[2].get_int(1)[0]=k_texture;
  return 0;
}
int Mesh_object::get_texture_number() const
{
  if(texture_type()!=0)
    return -1;
  
  return parameters[2].get_int(1)[0];
}

int Mesh_object::build_manifold()
{return connectivity.build_manifold();}
int Mesh_object::build_edges()
{return connectivity.build_edges();}


Connectivity& Mesh_object::get_connectivity() {return connectivity;}
const Connectivity& Mesh_object::get_connectivity() const {return connectivity;}


int Mesh_object::increase_alpha(double d_alpha)
{
  if(color_type()==-1)
    return -1;
  
  int k=0;
  int N=parameters[1].get_double(0).size();
  double current_alpha=0.0;
  double new_alpha;
  for(k=3;k<N;k+=4)
    {
      current_alpha=parameters[1].get_double(0)[k];
      new_alpha = current_alpha+d_alpha;
      if(new_alpha>1.0)
        new_alpha=1.0;
      else if(new_alpha<0.0)
        new_alpha=0.0;
      parameters[1].get_double(0)[k]=new_alpha;
    }
  return 0;
}
int Mesh_object::decrease_alpha(double d_alpha)
{return increase_alpha(-d_alpha);}

V_3D Mesh_object::get_closest_vertex(const V_3D& center) const
{return point_set.get_closest_vertex(center);}
V_3D Mesh_object::get_closest_vertex(double x,double y,double z) const
{return point_set.get_closest_vertex(x,y,z);}
V_3D Mesh_object::get_closest_vertex(const double *x) const
{return point_set.get_closest_vertex(x);}


// int Mesh_object::check_texture() const
// {
//   if(texture_type!=0)
//     return -1;

//   int N_texture_polygon=get_texture_polygon_number();
//   int N_polygon=polygon_number();
//   if(N_texture_polygon!=N_polygon)
//     {prinntf("Failed in check_texture in Mesh_object: N_texture_polygon=%d, N_polygon=%d\n",N_texture_polygon,N_polygon);exit(-1);}

//   int N_vertex_texture=parameters[2].get_double(0).size();
//   for(int k=0;k<N_texture_polygon;k++)
//     {
//       int k_polygon_acess=get_polygon_access(k_polygon);
//       int k_index = k_polygon_acess+k_vertex;
      
//     }


// }