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

Aller à la documentation de ce fichier.

#include "region.h"
#include "resollu.h"

//tempo pour le debug
#include <fstream>
#include <iostream> 

//on passe en parametre un vertex_handle sur un point 2D
double Region::Compute_nrj_sommet_sur_mesh (Vertex_handle pV) 
{
        double min;
        
        Vertex_handle pV3D; // le vertex handle 3D correspondant
        pV3D=pV;
        this->vh3Dfrom2D(pV3D); // pV3D est mis a jour
        
        // the normal at a vertex
        Vector nv = pV3D->normal(); // Normal est un Vector_3
        
        Halfedge_around_vertex_circulator h = pV3D->vertex_begin();
        
        min = nv * h->facet()->normal();
        
        do {
                h++;
                
                if( (nv * h->facet()->normal()) < min ) min = nv * h->facet()->normal();
                
        } while(h != pV3D->vertex_begin() );
        
        return min;
}

//on passe en parametre un vertex_handle sur un point 2D
double Region::Compute_nrj_sommet (Vertex_handle pV) 
{
        double nrj=0; // notre nrj
        int count=0;
        
        Vertex_handle pV3D; // le vertex handle 3D correspondant
        pV3D=pV;
        this->vh3Dfrom2D(pV3D); // pV3D est mis a jour
        
        Vertex_handle SommetMini = _mesh3D.vertices_begin(); // le premier sommet le plus proche de pV
        
        for(Polyhedron::Vertex_iterator j = _mesh3D.vertices_begin(); j != _mesh3D.vertices_end() ; ++j) {
                
                if (has_smaller_distance_to_point(pV3D->point(), j->point(), SommetMini->point())) SommetMini = j;
        }
        
        // ici SommetMini est le vertex le plus proche de pV
        
        Halfedge_around_vertex_circulator h = SommetMini->vertex_begin(); // on prend la halfedge incidente
        
        Facet_handle pF = h->facet(); // on prend la facet incidente a la halfedge
        
        double M[3][MAT_SIZE]; // car on ne travaille que sur des triangles
        double U[3];
        double B[3]; // pour notre sol
        
        Halfedge_around_facet_circulator h2 = pF->facet_begin(); 
        
        if (CGAL::circulator_size(h2)!=3) { std::cout << "ERROR in Compute_nrj_sommet: facet must be triangle" << std::endl;  exit(-1);  }
        
        do {
                M[count][0]=h2->vertex()->point()[0];
                M[count][1]=h2->vertex()->point()[1];
                M[count][2]=h2->vertex()->point()[2];
                
                U[count++] = Compute_nrj_sommet_sur_mesh( h2->vertex() );
                
        } while(h2 != pF->facet_begin() );
        
        if(! ResoudreSysteme(M, U, B, 3) )  { std::cout << "ERROR in Compute_nrj_sommet: problem occured during matrix resolution" << std::endl;  exit(-1);  }
        
        nrj= B[0]*pV->point()[0]+B[1]*pV->point()[1]+B[2]*pV->point()[2];
        
        return nrj;
}














class Border_Edge_circulator {
public:
        typedef Enriched_polyhedron<Enriched_kernel,Enriched_items> Polyhedron ;
        typedef Polyhedron::Vertex_handle Vertex_handle;
        typedef Polyhedron::Halfedge_handle Halfedge_handle;
        typedef Polyhedron::Halfedge_around_vertex_circulator Halfedge_around_vertex_circulator;
        typedef Polyhedron::Edge_iterator Edge_iterator;
        
private:
                Halfedge_handle prec ;
public:
        Border_Edge_circulator(Edge_iterator Ei) : prec(Ei) {
        }  
        
        void operator++() {
                Vertex_handle vh = prec->vertex() ;
                Halfedge_around_vertex_circulator h = vh->vertex_begin() ;
                bool DEBUGprecNonModifie = true ;
                do
                  {
                          if ((h->is_border()) && (h->opposite()->vertex() != prec->opposite()->vertex()))
                                {
                                  prec=h->opposite() ;
                                  DEBUGprecNonModifie = false ;
                                }
                          h++ ;
                  } while (h != vh->vertex_begin()) ;
                if (DEBUGprecNonModifie)
                  {
                        std::cout << "***dans Border_Edge_iterator : prec non modifie" << std::endl ;
                        exit(-1) ;
                  }
        }
        
        Halfedge_handle operator*() {
                return prec ;
        }

        bool operator==(Halfedge_handle h) {
                return (prec == h) ;
        }

        bool operator!=(Halfedge_handle h) {
                return (prec != h) ;
        }       
};







// Determination de la parametrisation 2D a partir de la mesh 3D 
void Region::Compute_Parameterization()
{
        typedef CGAL::Simple_cartesian<double> Kernel;
        typedef CGAL::Filtered_kernel_adaptor<Kernel> K; 
        
        // on commence par verifier que la region locale est bien ouverte
        assert ( this->valide3D() );
        
        K::Compute_squared_distance_3 squared_distance; 
        
        // Parametrisation sur un carre
        
        int nbv = _mesh3D.size_of_vertices();
        int bs = _mesh3D.size_of_border_edges(); // bs = nb border edges = nb border vertices
        int bhe = _mesh3D.size_of_border_halfedges();
        int nbf = _mesh3D.size_of_facets();  
        //      int iv = nbv - bs; // iv = nb de vertices internes
        
        std::cout << "nombre de sommets de la region locale " << nbv << std::endl;
        std::cout << "nombre d'aretes au bord " << bs << std::endl;
        std::cout << "nombre de demi-aretes au bord " << bhe << std::endl;
        std::cout << "nombre de facets " << nbf << std::endl;
        
        Point c1,c2,c3,c4; // les 4 coins
                                
        double d12=0,d23=0,d34=0,d41=0; // les longueurs des 4 cotes de la local region init zero
        double dist_totale=0;
        double cmax=0;
        double dist=0;
        
    //  double inter;
    //  double distiter;
        
        
        // pour separer la bordure de l'interieur (important pour le parcours des aretes)
        _mesh2D.normalize_border(); 
        
        
        // approx d'un carre de param interessant (a justifier dans le compte rendu) 
        if( nbv<801 ) cmax=6;
        else cmax=40;
        
        // valeurs des 4 points du carre
        c1=Point(0,cmax,0);
        c2=Point(0,0,0);
        c3=Point(cmax,0,0);
        c4=Point(cmax,cmax,0);
        
        
        
        Point p1,p2;
        
        //      glColor3f(0.8,0.2,0.2); // red
        //   glLineWidth(20.0);
        
        //glBegin(GL_LINES);    
        // glPointSize(18.0f);
    // glColor3f(1.0,0.2,0.4); 
        
        // calcul du perimetre de la region local       
        int count=0;
        Border_Edge_circulator h(_mesh2D.border_edges_begin());
        //Edge_iterator h = _mesh2D.border_edges_begin();
        
        p2= (*h)->vertex()->point(); // pour avoir la premiere dist a zero      
        do
          {
                  p1=p2;
                  p2= (*h)->vertex()->point();
                  
                  //std::cout << "point [" << p2[0] << "][" << p2[1] << "][" << p2[2] << "]" << std::endl;
                  
                  dist_totale += sqrt(squared_distance( p1, p2));
                
                  ++h ;
                  
          } while (h != _mesh2D.border_edges_begin()) ;
        
/*      
                for(; h != _mesh2D.edges_end(); h++)
 {      
                        p1=p2;
                        p2= (*h)->vertex()->point();
                        
                        //std::cout << "point [" << p2[0] << "][" << p2[1] << "][" << p2[2] << "]" << std::endl;
                        
                        dist_totale += sqrt(squared_distance( p1, p2));
                        
 }
        */
        
        
        h = _mesh2D.border_edges_begin();
        p2= (*h)->vertex()->point(); // pour avoir la premiere dist a zero      
        
        do
          {
                  p1=p2;
                  p2= (*h)->vertex()->point();
                  
                  dist = sqrt(squared_distance( p1, p2));
                  
                  // computation of edge distances
                  if(d12 + dist <= dist_totale/4) d12 += dist;
                  else if (d23 + dist <= dist_totale/4) d23 += dist;
                  else if (d34 + dist <= dist_totale/4) d34 += dist;
                  else d41 += dist;
                  
                  ++h;
                  
          } while (h != _mesh2D.border_edges_begin()) ;
        
/*      for(; h != _mesh2D.edges_end(); h++)
          {     
                p1=p2;
                p2= (*h)->vertex()->point();
                
                dist = sqrt(squared_distance( p1, p2));
                
                // computation of edge distances
                if(d12 + dist <= dist_totale/4) d12 += dist;
                else if (d23 + dist <= dist_totale/4) d23 += dist;
                else if (d34 + dist <= dist_totale/4) d34 += dist;
                else d41 += dist;
          }
*/      
        
        //assert(dist_totale == d12+d23+d34+d41);
        
        // pour Žviter les pbs de div par zero
        if(d12 <= 0) d12=0.00001;
        if(d23 <= 0) d23=0.00001;
        if(d34 <= 0) d34=0.00001;
        if(d41 <= 0) d41=0.00001;
        
        
        /*
         std::cout << "d12 " << d12 << std::endl;
         std::cout << "d23 " << d23 << std::endl;
         std::cout << "d34 " << d34 << std::endl;
         std::cout << "d41 " << d41 << std::endl;
         std::cout << "d12+d23+d34+d41 " << d12+d23+d34+d41 << std::endl;
         std::cout << "DISTANCE TOTALE " << dist_totale << std::endl;
         
         */     
        
        // maintenant on fixe les points du bord
        h = _mesh2D.border_edges_begin();
        
        count=0;
        dist=0;
        p2=(*h)->vertex()->point(); // pour avoir la premiere dist a zero
        while (1) {
                p1=p2;
                p2 = (*h)->vertex()->point();
                
                dist += sqrt(squared_distance( p1, p2 ));
                
                //std::cout << "debuggage :dist/12" << dist/d12 << "distinter :" <<  dist << std::endl;
                
                if (dist >= d12) break;
                
                (*h)->vertex()->point()=c1+(c2-c1)*dist/d12;
                std::cout << "Nouveau point [" << (*h)->vertex()->point()[0] << "][" << (*h)->vertex()->point()[1] << "][" << (*h)->vertex()->point()[2] << "]" << std::endl;
                
                ++h;
                
                count++;
        }
                                
        dist=0;
        p2=(*h)->vertex()->point(); // pour avoir la premiere dist a zero
        while (1) {
                p1=p2;
                p2 = (*h)->vertex()->point();
                dist += sqrt(squared_distance( p1, p2 ));
                
                //std::cout << "debuggage :dist/d23" << dist/d23 << "distinter :" <<  dist << std::endl;
                                        
                if (dist >= d23) break;

                (*h)->vertex()->point()=c2+(c3-c2)*dist/d23;
                std::cout << "Nouveau point [" << (*h)->vertex()->point()[0] << "][" << (*h)->vertex()->point()[1] << "][" << (*h)->vertex()->point()[2] << "]" << std::endl;                                   
                
                ++h;
                
                count++;
        }
                                
        dist=0;
        p2=(*h)->vertex()->point(); // pour avoir la premiere dist a zero
        while (1) {
                p1=p2;
                p2 = (*h)->vertex()->point();
                dist += sqrt(squared_distance( p1, p2 ));
                
                //std::cout << "debuggage :dist/d34" << dist/d34 << "distinter :" <<  dist << std::endl;
                
                if (dist >= d34) break;
                
                (*h)->vertex()->point()=c3+(c4-c3)*dist/d34;
                std::cout << "Nouveau point [" << (*h)->vertex()->point()[0] << "][" << (*h)->vertex()->point()[1] << "][" << (*h)->vertex()->point()[2] << "]" << std::endl;
                
                
                ++h;
                count++;
        }
        
        dist=0;
        p2=(*h)->vertex()->point(); // pour avoir la premiere dist a zero
        
        do {
                p1=p2;
                p2 = (*h)->vertex()->point();
                dist += sqrt(squared_distance( p1, p2 )) ;
                
                (*h)->vertex()->point()=c4+(c1-c4)*dist/d41;
                std::cout << "Nouveau point [" << (*h)->vertex()->point()[0] << "][" << (*h)->vertex()->point()[1] << "][" << (*h)->vertex()->point()[2] << "]" << std::endl;
                
                
                ++h;
                
                count++;
        } while (h != _mesh2D.border_edges_begin()) ; 
        
        std::cout << "nb de points modifies :" << count << std::endl;
        
        
        
        
        
        // maj
        count=0;
        dist=0; 
        
        Trace2file("AAAavantParamInt.off") ;
        
        //i sommet interieur
        //j sommet quelconque
        //calcul des lambda(i,j) de maniere uniforme :
        // qqsoit i,j, lambda(i,j)=adjacence(sommet i)
        //bool isVertexBorder ;
        int Dimension = _mesh2D.size_of_vertices() - _mesh2D.size_of_border_edges() ; // nb de sommets interieurs
        double Kmatx[Dimension][MAT_SIZE], Kmaty[Dimension][MAT_SIZE] ;
        double Ux[Dimension], Uy[Dimension] ;
        double Bx[Dimension], By[Dimension] ;
        int icount = 0 ;
        
        for(Vertex_iterator i = _mesh2D.vertices_begin(); i != _mesh2D.vertices_end(); i++)
          {
                std::cout << "un sommet ! " ;
                //an incident halfedge that starts from v.
                //i->halfedge().opposite() ;
                
                //determiner si i est sur la bordure ou non
                /*
                isVertexBorder=false ;
                Halfedge_around_vertex_circulator h = i->vertex_begin() ;
                do
                  {
                          if ((*h)->opposite()->is_border())
                                {
                                  isVertexBorder=true ;
                                  std::cout << "b" ; 
                                } 
                          else
                                  std::cout << "." ;
                          h++ ;
                  } while (h != i->vertex_begin()) ;            
                std::cout << i->vertex_degree() << " " << i->halfedge()->opposite()->is_border() << std::endl ;
                if (isVertexBorder)
                 */
                if (isBorderVertex(i))
                  {
                        
                  } 
                else
                  {
                        //i n est pas sur la bordure
                        double sommeUx = 0 ;
                        double sommeUy = 0 ;
                        //parcourir sommets adjacents a i (nommes vh)
                        Halfedge_around_vertex_circulator h = i->vertex_begin() ;
                        do
                          {
                                  Vertex_handle vh = h->opposite()->vertex() ;
                                  //determiner si vh est sur la bordure ou non
                                  std::cout << "   un sommet ! " ;
                                  /*
                                  bool isVHBorderVertex=false ;
                                  Halfedge_around_vertex_circulator hh = vh->vertex_begin() ;
                                  do
                                        {
                                                if (hh->opposite()->is_border())
                                                  {
                                                        isVHBorderVertex=true ;
                                                        std::cout << "b" ; 
                                                  } 
                                                else
                                                        std::cout << "." ;
                                                hh++ ;
                                        } while (hh != vh->vertex_begin()) ;
                                  if (isVHBorderVertex)
                                   */
                                  if (isBorderVertex(vh))
                                        {
                                          //si vh est sur la bordure on ajoute le point vh a la somme sommeU, avec un coeff
                                          //egal a 1/adjacence(i) cad a 1/lambda(i,vh)
                                          sommeUx += 1/((double) (i->vertex_degree())) * vh->point()[0] ;
                                          sommeUy += 1/((double) (i->vertex_degree())) * vh->point()[1] ;
                                          std::cout << sommeUx << " " << sommeUy << " " << vh->halfedge()->opposite()->is_border() ;
                                        }
                                  
                                  std::cout << std::endl ;
                                  h++;
                          } while (h != i->vertex_begin()) ;
                        Ux[icount] = sommeUx ;
                        Uy[icount] = sommeUy ;
                        
                        //remplissage de Kmat
                        //tres facile si on choisit la parametrisation uniforme
                        //cad lambda(i,j) = 1/adjacence(i)
                        //dans ce cas, la ligne i de Kmat vaut 1/adjacence(i) sauf a la colonne i ou le coeff de Kmat vaut 1.
                        std::cout<<"["<<icount<<"]   " ;
                        for (int iKmat = 0 ; iKmat < MAT_SIZE ; iKmat++)
                          {
                                if (iKmat != icount)
                                  {
                                        if (isEdge(icount,iKmat))
                                          {
                                                std::cout<<i->vertex_degree() ;
                                                Kmatx[icount][iKmat] = -1/((double) i->vertex_degree()) ;
                                                Kmaty[icount][iKmat] = -1/((double) i->vertex_degree()) ;
                                          }
                                        else
                                          {
                                                std::cout<<"0" ;
                                                Kmatx[icount][iKmat] = 0 ;
                                                Kmaty[icount][iKmat] = 0 ;
                                          }
                                  }
                                else
                                  {
                                        std::cout<<"1" ;
                                        Kmatx[icount][iKmat] = 1 ;
                                        Kmaty[icount][iKmat] = 1 ;
                                  }
                          }                                             
                        std::cout << std::endl ;
                        icount++ ;
                  }
          }
        std::cout << std::endl ;
        std::cout << "icount="<<icount<<", Dimension="<<Dimension<<std::endl ;
        
        // RESOLUTION SYSTEME !!!
        if (!ResoudreSysteme(Kmatx, Ux, Bx, Dimension))
                std::cout << "***ERREUR dans ResoudreSysteme en x" ;
        if (!ResoudreSysteme(Kmaty, Uy, By, Dimension))
                std::cout << "***ERREUR dans ResoudreSysteme en y" ;

        //MISE A JOUR POINTS 2D !!!
        int iB = 0 ;
        for(Vertex_iterator i2d = _mesh2D.vertices_begin(); i2d != _mesh2D.vertices_end(); i2d++)
          {
                //an incident halfedge that starts from v.
                //i->halfedge().opposite() ;
                
                //determiner si i est sur la bordure ou non
                /*
                isVertexBorder=false ;
                Halfedge_around_vertex_circulator h2d = i2d->vertex_begin() ;
                do
                  {
                          if (h2d->opposite()->is_border())
                                {
                                  isVertexBorder=true ;
                                } 
                          h2d++ ;
                  } while (h2d != i2d->vertex_begin()) ;                
                if (!isVertexBorder)
                 */
                if (!isBorderVertex(i2d))
                  {
                        Point p2d=Point(Bx[iB], By[iB], 0) ;
                        i2d->point()=p2d ;
                        iB++ ;
                  }
          }     
        
        Trace2file("AAAapresParamInt.off") ;
        
                                
                                // Maintenant il faut eliminer les concavites (si on utilise l'edge tweaking method)
                                
        
                                
}

// normalement, devrait etre membre de la classe Enriched_polyhedron
bool Region::isBorderVertex(Vertex_handle vh)
{
        Halfedge_around_vertex_circulator hh = vh->vertex_begin() ;
        do
          {
                  if ((hh->is_border()) || (hh->opposite()->is_border()))
                        {
                          return true ;
                        } 
                  hh++ ;
          } while (hh != vh->vertex_begin()) ;
        return false ;
}


// normalement, devrait etre membre de la classe Enriched_polyhedron
bool Region::isEdge(int i, int j)
{
        if (i==j)
                return false ;
        int count=0 ;
        Vertex_handle vi, vj ;
        for(Vertex_iterator v = _mesh2D.vertices_begin() ; v != _mesh2D.vertices_end() ; v++)
          {
                if (count==i)
                        vi=v ;
                if (count==j)
                        vj=v ;
                if (!isBorderVertex(v))
                  {
                        count++ ;
                  }
          }
        Halfedge_around_vertex_circulator hh = vi->vertex_begin() ;
        do
          {
                  if (hh->opposite()->vertex() == vj)
                          return true ;
                  hh++ ;
          } while (hh != vi->vertex_begin()) ;
        return false ;
}







void Region::Trace2file(char* filename) 
{
        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 << _mesh2D.size_of_vertices() << ' ' 
                << _mesh2D.size_of_facets() << " 0" << std::endl;
        std::copy( _mesh2D.points_begin(), _mesh2D.points_end(),
                           std::ostream_iterator<Point>( fichier, "\n"));
        for (  Facet_iterator i = _mesh2D.facets_begin(); i != _mesh2D.facets_end(); ++i) {
                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(_mesh2D.vertices_begin(), j->vertex());
                } while ( ++j != i->facet_begin());
                fichier << std::endl;
        }
        fichier.close();
        
}



// Methode qui renvoie un vertexhandle sur un point 2D a partir d'un vertex handle sur 
// un point 3D
void Region::vh2Dfrom3D (Vertex_handle& vh)
{
        Polyhedron::Vertex_iterator j2D = _mesh2D.vertices_begin();
        
        for(Polyhedron::Vertex_iterator j = _mesh3D.vertices_begin(); j != _mesh3D.vertices_end() ; j++) 
        {       
                /*
                std::cout << "point 3D[" << j->point()[0] << "][" << j->point()[1] << "][" << j->point()[2] << "]" << std::endl;
                std::cout << "point 2D[" << j2D->point()[0] << "][" << j2D->point()[1] << "][" << j2D->point()[2] << "]" << std::endl;
                 */
                if ( j->point() == vh->point() ){ vh=j2D; break; }
                
                ++j2D;
        }
}


// Methode qui renvoie un vertexhandle sur un point 3D a partir d'un vertex handle sur 
// un point 2D
void Region::vh3Dfrom2D(Vertex_handle& vh)
{
        Polyhedron::Vertex_iterator j3D = _mesh3D.vertices_begin();
        
        for(Polyhedron::Vertex_iterator j = _mesh2D.vertices_begin(); j != _mesh2D.vertices_end() ; j++) 
        {               
                if ( j->point() == vh->point()  ){ vh=j3D; break; }
                
                ++j3D;
        }
}



void Region::Point3D_from_2D(Point& p)
{
        int count=0;
        
        // on regarde d'abord si notre point p est sur la mesh 2D
        Polyhedron::Vertex_iterator j3D = _mesh3D.vertices_begin();
        for(Vertex_iterator pVertex = _mesh2D.vertices_begin(); pVertex != _mesh2D.vertices_end(); pVertex++)
        {
                // si notre point est un sommet de la mesh 2D
                if (p == pVertex->point() ) { p=j3D->point(); return; }
                j3D++;
        }
        
        Facet_handle pF2D;
        Facet_handle pF3D=_mesh3D.facets_begin();
        // Ici on est dans le cas ou le point p est sur une facet 
        for (  Facet_iterator i = _mesh2D.facets_begin(); i != _mesh2D.facets_end(); i++)
        {
                if( this->Appartient_facet(i,p) ) { pF2D=i;  break; }
                
                pF3D++;
        }       
                
/*      
        Vertex_handle SommetMini = _mesh2D.vertices_begin(); // le premier sommet 2Dle plus proche de p
        
        j3D = _mesh3D.vertices_begin();
        
        Vertex_handle SommetMini3D = _mesh3D.vertices_begin(); // le premier sommet 3D le plus proche de p (par correspondance)

        // determiner le sommet de _mesh3D qui correspond au plus proche de p.
        Polyhedron::Vertex_iterator j = _mesh2D.vertices_begin();
        while( j != _mesh2D.vertices_end()) 
        {               
                if (has_smaller_distance_to_point(p, j->point(), SommetMini->point()))
                { SommetMini = j; SommetMini3D = j3D;  }
                j++;
                j3D++;
        }
        
        std::cout << "SommetMini3D determine" << std::endl;
        
        // ici SommetMini3D est le vertex le plus proche du point p
        
        Halfedge_around_vertex_circulator h = SommetMini->vertex_begin(); // on prend une halfedge incidente
        

        Facet_handle pF = h->facet(); // on prend la facet incidente a la halfedge
        
  */
                
                
        double (*M)[200] = new double [2][200]; 
        
        double U[2];
        double B[2]; // pour notre sol  
        Point Mespt2D[3]; 
        
        Point Mespt3D[3];
        
        Halfedge_around_facet_circulator h2 = pF2D->facet_begin(); 
        
        if (CGAL::circulator_size(h2)!=3) { std::cout << "ERROR in Compute_nrj_sommet: facet must be triangle" << std::endl;  exit(-1);  }
        
        do {
                Mespt2D[count++]=h2->vertex()->point();
                h2++;
        } while(h2 != pF2D->facet_begin() );
        
        // on remplit notre matrice
        M[0][0]=Mespt2D[0][0]-Mespt2D[2][0];
        M[0][1]=Mespt2D[1][0]-Mespt2D[2][0];
        
        M[1][0]=Mespt2D[0][1]-Mespt2D[2][1];
        M[1][1]=Mespt2D[1][1]-Mespt2D[2][1];
        
        U[0]=p[0]-Mespt2D[2][0];
        U[1]=p[1]-Mespt2D[2][1];
        
        std::cout << "Avant appel a ResoudreSysteme" << std::endl;
        if(! ResoudreSysteme(M, U, B, 2) )  { std::cout << "ERROR in Point3D_from_2D: problem occured during matrix resolution" << std::endl;  exit(-1);  }     
        std::cout << "Apres appel a ResoudreSysteme" << std::endl;
        
        // maintenant on va se determiner les 3 points 3D correspondants a nos 3 points 2D
        
        //h = SommetMini3D->vertex_begin(); // on prend la halfedge incidente
        
        //pF = h->facet(); // on prend la facet incidente a la halfedge
        
        
        h2 = pF3D->facet_begin(); 
        
        count=0;
        
        if (CGAL::circulator_size(h2)!=3) { std::cout << "ERROR in Compute_nrj_sommet: facet must be triangle" << std::endl;  exit(-1);  }
        std::cout << "Avant le 2eme do" << std::endl;
        do {
                Mespt3D[count++]=h2->vertex()->point();
                h2++;
        } while(h2 != pF3D->facet_begin() );
        
        std::cout << "creation de x, y et z" << std::endl;
        double x = B[0]*Mespt3D[0][0]+B[1]*Mespt3D[1][0]+(1-B[0]+B[1])*Mespt3D[2][0];
        double y = B[0]*Mespt3D[0][1]+B[1]*Mespt3D[1][1]+(1-B[0]+B[1])*Mespt3D[2][1];   
        double z = B[0]*Mespt3D[0][2]+B[1]*Mespt3D[1][2]+(1-B[0]+B[1])*Mespt3D[2][2];   
        std::cout << "creation de x, y et z terminee" << std::endl;
        
        p = Point(x,y,z);
        
        delete [] M;// liberation
        
}
        
                        

---