/home/ann2.ext/paccaulf/projsem2/ContourActif/mesh.cpp

Aller à la documentation de ce fichier.

//
// Mesh.cpp : implementation of the Mesh class 
// 

#include "mesh.h" 

// STL stuff 
#include <iostream> 
#include <fstream>
//#include <list>
#include "math.h"

// CGAL (off) 
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Triangulation_2.h>

// OpenGL stuff
#include <GL/gl.h>

// Mesh construction/destruction 
Mesh::Mesh(const char* filename) 
{
        _isValid = false;
        
        if(loadFromFile(filename))
          {
                _mesh->normalize_border();
                
                if(!_mesh->is_closed())
                  {
                        std::cout << "*** Warning: mesh is not closed ***" << std::endl;
                  }
                if(!_mesh->is_pure_triangle())
                  {
                        std::cout << "*** Warning: mesh is not pure triangle ***" << std::endl;
                  }
                if(_mesh->empty())
                  {
                        std::cout << "*** Warning: mesh is empty ***" << std::endl;
                  }
                if(!_mesh->is_valid(false,3))
                  {
                        _message = "Mesh is not combinatorially consistent";
                  }
                else
                  {
                        _isValid = true;
                        _mesh->set_index_vertices();
                        _mesh->set_index_edges();
                        _mesh->compute_normals();
                  }
          }
        else _message = "Impossible to load file";
}; 

Mesh::~Mesh()
{ 
        delete _mesh; 
}; 

void Mesh::compute_normals()
{
        _mesh->compute_normals();
};

bool Mesh::isValid() const
{
        return _isValid;
};
char* Mesh::message() const
{
        return _message;
};

std::size_t Mesh::size_of_vertices() const
{
        return _mesh->size_of_vertices();
};

std::size_t Mesh::size_of_facets() const
{
        return _mesh->size_of_facets();
};

// draw using OpenGL commands (display lists) 
void Mesh::draw(bool smooth_shading, bool use_normals, bool invert_normals) const 
{   
        // draw polygons
        glColor3f(0.0,0.5,0.8); // bleu turquoise
    
        Facet_iterator pFacet = facets_begin();
        for(;pFacet != facets_end();pFacet++)
          {
                // begin polygon assembly
                glBegin(GL_POLYGON);
                gl_draw_facet(pFacet,smooth_shading,use_normals,invert_normals);
                glEnd(); // end polygon assembly
          }
        glFlush(); 
        
};

void Mesh::gl_draw_facet(Facet_handle pFacet, bool smooth_shading, bool use_normals, bool invert_normals) const
{
        // one normal per face
        if(use_normals && !smooth_shading)
          {
                const Vector& normal = pFacet->normal();
                if (invert_normals) glNormal3f(-normal[0],-normal[1],-normal[2]);
                else glNormal3f(normal[0],normal[1],normal[2]);
          }
        
        // revolve around current face to get vertices
        Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
        do
          {
                  // one normal per vertex
                  if(use_normals && smooth_shading)
                        {
                          const Vector normal = pHalfedge->vertex()->normal();
                          if (invert_normals) glNormal3f(-normal[0],-normal[1],-normal[2]);
                          else glNormal3f(normal[0],normal[1],normal[2]);
                        }
                  
                  // polygon assembly is performed per vertex
                  const Point& point  = pHalfedge->vertex()->point();
                  glVertex3d(point[0],point[1],point[2]);
          }while(++pHalfedge != pFacet->facet_begin());
};

// superimpose edges 
void Mesh::draw_edges() const
{
        glColor3f(0.0,0.0,0.0);
        glBegin(GL_LINES);
        for(Halfedge_iterator h = _mesh->edges_begin(); h != _mesh->edges_end(); h++)
          {
                // assembly and draw line segment
                const Point& p1 = h->prev()->vertex()->point();
                const Point& p2 = h->vertex()->point();
                glVertex3f(p1[0],p1[1],p1[2]);
                glVertex3f(p2[0],p2[1],p2[2]);
          }
        glEnd();
}; 

// superimpose vertices 
void Mesh::draw_vertices() const
{
        glDisable(GL_LIGHTING);
        glPointSize(10.0);
        glColor3f(0.7,1.0,0.1); // jaune fluo
        glBegin(GL_POINTS);
        for(Vertex_iterator pPoint = _mesh->vertices_begin(); pPoint != _mesh->vertices_end(); pPoint++)
                glVertex3f(pPoint->point()[0],pPoint->point()[1],pPoint->point()[2]);
        glEnd();
        glPointSize(1.0);
        glEnable(GL_LIGHTING);
};

void Mesh::BoundingBox(double * centre, double * diameter)
// Calcule la plus petite boite englobante contenant les sommets de la surface
{
        double w, h, k;
        
        Vertex_iterator pPoint = vertices_begin();
        
        double min[3] = { pPoint->point()[0],pPoint->point()[1],pPoint->point()[2] };
        double max[3] = { pPoint->point()[0],pPoint->point()[1],pPoint->point()[2] };
        for (Vertex_iterator pPoint = vertices_begin(); pPoint != vertices_end(); pPoint++) {
                if (pPoint->point()[0] > max[0]) max[0] = pPoint->point()[0];
                if (pPoint->point()[0] < min[0]) min[0] = pPoint->point()[0];
                if (pPoint->point()[1] > max[1]) max[1] = pPoint->point()[1];
                if (pPoint->point()[1] < min[1]) min[1] = pPoint->point()[1];
                if (pPoint->point()[2] > max[2]) max[2] = pPoint->point()[2];
                if (pPoint->point()[2] < min[2]) min[2] = pPoint->point()[2];
        }
        centre[0] = (max[0]+min[0])/2.0; w = (max[0]-min[0])/2.0;
        centre[1] = (max[1]+min[1])/2.0; h = (max[1]-min[1])/2.0;
        centre[2] = (max[2]+min[2])/2.0; k = (max[2]-min[2])/2.0;
        (*diameter) = (w>h)?w:h; if (k>(*diameter)) (*diameter) = k;
        
        if (*diameter < 0.1) (*diameter) = 0.1;
};

// open document file 
bool Mesh::loadFromFile(const char* filename) 
{ 
        char* file_ext = strstr(filename,".off"); 
        if (( file_ext == NULL ) || ( strcmp(file_ext,".off") != 0 ))
          { 
                std::cerr << " *** Unexpected file format (.off required)" << std::endl;        
                return false; 
          }
        
        // allocate a new mesh 
        _mesh = new Enriched_polyhedron<Enriched_kernel,Enriched_items>; 
        CGAL_assertion(_mesh != NULL); 
        
        // read from stream  
        std::ifstream stream(filename); 
        if(!stream) 
          { 
                std::cerr << "*** Error while opening file " << filename << std::endl;  
                return false; 
          } 
        stream >> *_mesh; 
        
        if ((_mesh->size_of_vertices() == 0)|| (_mesh->size_of_facets() == 0)) {
                
                std::cout << "*** Warning: mesh is not manifold ***" << std::endl;
                return false;
        }
        
        std::cout << _mesh->size_of_vertices() << " vertices, ";
        std::cout << _mesh->size_of_facets() << " faces" << std::endl;
        
        return true; 
};

void Mesh::saveToFile(const char* filename) const
{
        std::ofstream fichier;
        
        fichier.open(filename, std::ios::out | std::ios::trunc);
        if(fichier.bad()){
                return;
        }
        
        // Write polyhedron in Object File Format (OFF).
        CGAL::set_ascii_mode( fichier);
        fichier << "OFF" << std::endl << _mesh->size_of_vertices() << ' ' 
                << _mesh->size_of_facets() << " 0" << std::endl;
        std::copy( _mesh->points_begin(), _mesh->points_end(),
                           std::ostream_iterator<Mesh::Point>( fichier, "\n"));
        for (  Mesh::Facet_iterator i = _mesh->facets_begin(); i != _mesh->facets_end(); ++i) {
                Mesh::Halfedge_around_facet_circulator j = i->facet_begin();
                // Facets in polyhedral surfaces are at least triangles.
                CGAL_assertion( CGAL::circulator_size(j) >= 3);
                fichier << CGAL::circulator_size(j) << ' ';
                do {
                        fichier << ' ' << std::distance(_mesh->vertices_begin(), j->vertex());
                } while ( ++j != i->facet_begin());
                fichier << std::endl;
        }
        fichier.close();
};

// returns the nearest Vertex_handle of the facet F to the selectedPoint, which is supposedly on a facet named F.
Mesh::Vertex_handle Mesh::nearestVertexOnFacet(Mesh::Point selectedPoint) {
        Vertex_handle SommetMini = _mesh->vertices_begin();
        for(Mesh::Vertex_iterator j = _mesh->vertices_begin(); j != _mesh->vertices_end() ; ++j) {
                if (has_smaller_distance_to_point(selectedPoint,
                                                                                  j->point(),
                                                                                  SommetMini->point())) {
                        SommetMini = j;
                }
        }
        
        std::cout << "sommet le plus proche = "
            << SommetMini->point()[0] << ' '
            << SommetMini->point()[1] << ' '
            << SommetMini->point()[2] << std::endl ;
        return SommetMini ;
};

void Mesh::drawUserSnake()
{
        
        // on affiche ce sommet
        
        Mesh::Vertex_handle vtx1 = *(_snake_vertex.begin()) ;
        //if (_mesh->_snake_vertex.)
        
        glDisable(GL_LIGHTING) ;
        for(std::vector < Mesh::Vertex_handle >::iterator i=_snake_vertex.begin() ;
                i != _snake_vertex.end() ;
                ++i) 
          {      
                glPointSize(10.0) ;
                
                glBegin(GL_POINTS) ;
                
                glColor3f(1.0,0.0,0.0) ;
                
                glVertex3f((*i)->point()[0]-0.01*dir[0],
                                   (*i)->point()[1]-0.01*dir[1],
                                   (*i)->point()[2]-0.01*dir[2]) ;
                
                glEnd() ;
                
                glPointSize(1.0) ;
                
                drawSnakePath(vtx1, *i, 0) ;
                
                vtx1 = *i ;
                
          }
        glEnable(GL_LIGHTING) ;
        //glFlush() ;
        
        
        
        
}

void Mesh::drawSnakePath(Mesh::Vertex_handle source, Mesh::Vertex_handle target, int drawPathEnds = 3)
{
        
        glPointSize(6.0f);
        
        glColor3f(1.0,0.0,0.0);
        
        glBegin(GL_POINTS);     
        
        if (drawPathEnds == 1 || drawPathEnds == 3)
                // affichage du point de depart 
                glVertex3f(source->point()[0],source->point()[1],source->point()[2]);
        if (drawPathEnds == 2 || drawPathEnds == 3)     
                // affichage du point d'arrivee
                glVertex3f(target->point()[0],target->point()[1],target->point()[2]);   
        
        glEnd();
        
        glPointSize(1.0f);
        
        glColor3f(1.0,1.0,0.0);
        
        glLineWidth(5.0);
        
        glBegin(GL_LINES);
        
        //1. on parcourt l ensemble des chemins formant le snake
        // (ex : du sommet 0 au sommet 1, puis du sommet 1 au sommet 2, etc.)
        
        for (std::vector<Mesh::Vertex_chemin>::iterator vcv = _snake_total.begin() ;
                 vcv != _snake_total.end() ;
                 vcv++) 
          {
                
                Mesh::Point p1 =  (*(vcv->begin()))->point() ;
                //2. on parcourt l ensemble des points d un chemin
                
                for (Mesh::Vertex_chemin::iterator vc = vcv->begin() ;
                         vc != vcv->end() ;
                         vc++)
                  {
                        Mesh::Point p2 =  (*vc)->point() ;
                        
                        glVertex3f(p1[0]-0.001*dir[0],
                                           p1[1]-0.001*dir[1],
                                           p1[2]-0.001*dir[2]);
                        
                        glVertex3f(p2[0]-0.001*dir[0],
                                           p2[1]-0.001*dir[0],
                                           p2[2]-0.001*dir[0]);
                        
                        p1=p2;
                  } 
          }
        
        
        glColor3f(0.3,0.6,0.4);
        
        glLineWidth(15.0);
        
        glBegin(GL_LINES);
        
        for(Mesh::Halfedge_iterator h = _mesh->border_halfedges_begin(); h != _mesh->halfedges_end(); h++)
          {
                // assembly and draw line segment
                const Mesh::Point& p1 = h->prev()->vertex()->point();
                const Mesh::Point& p2 = h->vertex()->point();
                
                glVertex3f(p1[0],p1[1],p1[2]);
                glVertex3f(p2[0],p2[1],p2[2]);
          }
        
        glEnd();
        
        glLineWidth(1.0);
}

/***********************************************************************************/
/*                                                                                 */
/*                   l e s   f o n c t i o n s   d e   v i n c e                   */
/*                                                                                 */
/***********************************************************************************/

/* drawLocalRegion
dessine la region locale en superposition au mesh entier 
c'est une copie de la procedure draw de Mesh
*/
void Mesh::drawLocalRegion() 
{
        Polyhedron* P = _snake_fini.getRegion()->mesh3D() ;
        glColor3f(0.8,0.8,0.8); // gris
        glLineWidth(2.0);
        for(Polyhedron::Facet_iterator pFacet = P->facets_begin(); pFacet != P->facets_end(); pFacet++)
          {
                glBegin(GL_POLYGON);
                
                Polyhedron::Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
                do
                  {
                          const Vector& normal = orthogonal_vector ( pHalfedge->vertex()->point(), 
                                                                                                                 pHalfedge->next()->vertex()->point(), 
                                                                                                                 pHalfedge->prev()->vertex()->point());
                          //computes an orthogonal vector of the plane defined by p, q and r, which is directed to the positive side of this plane.
                          
                          
                          //CGAL::cross_product( pHalfedge->next()->vertex()->point() - pHalfedge->vertex()->point(), pHalfedge->prev()->vertex()->point() - pHalfedge->vertex()->point() ); 
                          
                          glNormal3f(-normal[0],-normal[1],-normal[2]);
                          
                          // polygon assembly is performed per vertex
                          const Point& point  = pHalfedge->vertex()->point();
                          
                          glVertex3d(point[0]-0.01*normal[0],point[1]-0.01*normal[1],point[2]-0.01*normal[2]);
                          
                  }while(++pHalfedge != pFacet->facet_begin());
                
                glEnd();
                
          }
        glLineWidth(1.0);
        
}

/* drawLocalRegion
dessine la region locale parametrisee (2D) toute seule (sans superposition au mesh entier)
c'est une copie de la procedure draw de Mesh
*/
void Mesh::draw2dRegion() {
        Polyhedron* P = _snake_fini.getRegion()->mesh2D() ;
        glColor3f(0.8,0.8,0.8); // gris
        glLineWidth(2.0);
        for(Polyhedron::Facet_iterator pFacet = P->facets_begin(); pFacet != P->facets_end(); pFacet++)
          {
                glBegin(GL_POLYGON);
                
                Polyhedron::Halfedge_around_facet_circulator pHalfedge = pFacet->facet_begin();
                do
                  {
                          const Vector& normal = orthogonal_vector ( pHalfedge->vertex()->point(), 
                                                                                                                 pHalfedge->next()->vertex()->point(), 
                                                                                                                 pHalfedge->prev()->vertex()->point());
                          //computes an orthogonal vector of the plane defined by p, q and r, which is directed to the positive side of this plane.
                          
                          
                          //CGAL::cross_product( pHalfedge->next()->vertex()->point() - pHalfedge->vertex()->point(), pHalfedge->prev()->vertex()->point() - pHalfedge->vertex()->point() ); 
                          
                          glNormal3f(-normal[0],-normal[1],-normal[2]);
                          
                          // polygon assembly is performed per vertex
                          const Point& point  = pHalfedge->vertex()->point();
                          
                          glVertex3d(point[0],point[1],point[2]);
                          
                  }while(++pHalfedge != pFacet->facet_begin());
                
                glEnd();
                
          }
        glLineWidth(1.0);
        
}


// superimpose bounding edges
/* approx_shortest_path
calcul du plus court chemin entre source et target
resultat mis dans chemin
*/
void Mesh::approx_shortest_path(Mesh::Vertex_chemin &chemin, Vertex_handle source, Vertex_handle target) const
{
        if (source == target) 
          {
                std::cout<<"source == target !"<<std::endl ;
                exit(-1) ;
          }
        else
          {
                //std::cout<<"source != target"<<std::endl ;  
          }
        
        // maintenant on va se choisir 2 sommets du maillage et on va tracer un chemin approximativement optimal
        // de la source a la cible      
        Point prochainmin, prochainp;
        
        // nombre d iterations
        int count=0;
    
        // sommet source rajoute
        chemin.push_back(source) ;
        
        Vertex_handle pVertexCour = source ;
        
        while( (pVertexCour!=target) && (count<=2000) ){// tant que l'on n'a pas atteint la cible
        
                // on regarde autour de ce point
                Halfedge_around_vertex_circulator curr; 
                
                //vertices_circulator
                
                curr = pVertexCour->vertex_begin(); // on se place sur la premiere halfedge
        
                Vertex_handle pVermin = curr->prev()->vertex(); // le vertex courant autour
                
                prochainmin=pVermin->point();// le prochain point choisi autour pour le moment
                        
                        curr++;
                        do { 
                                prochainp = curr->prev()->vertex()->point(); 
                                
                                if( has_smaller_distance_to_point(target->point(), prochainp, prochainmin) ) {
                                        pVermin=curr->prev()->vertex(); 
                                        prochainmin=prochainp; 
                                }
                                
                        }  while (curr++ != pVertexCour->vertex_begin()); // tant que l'on n'a pas fait un tour complet
                        
                        // le sommet le plus proche de la destination est rajoute au chemin
                        chemin.push_back(pVermin) ;
                        
                        pVertexCour=pVermin; //maj
                        
                        count++;
                        
        }
        if (count >= 2000)
                std::cout << "*** Warning ***: aucun plus court chemin n'a été trouvé" << std::endl ;
        
        count=0;
        
        //  return chemin ;  
}



/*
 Structures necessaires pour determiner la local region
 Ces structures permettent de selectionner les sommets et facettes qui feront partie de la local region.
 */
// Une structure toute bete qui stocke 3 entiers (ces entiers seront des indices de facettes)
class Triangle {
public:
        int p1, p2, p3 ;
        Triangle(int p1, int p2, int p3) : p1(p1), p2(p2), p3(p3) {};
};

// operateur d affichage
std::ostream& operator<< (std::ostream& out, const Triangle& T) {
        return out << T.p1 << ' ' << T.p2 << ' ' << T.p3;
}

// classe pour associer un point a un indice
// possibilite d optimiser en etendant la classe Point_3 au lieu de creer une nouvelle classe...
class PtIndex {
public:
        Mesh::Point p ;
        int i ; 
        PtIndex() {}
        PtIndex(Mesh::Point p, int i): p(p), i(i) {}
};

// PtVector rajoute une fonction a la classe vector<PtIndex>
// On aurait pu se servir d'un set plutot qu'un vector, mais les indices des elements d'un set sont susceptibles de changer
class PtVector : public std::vector<PtIndex> {
public:
        //isPresent : renvoie -1 si le point p n'existe pas dans le vector, renvoie son indice sinon
        int isPresent(const Mesh::Point p) {
                for(std::vector<PtIndex>::iterator i = this->begin() ; i != this->end() ; i++) {
                        if (i->p == p) return i->i ;
                }
                return -1 ;
        }
};

//TrVector rajoute une fonction a la classe vector<Triangle>
class TrVector : public std::vector<Triangle> {
public:
        //isPresent : si un triangle avec les memes indices de sommets est present dans le vector, renvoie 1, sinon -1. 
        int isPresent(const Triangle Tri) {
                for(std::vector<Triangle>::iterator i = this->begin() ; i != this->end() ; i++) {
                        if (((i->p1 == Tri.p1) && (i->p2 == Tri.p2) && (i->p3 == Tri.p3)) ||
                                ((i->p1 == Tri.p2) && (i->p2 == Tri.p3) && (i->p3 == Tri.p1)) ||
                                ((i->p1 == Tri.p3) && (i->p2 == Tri.p1) && (i->p3 == Tri.p2)) ||
                                ((i->p1 == Tri.p3) && (i->p2 == Tri.p2) && (i->p3 == Tri.p1)) ||
                                ((i->p1 == Tri.p2) && (i->p2 == Tri.p1) && (i->p3 == Tri.p3)) ||
                                ((i->p1 == Tri.p1) && (i->p2 == Tri.p3) && (i->p3 == Tri.p2)))
                                return 1 ;
                }
                return -1 ;
        }
};

/*
 Local_region
 calcul de la region locale associee au snake Snake
 Selectionne les sommets et facettes du mesh d origine, les stocke en memoire, puis ecrit un fichier au format off.
 Ensuite on convertit simplement ce fichier en polyhedron grace a l operateur >>.
 Il y a evidemment une perte de performance, et un risque sur la stabilite du programme 
 (par ex si le dossier n est pas libre en ecriture, ou si le disque est plein).
 /!\ les facettes autres que triangles ne sont pas gerees !!!
 */
Mesh::Polyhedron Mesh::Local_region( std::vector< Vertex_chemin >& Snake ) // only first release
{
        int vecsize = 0 ;
        // parcours de l'ensemble des sommets de la mesh jusqu'a trouver le premier sommet du snake
        
        PtVector S;
        TrVector T ;
        PtIndex pind1, pind2, pind3 ;
        int indice, ind1, ind2, ind3 ;
        
        for(std::vector< Vertex_chemin >::const_iterator chem_it=Snake.begin(); chem_it!=Snake.end(); chem_it++)
          {
                for(std::vector< Vertex_handle >::const_iterator it=chem_it->begin(); it!=chem_it->end(); it++)
                  {
                        //std::cout <<  P.size_of_vertices()  <<std::endl;

                        //vertices_circulator
                        Halfedge_around_vertex_circulator curr = (*it)->vertex_begin(); // on se place sur la premiere halfedge
                        
                        Point p1(curr->vertex()->point()); // point du snake
                        indice = S.isPresent(p1) ;
                        if (indice == -1) 
                          {
                                ind1 = vecsize ;
                                pind1 = PtIndex(p1, vecsize++) ;
                                S.push_back(pind1) ;
                                //std::cout<<"SOMMET AJOUTE, ind = "<< vecsize-1 <<std::endl ;
                          } 
                        else
                          {
                                ind1 = indice ;
                                //std::cout<<"SOMMET PRESENT, ind = "<< indice <<std::endl ;
                          }
                        
                        do
                          {     
                                  
                                  Point p2(curr->prev()->vertex()->point()); // un sommet incident
                                  indice = S.isPresent(p2) ;
                                  if (indice == -1) 
                                        {
                                          ind2 = vecsize ;
                                          pind2 = PtIndex(p2, vecsize++) ;
                                          S.push_back(pind2) ;
                                          //std::cout<<"SOMMET AJOUTE, ind = "<< vecsize-1 <<std::endl ;
                                        } 
                                  else
                                        {
                                          ind2 = indice ;
                                          //std::cout<<"SOMMET PRESENT, ind = "<< indice <<std::endl ;
                                        }
                                  
                                  
                                  curr++;
                                  
                                  Point p3(curr->prev()->vertex()->point()); // un autre sommet incident
                                  indice = S.isPresent(p3) ;
                                  if (indice == -1) 
                                        {
                                          ind3 = vecsize ;
                                          pind3 = PtIndex(p3, vecsize++) ;
                                          S.push_back(pind3) ;
                                          //std::cout<<"SOMMET AJOUTE, ind = "<< vecsize-1 <<std::endl ;
                                        } 
                                  else
                                        {
                                          ind3 = indice ;
                                          //std::cout<<"SOMMET PRESENT, ind = "<< indice <<std::endl ;
                                        }
                                  //si a ce niveau on a encore faceMayBePresent == true, alors la face est deja presente.
                                  
                                  Triangle facette = Triangle(ind1, ind2, ind3) ;
                                  if (T.isPresent(facette) != -1) 
                                        {
                                          /*
                                           std::cout<<"facette deja presente !!"<<std::endl ;
                                           std::cout<<"barycentre == "
                                           <<(p1.x()+p2.x()+p3.x())/3<<" "
                                           <<(p1.y()+p2.y()+p3.y())/3<<" "
                                           <<(p1.z()+p2.z()+p3.z())/3<<std::endl ;
                                           /**/
                                        } 
                                           else
                                           // le triangle n'est pas present, on peut l'ajouter.
                                           T.push_back(Triangle(ind1, ind2, ind3)) ;
                                           
                                           
                          }  while (curr != (*it)->vertex_begin() ); // tant que l'on n'a pas fait un tour complet      
                  }
          }
                                           // On a rempli les structures de stockage des points et des facettes. Maintenant on peut envoyer tout ca dans un fichier
                                           
                                           const char* filename = "SNAKETEMP.off" ;
                                           std::ofstream fichier;
                                           
                                           fichier.open(filename, std::ios::out | std::ios::trunc);
                                           if(fichier.bad()){
                                                   std::cout<<"erreur lecture fichier"<<std::endl ;
                                                   exit(-1) ;
                                           }
                                           
                                           // Write polyhedron in Object File Format (OFF).
                                           CGAL::set_ascii_mode( fichier);
                                           fichier << "OFF" << std::endl << S.size() << ' ' 
                                           << T.size() << " 0" << std::endl;
                                           for(PtVector::iterator ptv = S.begin() ; ptv != S.end() ; ptv++) {
                                                   fichier<<ptv->p<<std::endl ;
                                           }
                                           for (  std::vector<Triangle>::iterator i = T.begin(); i != T.end(); i++) {
                                                   // Alors la si jamais la facette n'est pas un triangle, on est mal. 
                                                   // Mais c'est plus tot qu'il fallait le verifier.
                                                   fichier << '3' << ' ';
                                                   
                                                   //std::cout << *i << std::endl ;
                                                   fichier << *i;
                                                   fichier << std::endl;
                                           }
                                           fichier.close();
                                           
                                           // lecture du fichier (utilisation de l operateur >> vers Polyhedron)
                                           Mesh M(filename) ;
                                           return *M._mesh ;
}

// fonctionne pour des Polyhedron (pas seulement pour des Enriched_polyhedron)
/*
 void Mesh::Trace_Local_region( Mesh::Polyhedron P ) // OpenGL 
{
        std::cout << "obsolete !!" << std::endl ;
}
/**/

---