triangulation.cpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Pascal Have

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

//  $Id: triangulation.cpp,v 1.6 2007/05/20 23:02:46 delpinux Exp $

/* 
T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   
Dès que Find_P_in_elt sera en quadtree virer tous les 'last_find' d'accélération relative
Dérécursivation de insertPoint et surtout scanNeigh
Optimiser la détection d'une frontière dans le coloriage
Mettre les lignes de contraintes à NULL quand les Quadtrees seront fonctionnels
T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O   T O  D O 
*/
#include <stack>
#include <fstream>
#include "triangulation.hpp"

#define ALTERNATE_INCIRCLE
#define DEBUG_FIND
//#define STEP_DETAILS
#define STEP_WRITE
#define STANDALONE

#ifdef STANDALONE 
template<class FLT>
inline FLT sqr(const FLT &x) { return x * x; }
#else /* STANDALONE */
//#include "timer.hpp"
//#include "communicator.hpp"
//#include "local_math.hpp"
//#include "utils.hpp"
#endif /* STANDALONE */

#include <cstdlib>

int myrand(const int width) {
    return (int)(static_cast<double>(width)*rand()/(RAND_MAX+1.0));
}

#include <map>
#include <set>
/************************************************************/
/******************** Fonctions locales *********************/
/************************************************************/

inline double SDet(const Vertex & V1, const Vertex & V2, const Vertex & V3) {
  return (V2[0]-V1[0]) * (V3[1]-V1[1]) -  (V2[1]-V1[1]) * (V3[0]-V1[0]);
}

// inline double SDet(const Vertex & V1, const Vertex & V2, const Vertex & V3) {
//   Triangle tr(&V1,&V2,&V3);
//   return Vertex((tr[1]-tr[0]) ^ (tr[2]-tr[0]))[2];
// }

inline bool checkInter(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) {
    // retourne true si P1P2 inter Q1Q2 != vide
    ASSERT(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2));
    return (SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0); // version optimisée
}

inline double Angle(const Vertex &V1, const Vertex &V2, const Vertex &V3) { // cos de l'angle * ||V2-V1|| (V1 point pivot)
  const double distance = sqrt(sqr(V3[0]-V1[0])+sqr(V3[1]-V1[1]));
  return (((V2[0]-V1[0])*(V3[0]-V1[0]))+((V2[1]-V1[1])*(V3[1]-V1[1])))/distance;
}

/************************************************************/
/********************** Triangulation ***********************/
/************************************************************/

void Triangulation::setBox(const PointNumType xmin, const PointNumType ymin,
                           const PointNumType xmax, const PointNumType ymax) {
  __triangles.clear(); // pour des bases saines
  __corners[0] = Vertex(xmin,ymax,0);
  __corners[1] = Vertex(xmin,ymin,0);
  __corners[2] = Vertex(xmax,ymin,0);
  __corners[3] = Vertex(xmax,ymax,0);
  __triangles.push_back(ConnectedTriangle(&__corners[0],&__corners[2],&__corners[3]));
  __triangles.push_back(ConnectedTriangle(&__corners[0],&__corners[1],&__corners[2]));
  __triangles.front().neigh(2) = (&__triangles.back());
  __triangles.back().neigh(1) = (&__triangles.front());

#ifdef STEP_WRITE
  { ASSERT(__checkAll()); std::ofstream o("init.mesh"); export_mesh(o); }
#endif /* STEP_WRITE */
}

#ifndef ALTERNATE_INCIRCLE
bool Triangulation::isInCircle(const Vertex & position, const Triangle & tri) const {
  double x1,x2,x3,y1,y2,y3;
  double a,b,c,d,e,f;
  double centre_cercle_x,centre_cercle_y;
  double rayon_cercle,distance;
  double tolerance = 1.e-13;
  x1=tri[0][0];
  x2=tri[1][0];
  x3=tri[2][0];
  y1=tri[0][1];
  y2=tri[1][1];
  y3=tri[2][1];

  a=x2-x3;
  b=y2-y3;
  c=((y2+y3)*(y2-y3)+(x2+x3)*(x2-x3))/2.;
  d=x2-x1;
  e=y2-y1;
  f=((y2+y1)*(y2-y1)+(x2+x1)*(x2-x1))/2.;
        
  centre_cercle_y = (d*c-a*f)/(b*d-a*e);
  centre_cercle_x = (c*e-f*b)/(a*e-b*d);
  //printf("centre du cercle :%g %g\n",(d*c-a*f)/(b*d-a*e),(c*e-f*b)/(a*e-b*d));
  distance = sqr(centre_cercle_x-position[0]) + sqr(centre_cercle_y-position[1]);
  rayon_cercle = sqr(centre_cercle_x-x1) + sqr(centre_cercle_y-y1);
  return (distance<=(rayon_cercle+tolerance));
}

#else /* ALTERNATE_INCIRCLE */

bool Triangulation::__isInCircle(const Vertex & position, const Triangle & tri) const {
  const double eps = 1e-13;
  // On peut éviter quelques indirections par des variables locales
  /* 9 affectations : 3+3+3*3 (+) et 0+0+3*2 (*) */
  const double
    BxmAx=tri(1)[0]-tri(0)[0],
    CxmAx=tri(2)[0]-tri(0)[0],
    DxmAx=position[0]-tri(0)[0];
  const double
    BymAy=tri(1)[1]-tri(0)[1],
    CymAy=tri(2)[1]-tri(0)[1],
    DymAy=position[1]-tri(0)[1];
  const double
    BA2=BxmAx*(tri(1)[0]+tri(0)[0])+BymAy*(tri(1)[1]+tri(0)[1]),
    CA2=CxmAx*(tri(2)[0]+tri(0)[0])+CymAy*(tri(2)[1]+tri(0)[1]),
    DA2=DxmAx*(position[0]+tri(0)[0])+DymAy*(position[1]+tri(0)[1]);
  const double r= // 5 (+) , 3*3 (*)
    (BxmAx*CymAy-CxmAx*BymAy)*DA2+(CxmAx*DymAy-DxmAx*CymAy)*BA2+(DxmAx*BymAy-BxmAx*DymAy)*CA2;
  // Total 20 (+) et 12 (*)
  return (r<eps); // ou <=0
}

#endif /* ALTERNATE_INCIRCLE */

void Triangulation::__permutation(TriangleIndex num_tri1, const unsigned k1) {
    // realise la permutation d'une arête commune a deux triangles.
    // TriangleIndex = pointeur du triangle.
    // k1 = numero du triangle voisin.
    // S = numero du sommet
    // V = numero des triangles voisins.

    //               S1                                         S1
    //         *                                                *
    //   V2   * *  V1                     V2   ***      V1
    //       *   *                                * * *
    //      *  T1 *                              *  *  *
    //   S2*********S4       S2 * T1*T2 *S4
    //      *  T2 *                              *  *  *
    //       *   *                                * * *
    //   V3   * *  V4                     V3   ***      V4
    //         *                                                *
    //                   S3                                 S3

    // => modification de V1 qui deivient T2 et V3 qui devient T1

    const TriangleIndex num_tri2 = num_tri1->neigh(k1);
    ASSERT(num_tri2 != NULL);
    const unsigned k2 = num_tri2->localizeNeigh(num_tri1);

    const PointIndex
        S1=num_tri1->base()[k1],
        S2=num_tri1->base()[(k1+1)%3],
        S4=num_tri1->base()[(k1+2)%3];

    const unsigned tag01 = num_tri1->tag;
    const unsigned tag02 = num_tri2->tag;
    //   say() << *num_tri1 << " " << k1 << " " << tag01 << " " << (tag01 & HSide[k1]) << "\n";
    //   say() << *num_tri2 << " " << k2 << " " << tag02 << " " << (tag02 & HSide[k2]) << "\n";

    ASSERT((tag01 & HSide[k1]) == 0); // Ne flip pas des arêtes marquées
    ASSERT((tag02 & HSide[k2]) == 0);
    num_tri1->tag = ((tag01&HSide[(k1+2)%3])?HSide[2]:0)|((tag02&HSide[(k2+1)%3])?HSide[0]:0);
    num_tri2->tag = ((tag01&HSide[(k1+1)%3])?HSide[1]:0)|((tag02&HSide[(k2+2)%3])?HSide[0]:0);

    const TriangleIndex
        V1=num_tri1->neigh((k1+1)%3),
        V2=num_tri1->neigh((k1+2)%3);

    PointIndex S3;
    TriangleIndex V3,V4;

    // recherche du sommet dans le triangle num_tri2 opposé au coté commun (ie opposé à k1 dans num_tri1)
    for(unsigned k=0;k<3;++k)
        if(num_tri2->neigh(k) == num_tri1) {
            S3=num_tri2->base()[k];
            V3=num_tri2->neigh((k+1)%3); // Ici, un bug est mort dans d'attroces souffrances
            V4=num_tri2->neigh((k+2)%3);
            break;
        }

            for(unsigned k=0;k<3;++k) {
                if(V1 != NULL)
                    if(V1->neigh(k) == num_tri1)
                        V1->neigh(k) = num_tri2;
                if(V3 != NULL)
                    if(V3->neigh(k) == num_tri2)
                        V3->neigh(k) = num_tri1;
            }

            num_tri1->setVertices(S1,S2,S3);
    num_tri2->setVertices(S1,S3,S4);
    num_tri1->setNeighs(V3,num_tri2,V2);
    num_tri2->setNeighs(V4,V1,num_tri1);
}

        
void Triangulation::__elimination_tri(TriangleIndex num_tri) {
  // elimination fu triangle pointe par neigh_index
  for(unsigned k=0;k<3;++k) {
    TriangleIndex neigh = num_tri->neigh(k);
    if(neigh != NULL) {
      for(unsigned k2=0;k2<3;++k2)
        if(neigh->neigh(k2) == num_tri) {
          neigh->neigh(k2) = NULL;
          break;
        }
    }
    num_tri->neigh(k) = NULL;
  }
  //  num_tri->tag = -1;
  __deleteTriangle(num_tri);
}


// recherche du triangle comportant l'element P
TriangleIndex Triangulation::__find_P_in_elt(const Vertex & P) {

#ifdef DEBUG_FIND
  // Version non optimisée de debug (en cas de trous de maillages supposé convexe)
  for(Triangles::iterator i = __triangles.begin();true;++i) {
    ASSERT(i != __triangles.end());
    const Vertex & P1 = (*i)(0);
    const Vertex & P2 = (*i)(1);
    const Vertex & P3 = (*i)(2);
    ASSERT(SDet(P1,P2,P3)>=0);
//    ffout(4)<<P1<<"\n";
//    ffout(4)<<P2<<"\n";
//    ffout(4)<<P3<<"\n";
//    ffout(4)<<"dans triangles?\n";
//    ffout(4)<<SDet(P1,P2,P)<<" "<<SDet(P2,P3,P)<<" "<<SDet(P3,P1,P)<<"\n";
    if ( (SDet(P1,P2,P)>=-1e-15) && (SDet(P2,P3,P)>=-1e-15) && (SDet(P3,P1,P)>=-1e-15) ) {
      return &*i;
    }
  }

#else /* DEBUG_FIND */

//   static TriangleIndex lastFind = NULL;
//   if (lastFind == NULL) lastFind = &*__triangles.begin();
  
  for(TriangleIndex i = &*__triangles.begin();;) {
    ASSERT(i != NULL);

    const Vertex & P2 = (*i)(1);
    const Vertex & P3 = (*i)(2);
    
    if(SDet(P2,P3,P)<0) {
      i=i->neigh(0);
    } else {
      const Vertex & P1 = (*i)(0); // on en a pas besoin avant
      if(SDet(P3,P1,P)<0) {
        i=i->neigh(1);
      } else {
        if(SDet(P1,P2,P)<0) {
          i=i->neigh(2);
        } else {
          // lastFind = i;  // Lié à Find_P_in_elt
          return i;
        }
      }
    }
  }

#endif /* DEBUG_FIND */
}


void Triangulation::__insertPoint(Vertex * newPoint, const TriangleIndex & curTriangle) {
  // triangle forcant cette position; évite les problèmes dus à l'éffacement de triangles
  __newTriangle(ConnectedTriangle(&__corners[0],&__corners[0],&__corners[0]));
  const Triangles::iterator first = --__triangles.end(); 

  Vertex * last = curTriangle->base()[1];
  for(unsigned i=0;i<3;++i) {
    if (__scanNeigh(curTriangle, i, newPoint, last) ) {
      // au bord de la bulle
      TriangleIndex tmpNeigh = curTriangle->neigh(i);
      __newTriangle(ConnectedTriangle(last,curTriangle->base()[(i+2)%3],newPoint,NULL,NULL,tmpNeigh));
      if (tmpNeigh != NULL)
        tmpNeigh->neigh(tmpNeigh->localizeNeigh(curTriangle)) = &__triangles.back();
      last = curTriangle->base()[(i+2)%3];
    }
  }
  __deleteTriangle(curTriangle);

  Triangles::iterator
    LTi = first,
    LTf = __triangles.end(),  
    LTaux;
  ++LTi;

  LTi->neigh(1) = &*(--LTf); // Cas particulier
  LTf->neigh(0) = &*LTi;

  LTaux = LTi;
  ++LTaux;

  while(LTi != LTf) { // ici LTf = --triangle.end()
    LTi->neigh(0) = &*LTaux;
    LTaux->neigh(1) = &*LTi;
    ++LTi; ++LTaux;
  }
  __triangles.erase(first); // plus utile 
}

bool Triangulation::__scanNeigh(const TriangleIndex & start,
                                const unsigned in,
                                Vertex * newPoint,
                                Vertex * &last) {
  // Scan à partir du triangle Start par le voisin in, à la recherche de triangles voisins contenant dans
  // leur cercle circonscrit le point newPoint; si OK met à jour la liste L entre LP1 et LP2
  const TriangleIndex curNeigh = start->neigh(in);

  if (curNeigh != NULL && __isInCircle(*newPoint,*curNeigh)) {
    const unsigned k = curNeigh->localizeNeigh(&*start); // sommet opposé au coté d'entrée
    
    if (__scanNeigh(curNeigh, (k+1)%3, newPoint, last) ) {
      // au bord de la bulle
      TriangleIndex tmpNeigh = curNeigh->neigh((k+1)%3);
      __newTriangle(ConnectedTriangle(last,curNeigh->base()[k],newPoint,NULL,NULL,tmpNeigh));
      if (tmpNeigh != NULL)
        tmpNeigh->neigh(tmpNeigh->localizeNeigh(curNeigh)) = &__triangles.back();
      last = curNeigh->base()[k];
    }

    // Idem que précédemment, mais pour pas de boucle, code dupliqué
    if (__scanNeigh(curNeigh, (k+2)%3, newPoint, last)) {
      // au bord de la bulle
      TriangleIndex tmpNeigh = curNeigh->neigh((k+2)%3);
      __newTriangle(ConnectedTriangle(last,curNeigh->base()[(k+1)%3],newPoint,NULL,NULL,tmpNeigh));
      if (tmpNeigh != NULL)
        tmpNeigh->neigh(tmpNeigh->localizeNeigh(curNeigh)) = &__triangles.back();
      last = curNeigh->base()[(k+1)%3];
    }
    __deleteTriangle(curNeigh);
    return false;
  }
  return true;
}

void Triangulation::__resetTags() {
  for(Triangles::iterator i = __triangles.begin(); i != __triangles.end(); ++i)
    i->tag = 0;
}


bool Triangulation::__checkLine(const CurveVertex &L, const bool closed) {
    ASSERT(L.size() > 1);
    ASSERT(closed); // !closed not implemented
    ASSERT(!closed || L.front() == L.back());

    typedef TinyVector<2,PointIndex> Edge;
    typedef std::map<Edge,TinyVector<2,TriangleIndex> > EdgeMapping;

    typedef std::set<Edge> BorderEdges;
    BorderEdges borderEdges;
    PointIndex lastPoint = L.front();

    for(CurveVertex::const_iterator i=L.begin(); ++i != L.end();) {
        Edge edge;
        edge[0]=lastPoint;
        edge[1]=*i;
        if(edge[0] > edge[1]) std::swap(edge[0],edge[1]); // arithmétique des pointeurs ou entiers
        borderEdges.insert(edge);
        lastPoint = *i;
    }

    EdgeMapping edgeMapping;
    do {
        edgeMapping.clear();
        for(Triangles::iterator tr = __triangles.begin();tr !=__triangles.end(); ++tr) {
            for(unsigned j=0; j<3; ++j) {
                PointIndex p1 = tr->base()[(j+1)%3];
                PointIndex p2 = tr->base()[(j+2)%3];
                if(p1 > p2) std::swap(p1,p2); // arithmétique des pointeurs ou des entiers
                Edge edge;
                edge[0]=p1;
                edge[1]=p2;
                const TriangleIndex tr2 = tr->neigh(j);
                
                if ((tr2 != NULL) && !(tr->tag & HSide[j]) && edgeMapping.find(edge) == edgeMapping.end()) {
                    Edge edgetemp;
                    edgetemp[0]=p1;
                    edgetemp[1]=p2;
                    TinyVector<2,TriangleIndex> trtemp;
                    trtemp[0]=&*tr;
                    trtemp[1]=tr2;
                    edgeMapping.insert(EdgeMapping::value_type(edgetemp,trtemp));
                }
            }
        }

        for(BorderEdges::const_iterator i = borderEdges.begin(); i != borderEdges.end();) {
EdgeMapping::iterator foundi = edgeMapping.find(*i);
            if (foundi != edgeMapping.end()) {
                // say() << "Edge " << *(*i)[0] << " " << *(*i)[1] << " found\n";
                TriangleIndex tr0 = foundi->second[0];
                TriangleIndex tr1 = foundi->second[1];
                ASSERT(tr0 != NULL && tr1 != NULL);
                tr0->tag |= HSide[tr0->localizeNeigh(tr1)];
                tr1->tag |= HSide[tr1->localizeNeigh(tr0)];
                edgeMapping.erase(foundi);
BorderEdges::const_iterator i_to_rm = i++;
                borderEdges.erase(i_to_rm);
            } else ++i;
        }

        if (!borderEdges.empty()) {
            bool done = false;
            while (!done) {
EdgeMapping::iterator e = edgeMapping.begin();
                for(unsigned n = myrand(edgeMapping.size());n>0;--n) ++e;
                const PointIndex P1 = e->first[0];
                const PointIndex P2 = e->first[1];
                TriangleIndex tr0 = e->second[0];
                TriangleIndex tr1 = e->second[1];
                ASSERT(tr0 != NULL && tr1 != NULL);
                const unsigned l0 = tr0->localizeNeigh(tr1);
                const unsigned l1 = tr1->localizeNeigh(tr0);
                const PointIndex Q1 = tr0->base()[l0];
                const PointIndex Q2 = tr1->base()[l1];
                if (checkInter(*P1,*P2,*Q1,*Q2)) {
                    // say() << "Flip edge " << *P1 << " " << *P2 << "\n";
                    done = true;
                    __permutation(tr0,tr0->localizeNeigh(tr1));
                }
            }
        }
    } while (!borderEdges.empty());

    return false;
}

bool Triangulation::__colorize(const Points & points) {
  // les variables pour les différents parcours
  std::stack<TriangleIndex> Stack;
  
  // On parcours les triangles pour propager la couleur de l'extérieur : 1
  // Recherche de l'élément initial : ie un triangle de coin de la box englobante
  Stack.push(__find_P_in_elt(__corners[0]));

  while(!Stack.empty()) {
    TriangleIndex itr = Stack.top(); Stack.pop();
    if (!(itr->tag & HMask)) { // pas de couleur
      itr->tag |= 1; // 1 : Couleur de remplissage
      // parcours des voisins (il faudrait peut être pas de boucle pour être plus rapide)
      for(unsigned j=0;j<3;j++) {
        TriangleIndex ktr = itr->neigh(j);
        if (ktr != NULL && // hors domaine
            !(itr->tag & HSide[j]) && 
            !(ktr->tag & HMask) // Déjà colorié
            ) {
          Stack.push(ktr); 
        }
      }
    }
  }

  Triangles::iterator llf = __triangles.begin();
  const Triangles::iterator lle = __triangles.end();
  while(llf != lle && !((llf->tag & HUMask) && 
                        !(llf->tag & HMask) &&
                        (llf->neigh(0)->tag & HMask ||
                         llf->neigh(1)->tag & HMask ||
                         llf->neigh(2)->tag & HMask))) {
    ++llf; 
  }

  if (llf == lle) {
#ifndef STEP_DETAILS
    ffout(4) << "Bad Domain\n"; // c'est pour cela qu'il faut commencer par colorier l'extérieur
#endif /* STEP_DETAILS */
    return false;
  }
  
  // On parcours les triangles pour propager une couleur
  // Initialisation
  Stack.push(&*llf);

  while(!Stack.empty()) {
    TriangleIndex itr = Stack.top();
    Stack.pop();
    if (!(itr->tag & HMask)) { // pas de couleur
        itr->tag |= 2; // colorisation du domaine avec la couleur 2
        // parcours des voisins (il faudrait peut être pas de boucle pour être plus rapide)               
        for(unsigned j=0;j<3;j++) {
          TriangleIndex ktr = itr->neigh(j);
          if (ktr != NULL && 
              !(itr->tag & HSide[j]) &&
              !(ktr->tag & HMask)) {
            Stack.push(ktr);
          }
        }
    }
  }
  return true;
}

TriangleIndex Triangulation::__findVertex(const Vertex * const p) {
  for(Triangles::iterator i = __triangles.begin(); i != __triangles.end(); ++i) {
    for(unsigned j=0;j<3;++j)
      if (i->base()[j] == p)
        return &*i;
  }
  ASSERT(false);
  return NULL;
}

bool Triangulation::triangulize(Points & points,
                                const CurveVertex & envVertex,
                                const std::vector<CurveVertex> & holes,
                                const std::vector<CurveVertex> & curves) {

    ASSERT(points.size() > 0);
  __points = &points;
  const Vertex & firstPoint = points[0];

  PointNumType
    xmin=firstPoint[0],xmax=firstPoint[0],
    ymin=firstPoint[1],ymax=firstPoint[1];

  for(unsigned i=0;i<points.size();++i) {
    const PointNumType x = points[i][0];
    const PointNumType y = points[i][1];
    xmin = std::min(xmin,x);
    xmax = std::max(xmax,x);
    ymin = std::min(ymin,y);
    ymax = std::max(ymax,y);
  }
  
  PointNumType
    dx=(xmax-xmin)/10,
    dy=(ymax-ymin)/10;

  setBox(xmin-2.*dx,ymin-dy,xmax+2.*dx,ymax+dy);
#ifdef STEP_DETAILS
  ffout(4) << "Box : (" << xmin-dx << " , " << ymin-dy << " ) , ( " << xmax+dx << " , " << ymax+dy << " )" << "\n";
#endif /* STEP_DETAILS */

#ifdef STEP_DETAILS
  //Timer timer(true);
#endif /* STEP_DETAILS */
  for(unsigned i=0;i<points.size();++i) {
    ASSERT(__checkAll());
    __findInsert(points[i]);
  }
#ifdef STEP_DETAILS
//  ffout(4) << "Insertion " << points.size() << " points [" << __triangles.size() << " triangles] in " << timer.//stop() << "_s" << "\n";
#endif /* STEP_DETAILS */

#ifdef STEP_WRITE
  { ASSERT(__checkAll()); std::ofstream o("insert.mesh"); export_mesh(o); }
#endif /* STEP_WRITE */

  // Important car c'est les CheckLine qui marque les cotés frontières dans les triangles
  __resetTags();

#ifdef STEP_DETAILS
  ffout(4) << "Check Enveloppe" << "\n";
  for(unsigned i=1;checkLine(envVertex,true);)
    
  {
      ffout(4) << "\tPasse " << ++i << "\n";
      ASSERT(checkAll()); std::ofstream o("passe.mesh"); export_mesh(o);
  }
#else /* STEP_DETAILS */
  while(__checkLine(envVertex,true));
#endif /* STEP_DETAILS */

//#ifdef STEP_DETAILS
//  ffout(4) << "Check Trou(s)" <<  "\n";
//  for(std::vector<CurveVertex>::iterator lf = holes.begin(); lf != holes.end(); ++lf)
//    for(unsigned i=1;checkLine(*lf,true);)
//      ffout(4) << "\tPasse " << ++i << "\n";
//#else /* STEP_DETAILS */
  for(std::vector<CurveVertex>::const_iterator lf = holes.begin(); lf != holes.end(); ++lf)
    while(__checkLine(*lf,true));
//#endif /* STEP_DETAILS */

#ifdef STEP_DETAILS
  ffout(4) << "Check Courbe(s)" << "\n";
  for(std::vector<CurveVertex>::const_iterator lf = curves.begin(); lf != curves.end(); ++lf)
    for(unsigned i=1;checkLine(*lf,false);)
      ffout(4) << "\tPasse " << ++i << "\n";
#else /* STEP_DETAILS */
  for(std::vector<CurveVertex>::const_iterator lf = curves.begin(); lf != curves.end(); ++lf)
    while(__checkLine(*lf,false));
#endif /* STEP_DETAILS */

#ifdef STEP_WRITE
  { ASSERT(__checkAll()); std::ofstream o("check.mesh"); export_mesh(o); }
#endif /* STEP_WRITE */

#ifdef STEP_DETAILS
  ffout(4) << "Colorize"  << "\n";
#endif /* STEP_DETAILS */
  __colorize(points);

#ifdef STEP_WRITE
  { ASSERT(__checkAll()); std::ofstream o("color.mesh"); export_mesh(o); }
#endif /* STEP_WRITE */

  unsigned count_tri = 0;
  const Triangles::iterator e = __triangles.end();
  for(Triangles::iterator ff,f = __triangles.begin(); f != e;) {
    ff = f++;
    if ((ff->tag & HMask) != 2) {
      __elimination_tri(&*ff); 
      ++count_tri;
    }
  }
#ifdef STEP_DETAILS
  ffout(4) << count_tri << " Triangles Deleted " << __triangles.size() << " triangles left" << "\n";
#endif /* STEP_DETAILS */

#ifdef STEP_WRITE
  { ASSERT(__checkAll()); std::ofstream o("final.mesh"); export_mesh(o); }
#endif /* STEP_WRITE */

  return true;
}

void Triangulation::export_mesh(std::ostream &o) const {
  const Points & points = *__points;

  o << "MeshVersionFormatted 1\n"
    "Dimension\n"
    "\t2\n"
    "Vertices\n"
    "\t" << points.size()+4
    << "\n";

  for(unsigned i=0;i<4;++i)
    o << __corners[i][0] << " " << __corners[i][1] << " 1\n";    
  for(Points::const_iterator i = points.begin(); i != points.end(); ++i) {
    o << i->x() << " " << i->y() << " 1\n";
  }

  class BestRef {
  private:
  public:
    const Vertex *points, *__corners;
  public:
//     BestRef(const Points * _points, const Point * ___corners) 
//       : points(_points), 
//      __corners(___corners) { }
    unsigned operator()(const Vertex & p) {
      const unsigned n = &p-__corners;
      if (n < 4)
        return n+1;
      else
        return &p-points+5;
    }
  } ref; // (firstPointer(points),firstPointer(__corners));
  ref.points = &points[0];
  ref.__corners = __corners;

  o << "Triangles\n\t" << __triangles.size() << "\n";
  for(Triangles::const_iterator i = __triangles.begin(); i != __triangles.end(); ++i) {
    o << " " << ref((*i)(0))
      << " " << ref((*i)(1))
      << " " << ref((*i)(2)) << " " << (i->tag & HMask) << "\n";
  }
}


bool Triangulation::__checkAll() {
  bool ok = true;
  
  unsigned jj = 0;
  for(Triangles::iterator i = __triangles.begin(); i != __triangles.end(); ++i,++jj) {
    // Check Orientation
    const double S = SDet((*i)(0),(*i)(1),(*i)(2));
    const bool check_point = 
      ((i->base())[0] != (i->base())[1]) && 
      ((i->base())[0] != (i->base())[2]) && 
      ((i->base())[1] != (i->base())[2]);
    
    if (!check_point) {
      ffout(4) << "Warning  : Triangle " << jj << " Plat : " /*<< *i */<< "\n";
      ok = false;
    }
    ASSERT(check_point);
    if (S <= 0) {
      ffout(4) << "< triangle " << jj << " " /*<< *i */<< " : " << S << "\n";
      ok = false;
    }

    // Check Neigh
    for(unsigned j=0;j<3;j++) {
      const TriangleIndex k = i->neigh(j);
      if (k != NULL) {
        const unsigned l = k->localizeNeigh(&*i);
        if (!(k->base()[(l+2)%3] == i->base()[(j+1)%3] && 
              k->base()[(l+1)%3] == i->base()[(j+2)%3])) {
          ffout(4) << "Bad edge " << j << " for triangle : " /*<< *i */<< "\n";
          ok = false;
        }
      }
    }
  }
  return ok;
}

---