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

Aller à la documentation de ce fichier.

// STL stuff 
#include <iostream> 
#include <fstream>

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

#include "resollu.h"


#include "snake.h"



Snake::Snake() //Default constructor (ca ne sert a rien de reserver de la memoire pour rien...)
{
        Pas_Gamma=1;
        
        Snake_VH = std::vector< Vertex_handle >();
        dim=0;
        
        WeightsA = NULL;
        WeightsB = NULL;
        
        A = NULL; 
        
    R_is_initialized=false; // la region associee au snake n'est pas valide pour le moment
}
        
    
/*
  The user input is utilized for initially approximating the features 
  detection process
*/
Snake::Snake( std::vector< Vertex_handle > S , Polyhedron localRegion = Polyhedron())
{
        
        if (!localRegion.empty())
        {
                this->R=Region(localRegion) ; // contruction de la region
                R.Compute_Parameterization() ; // parametrisation
                R_is_initialized=true;  
                
                Snake_VH = std::vector< Vertex_handle >(S); // utilisation du constructeur par copie
                
                // Passage aux points 2D
                for(std::vector< Vertex_handle >::iterator it=Snake_VH.begin(); it != Snake_VH.end(); it++)
                             R.vh2Dfrom3D(*it); 
        
                dim=Snake_VH.size();
                
                WeightsA = new double[dim];
                WeightsB = new double[dim];
                
                
                for(int i=0; i<dim; i++)
                {
                        WeightsA[i]=1;
                        WeightsB[i]=1;                     
                }
                
                
                A = new double [dim][200]; // allocation memoire
                this->Maj_A();
        }
        else
        {
                R_is_initialized=false; // la region associee au snake n'est pas valide pour le moment
                
                Snake_VH = std::vector< Vertex_handle >();
                dim=0;
                
                WeightsA = NULL;
                WeightsB = NULL;
                
                A = NULL; 
        }
        
}

/*
  The user input is utilized for initially approximating the features 
  detection process
*/

Snake::Snake( std::vector< std::vector< Vertex_handle > > S, Polyhedron localRegion = Polyhedron() )
{
        //defensif
        
        Snake_VH.clear() ;
        
        Snake_VH = std::vector< Vertex_handle >();
        
        std::vector< Vertex_handle > Sninter = std::vector< Vertex_handle >();
        
        // besoin de sortir unSommet des deux for imbriques pour pouvoir faire le push_back a la suite
        std::vector< Vertex_handle >::iterator unSommet ;
        for (std::vector< std::vector< Vertex_handle > >::iterator unSegment = S.begin() ; 
                 unSegment != S.end() ; 
                 unSegment++) 
          {
                unSommet = unSegment->begin() ;
                for ( ; unSommet != unSegment->end() ; unSommet++) 
                  {
                        if (*unSommet != *(unSegment->end())) Sninter.push_back(*unSommet) ;
                  }
          }
        
        
        if (!localRegion.empty())
        {
                this->R=Region(localRegion) ; // contruction de la region
                R.Compute_Parameterization() ; // parametrisation
                R_is_initialized=true;  
                
                Vertex_handle vh;
                // Passage aux points 2D
                for(std::vector< Vertex_handle >::iterator it=Sninter.begin(); it != Sninter.end(); it++)
                {
                        vh=*it;
                        R.vh2Dfrom3D(vh); 
                        Snake_VH.push_back(vh);
                }
                
                dim=Snake_VH.size();
                
                WeightsA = new double[dim];
                WeightsB = new double[dim];
                
                
                for(int i=0; i<dim; i++)
                {
                        WeightsA[i]=1;
                        WeightsB[i]=1;                     
                }
                
                
                A = new double [dim][200]; // allocation memoire
                this->Maj_A();
                
                this->Energy_minimization(20); 

        }
        else
        {
                R_is_initialized=false; // la region associee au snake n'est pas valide pour le moment
                
                Snake_VH.clear() ;
                Snake_VH = std::vector< Vertex_handle >();
                
                dim=0;
                
                WeightsA = NULL;
                WeightsB = NULL;
                
                A = NULL; 
        }
        
        
}

Snake::~Snake() //Destructor
{
  delete [] WeightsA;  
  delete [] WeightsB; 
  delete [] A;  
}
        
// Methods
void Snake::Affiche_snake() const
{   
        if(Snake_VH.size()==0) std::cout << "Snake vide" << std::endl;
        else {
                int count=1;
                for(std::vector< Vertex_handle >::const_iterator it=Snake_VH.begin(); it != Snake_VH.end(); it++)
                {
                        std::cout << "Point numero " << count++ << std::endl;
                        std::cout << "   [" << (*it)->point()[0] << "][" << (*it)->point()[1] << "][" << (*it)->point()[2] << "]" << std::endl;
                }
        }
        
}


void Snake::draw_snake2D() const // Pour dessiner le snake 2D (OpenGl)
{
        if(Snake_VH.size()==0)          
        {               
                std::cout << "WARNING in draw_snake: the snake is empty." << std::endl; 
        }       
        if(Snake_VH.size()==1)          
        {               
                std::cout << "WARNING in draw_snake: the snake contains only one point." << std::endl;  
        }       
        
        glClear(1);
        
        // draw lines
        glColor3f(0.8,0.2,0.2); // red
        glLineWidth(10.0);
        
        std::vector< Vertex_handle >::const_iterator it; 
        it=Snake_VH.begin();
        
        Point p1, p2;
        
        glDisable(GL_LIGHTING) ;
        
        glBegin(GL_LINES);
        
        p2==(*it)->point();
        while( (it!=Snake_VH.end()) && (++it!=Snake_VH.end()) ){
                p1=p2;
                p2=(*it)->point();
                // opengl
                glVertex3f(p1[0],p1[1],p1[2]);
                glVertex3f(p2[0],p2[1],p2[2]);
                
        }
        
        glEnd(); 
        glEnable(GL_LIGHTING) ;
        glFlush(); // pour la mise a jour
                
        glLineWidth(1.0); // on redonne une largeur normale a la ligne
}


void Snake::draw_snake3D() const // Pour dessiner le snake 3D (OpenGl)
{
        if(Snake_VH.size()==0)          
        {               
                std::cout << "WARNING in draw_snake: the snake is empty." << std::endl; 
        }       
        if(Snake_VH.size()==1)          
        {               
                std::cout << "WARNING in draw_snake: the snake contains only one point." << std::endl;  
        }       
        
        
        
        // draw lines
        glColor3f(0.8,0.2,0.2); // red
        glLineWidth(20.0);
        
        std::vector< Vertex_handle >::const_iterator it; 
        it=Snake_VH.begin();
        
        Point p1, p2;
        
        glDisable(GL_LIGHTING) ;
        
        glBegin(GL_LINES);
        
        p2==(*it)->point();
        R.Point3D_from_2D(p2);
        while( (it!=Snake_VH.end()) && (++it!=Snake_VH.end()) ){
                p1=p2;
                p2=(*it)->point();
                R.Point3D_from_2D(p2);
                // opengl
                glVertex3f(p1[0],p1[1],p1[2]);
                glVertex3f(p2[0],p2[1],p2[2]);
                
        }
        
        glEnd(); 
        glEnable(GL_LIGHTING) ;
        
        glFlush(); // pour la mise a jour
        
        glLineWidth(1.0); // on redonne une largeur normale a la ligne
}
                

// operationnel (teste) 
void Snake::Maj_A() // Maj de A en fonction des poids (attention a ne pas le faire pour rien)
{
        
        if ( WeightsA!= NULL && WeightsB!= NULL )
          {
                
        for(int i=0; i<dim; i++) //parcours des lignes de A
                  {
                for(int j=0; j<dim; j++) //parcours des colonnes de A
                          {
                                if (i == 0) 
                                        switch (j)
                                          {
                        case 0: 
                                                        A[i][j]= WeightsA[i]+WeightsA[i+1]+4*WeightsB[i]+WeightsB[i+1];
                                                        break;
                        case 1: 
                                                        A[i][j]=-WeightsA[i+1]-2*WeightsB[i]-2*WeightsB[i+1];
                                                        break;
                        case 2: 
                                                        A[i][j]=WeightsB[i+1];
                                                        break;
                        default:
                                                        A[i][j]=0 ;
                                                        break;
                                          }
                                                else if (i == dim-1)        
                                                        if (j == i-2)
                                A[i][j]=WeightsB[i-1];
                                                        else if (j==i-1)
                                A[i][j]=-WeightsA[i]-2*WeightsB[i-1]-2*WeightsB[i];
                                                        else if (j==i)
                                A[i][j]=WeightsA[i]+WeightsB[i-1]+4*WeightsB[i];
                                                        else A[i][j]=0;           
                                        else
                        if (j == i-2) 
                                                        A[i][j]=WeightsB[i-1];
                                        else if (j == i-1)
                                                A[i][j]=-WeightsA[i]-2*WeightsB[i-1]-2*WeightsB[i];
                                        else if (j == i)
                                                A[i][j]=WeightsA[i]+WeightsA[i+1]+WeightsB[i-1]+4*WeightsB[i]+WeightsB[i+1];
                                        else if (j == i+1)
                                                A[i][j]=-WeightsA[i+1]-2*WeightsB[i]-2*WeightsB[i+1];
                                        else if (j == i+2)
                                                A[i][j]=WeightsB[i+1];
                                        else A[i][j]=0;
                                
                          }
            
                  }   
                
          }else {
                  std::cout << "A cannot be modified without initializing WeightsA and WeightsB. " << std::endl; 
          }
        
}


// Modifieur et accesseur de la matrice A (pour des tests au cas ou...)
void Snake::Set_A(int i,int j, double f)
{
  A[i][j]=f;
}
        
void Snake::Affiche_A()
{
        if (A==NULL) std::cout << "La matrice A est vide." << std::endl;
        else {
                
                for (int i=0; i<dim; i++) // affichage de la ligne i
                  {
                        for (int j=0; j<dim; j++) std::cout << "["<< A[i][j] << "]";
                        
                        std::cout << std::endl;
                  }
                
        }
        
        
}

double Snake::Get_A(int i,int j)
{
  return A[i][j];
}
                
                                                        
                                                      
// Pour la mise a jour des points 2D du snake dans le processus de minimisation d'NRJ
// n'est adaptee que pour la parametrisation sur un carre de cote 6 unites
void Snake::Energy_minimization(int n_iterations)       
{   
        
        if (R_is_initialized) {
                
        double B[dim];  
        int count=0;
                
                // the current positions of the snake
                // of course we must have the same number n of 2 and 3D positions
                double Snake_PosX[dim];
                double Snake_PosY[dim];
                
                // pour le calcul de l'nrj feature des points du snake
                double fx[dim];  
                double fy[dim];
                // comment determiner l'nrj du premier point du snake?????????
                fx[0]=0;
                fy[0]=0;        
                
                Vertex_handle vh;
                
                std::vector< Vertex_handle >::const_iterator it; // notre iterator de vecteur
                for(it=Snake_VH.begin(); it!=Snake_VH.end(); it++)        
                  {
                        Snake_PosX[count++]=(*it)->point()[0];
                        Snake_PosY[count++]=(*it)->point()[1];
                        
                  }                                       
                
                
        for(int n=0; n<n_iterations; n++) // on a un nb d'iteration fixe pour l'evolution de notre snake dans le polyhedre 2D
                  {             
                        
                        //this->Affiche_snake(); // debug
                        
                        this->Maj_A(); // Calcul de la matrice A ici car A est modifiee par ResoudreSysteme (ou premiere it)
                        
                        // modif de A 
                        for (int i=0; i<dim; i++) this->Set_A(i, i, this->Get_A(i,i) + this->Get_Pas_Gamma() ); 
                        
                        
                        it=Snake_VH.begin(); 
                        // on calcule les nouvelles derivees partielles
                        for(int i=1; i<dim; i++) 
                          {
                                vh=*it; // vertex handle precedent
                                it++;                           
                                
                                if ( Snake_PosX[i]-Snake_PosX[i-1] != 0 ) 
                                        fx[i]=( R.Nrj_sommet(*it) - R.Nrj_sommet(vh) )/(Snake_PosX[i]-Snake_PosX[i-1]); 
                                
                                if ( Snake_PosY[i]-Snake_PosY[i-1] != 0 ) 
                                        fy[i]=( R.Nrj_sommet(*it) - R.Nrj_sommet(vh) )/(Snake_PosY[i]-Snake_PosY[i-1]); 
                                }
                        
                        // on prend en compte les dp
                        for(int i=0; i<dim; i++) Snake_PosX[i] += -fx[i];
                        
                        for(int i=0; i<dim; i++) Snake_PosY[i] += -fy[i];
                        
                        // on resoud notre systeme
                        if ( !ResoudreSysteme(A, Snake_PosX, B, dim) )
                          {
                                std::cerr << "ERROR IN ResoudreSysteme occured with Snake_PosX." << std::endl;
                                exit(-1);
                          }
                        
                        // maintenant il faut remettre a jour Snake_PosX                
                        for(int i=0; i<dim; i++) Snake_PosX[i] = B[i]; 
                        
                        this->Maj_A(); // Calcul de la matrice A ici car A est modifiee par ResoudreSysteme
                        
                        // modif de A 
                        for (int i=0; i<dim; i++) this->Set_A(i, i, this->Get_A(i,i) + this->Get_Pas_Gamma() );         
                        
                        
                        if ( !ResoudreSysteme(A, Snake_PosY, B, dim) )
                          {
                                std::cerr << "ERROR IN ResoudreSysteme occured with Snake_PosY." << std::endl;
                                exit(-1);        
                          }
                        
                        // maintenant il faut remettre a jour Snake_PosY                
                        for(int i=0; i<dim; i++) Snake_PosY[i] = B[i];
                        
                        
                        // on met a jour notre snake de telle sorte 
                        // que les points du snake reste dans le carre
                        count=0;
                        
                        for(it=Snake_VH.begin(); it!=Snake_VH.end(); it++)        
                          {     
                                if(0 <= Snake_PosX[count] && Snake_PosX[count]<=6 && 0 <= Snake_PosY[count] && Snake_PosY[count]<=6 )
                                {Point p(Snake_PosX[count],Snake_PosY[count++],0);
                                
                                (*it)->point() = p;
                                
                                }
                                
                                count++;
                          }
                        
                  }
                
                
        } else std::cerr << "The local region associated with the snake must be initialized to compute the energy minimization process." << std::endl; 
        
        
        
}
        
// Modifieur et accesseur du pas gamma
void Snake::Set_Pas_Gamma(double f)     
{
  Pas_Gamma=f;
}
        
double Snake::Get_Pas_Gamma()   
{
  return Pas_Gamma;
}
       
                
int Snake::Get_dim()
{
  return dim;
}        
                   
// Modifieur et accesseur des poids
void Snake::Set_Snake_WeightsA( double* Vec, int taille )       
{
  if (dim==taille) 
    {
      for(int i=0; i<dim; i++)
        {
          WeightsA[i]= Vec[i];      
        }            
    }
  else std::cerr << "ERROR IN Set_Snake_WeightsA( double* Vec, int taille ): taille must have the same value than dim. " << std::endl;
                          
}
        
void Snake::Set_Snake_WeightsB( double* Vec, int taille )
{
  if (dim==taille) 
    {
      for(int i=0; i<dim; i++)
        {
          WeightsB[i]= Vec[i];      
        }            
    }
  else std::cerr << "ERROR IN Set_Snake_WeightsB( double* Vec, int taille ): taille must have the same value than dim." << std::endl;
                
}
        
double* Snake::Get_Snake_WeightsA()     
{
  return WeightsA;
}
        
double* Snake::Get_Snake_WeightsB()
{
  return WeightsB;
}
            
            

---