/[ff3d]/ff3d/solver/triangulation.cpp
ViewVC logotype

Diff of /ff3d/solver/triangulation.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by dobj, Wed Oct 29 10:20:43 2003 UTC revision 1.3 by dobj, Tue Nov 4 15:35:27 2003 UTC
# Line 26  inline FLT sqr(const FLT &x) { return x Line 26  inline FLT sqr(const FLT &x) { return x
26  //#include "utils.hpp"  //#include "utils.hpp"
27  #endif /* STANDALONE */  #endif /* STANDALONE */
28    
29    #include <cstdlib>
30    
31    int myrand(const int width) {
32        return (int)(static_cast<double>(width)*rand()/(RAND_MAX+1.0));
33    }
34    
35    #include <map>
36    #include <set>
37  /************************************************************/  /************************************************************/
38  /******************** Fonctions locales *********************/  /******************** Fonctions locales *********************/
39  /************************************************************/  /************************************************************/
# Line 40  inline double SDet(const Vertex & V1, co Line 48  inline double SDet(const Vertex & V1, co
48  //   return Vertex((tr[1]-tr[0]) ^ (tr[2]-tr[0]))[2];  //   return Vertex((tr[1]-tr[0]) ^ (tr[2]-tr[0]))[2];
49  // }  // }
50    
51   inline bool checkInter(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) {  inline bool checkInter(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) {
52    // retourne true si P1P2 inter Q1Q2 != vide      // retourne true si P1P2 inter Q1Q2 != vide
53    assert(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2));      assert(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2));
54        return (SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0); // version optimisée
 //   const double  
 //     s1 = SDet(P1,P2,Q1)*SDet(P1,P2,Q2),  
 //     s2 = SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2);  
 //   if (s1 == 0 || s2 == 0) ffout(4) << "warning s=0\n";  
 //   return (s1 <= 0 && s2 <= 0);  
   
     const double eps = 1e-15;  
     if (std::abs(SDet(P1,P2,Q1)) < eps ||  
         std::abs(SDet(P1,P2,Q2)) < eps ||  
         std::abs(SDet(Q1,Q2,P1)) < eps ||  
         std::abs(SDet(Q1,Q2,P2)) < eps) {  
         //ffout(0) << " !!!!! NO !!!!!!!" << std::endl;  
         return false;  
     }  
       
   return (SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0); // version optimisée  
55  }  }
56    
57  inline double Angle(const Vertex &V1, const Vertex &V2, const Vertex &V3) { // cos de l'angle * ||V2-V1|| (V1 point pivot)  inline double Angle(const Vertex &V1, const Vertex &V2, const Vertex &V3) { // cos de l'angle * ||V2-V1|| (V1 point pivot)
# Line 143  bool Triangulation::isInCircle(const Ver Line 135  bool Triangulation::isInCircle(const Ver
135    
136  #endif /* ALTERNATE_INCIRCLE */  #endif /* ALTERNATE_INCIRCLE */
137    
138  void Triangulation::permutation(TriangleIndex num_tri1, const unsigned k) {  void Triangulation::permutation(TriangleIndex num_tri1, const unsigned k1) {
139    // realise la permutation d'une arête commune a deux triangles.      // realise la permutation d'une arête commune a deux triangles.
140    // TriangleIndex = pointeur du triangle.      // TriangleIndex = pointeur du triangle.
141    // k = numero du triangle voisin.      // k1 = numero du triangle voisin.
142    // S = numero du sommet      // S = numero du sommet
143    // V = numero des triangles voisins.      // V = numero des triangles voisins.
144      
145    //         S1                                         S1                            //               S1                                         S1
146    //         *                                          *                            //         *                                                *
147    //   V2   * *  V1                       V2   ***      V1            //   V2   * *  V1                     V2   ***      V1
148    //       *   *                                  * * *                      //       *   *                                * * *
149    //      *  T1 *                                *  *  *                      //      *  T1 *                              *  *  *
150    //   S2*********S4         S2 * T1*T2 *S4                  //   S2*********S4       S2 * T1*T2 *S4
151    //      *  T2 *                                *  *  *                      //      *  T2 *                              *  *  *
152    //       *   *                                  * * *                      //       *   *                                * * *
153    //   V3   * *  V4                       V3   ***      V4            //   V3   * *  V4                     V3   ***      V4
154    //         *                                          *      //         *                                                *
155    //             S3                                 S3                        //                   S3                                 S3
156    
157    // => modification de V1 qui deivient T2 et V3 qui devient T1      // => modification de V1 qui deivient T2 et V3 qui devient T1
158    
159        const TriangleIndex num_tri2 = num_tri1->neigh(k1);
160        assert(num_tri2 != NULL);
161        const unsigned k2 = num_tri2->localizeNeigh(num_tri1);
162    
163        const PointIndex
164            S1=num_tri1->base()[k1],
165            S2=num_tri1->base()[(k1+1)%3],
166            S4=num_tri1->base()[(k1+2)%3];
167    
168        const unsigned tag01 = num_tri1->tag;
169        const unsigned tag02 = num_tri2->tag;
170        //   say() << *num_tri1 << " " << k1 << " " << tag01 << " " << (tag01 & HSide[k1]) << "\n";
171        //   say() << *num_tri2 << " " << k2 << " " << tag02 << " " << (tag02 & HSide[k2]) << "\n";
172    
173        assert((tag01 & HSide[k1]) == 0); // Ne flip pas des arêtes marquées
174        assert((tag02 & HSide[k2]) == 0);
175        num_tri1->tag = ((tag01&HSide[(k1+2)%3])?HSide[2]:0)|((tag02&HSide[(k2+1)%3])?HSide[0]:0);
176        num_tri2->tag = ((tag01&HSide[(k1+1)%3])?HSide[1]:0)|((tag02&HSide[(k2+2)%3])?HSide[0]:0);
177    
178        const TriangleIndex
179            V1=num_tri1->neigh((k1+1)%3),
180            V2=num_tri1->neigh((k1+2)%3);
181    
182        PointIndex S3;
183        TriangleIndex V3,V4;
184    
185        // recherche du sommet dans le triangle num_tri2 opposé au coté commun (ie opposé à k1 dans num_tri1)
186        for(unsigned k=0;k<3;++k)
187            if(num_tri2->neigh(k) == num_tri1) {
188                S3=num_tri2->base()[k];
189                V3=num_tri2->neigh((k+1)%3); // Ici, un bug est mort dans d'attroces souffrances
190                V4=num_tri2->neigh((k+2)%3);
191                break;
192            }
193    
194    const TriangleIndex num_tri2 = num_tri1->neigh(k);              for(unsigned k=0;k<3;++k) {
195    assert(num_tri2 != NULL);                  if(V1 != NULL)
196                          if(V1->neigh(k) == num_tri1)
197    PointIndex                          V1->neigh(k) = num_tri2;
198      S1=num_tri1->base()[k],                  if(V3 != NULL)
199      S2=num_tri1->base()[(k+1)%3],                      if(V3->neigh(k) == num_tri2)
200      S4=num_tri1->base()[(k+2)%3];                          V3->neigh(k) = num_tri1;
201                }
   TriangleIndex  
     V1=num_tri1->neigh((k+1)%3),  
     V2=num_tri1->neigh((k+2)%3);  
     
   PointIndex S3;  
   TriangleIndex V3,V4;  
202    
203    // recherche du sommet dans le triangle num_tri2 opposé au coté commun (ie opposé à k dans num_tri1)              num_tri1->setVertices(S1,S2,S3);
204    for(unsigned k=0;k<3;++k)      num_tri2->setVertices(S1,S3,S4);
205      if(num_tri2->neigh(k) == num_tri1) {      num_tri1->setNeighs(V3,num_tri2,V2);
206        S3=num_tri2->base()[k];      num_tri2->setNeighs(V4,V1,num_tri1);
       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);  
207  }  }
208    
209                    
210  void Triangulation::elimination_tri(TriangleIndex num_tri) {  void Triangulation::elimination_tri(TriangleIndex num_tri) {
211    // elimination fu triangle pointe par neigh_index    // elimination fu triangle pointe par neigh_index
# Line 352  void Triangulation::resetTags() { Line 356  void Triangulation::resetTags() {
356      i->tag = 0;      i->tag = 0;
357  }  }
358    
359    
360  bool Triangulation::checkLine(const CurveVertex &L, const bool closed) {  bool Triangulation::checkLine(const CurveVertex &L, const bool closed) {
361    assert(L.size() > 1);      assert(L.size() > 1);
362    if (closed)      assert(closed); // !closed not implemented
363      assert(L.front() == L.back());      assert(!closed || L.front() == L.back());
364    
365        typedef TinyVector<2,PointIndex> Edge;
366        typedef std::map<Edge,TinyVector<2,TriangleIndex> > EdgeMapping;
367    
368        typedef std::set<Edge> BorderEdges;
369        BorderEdges borderEdges;
370        PointIndex lastPoint = L.front();
371    
372        for(CurveVertex::const_iterator i=L.begin(); ++i != L.end();) {
373            Edge edge;
374            edge[0]=lastPoint;
375            edge[1]=*i;
376            if(edge[0] > edge[1]) std::swap(edge[0],edge[1]); // arithmétique des pointeurs ou entiers
377            borderEdges.insert(edge);
378            lastPoint = *i;
379        }
380    
381    bool change = false;      EdgeMapping edgeMapping;
382        do {
383            edgeMapping.clear();
384            for(Triangles::iterator tr = triangles.begin();tr !=triangles.end(); ++tr) {
385                for(unsigned j=0; j<3; ++j) {
386                    PointIndex p1 = tr->base()[(j+1)%3];
387                    PointIndex p2 = tr->base()[(j+2)%3];
388                    if(p1 > p2) std::swap(p1,p2); // arithmétique des pointeurs ou des entiers
389                    Edge edge;
390                    edge[0]=p1;
391                    edge[1]=p2;
392                    const TriangleIndex tr2 = tr->neigh(j);
393                    
394                    if ((tr2 != NULL) && !(tr->tag & HSide[j]) && edgeMapping.find(edge) == edgeMapping.end()) {
395                        Edge edgetemp;
396                        edgetemp[0]=p1;
397                        edgetemp[1]=p2;
398                        TinyVector<2,TriangleIndex> trtemp;
399                        trtemp[0]=&*tr;
400                        trtemp[1]=tr2;
401                        edgeMapping.insert(EdgeMapping::value_type(edgetemp,trtemp));
402                    }
403                }
404            }
405    
406    CurveVertex::const_iterator position = L.begin();          for(BorderEdges::const_iterator i = borderEdges.begin(); i != borderEdges.end();) {
407    const CurveVertex::const_iterator end = L.end();  EdgeMapping::iterator foundi = edgeMapping.find(*i);
408                  if (foundi != edgeMapping.end()) {
409    PointIndex                  // say() << "Edge " << *(*i)[0] << " " << *(*i)[1] << " found\n";
410      curVertex = *position++,                  TriangleIndex tr0 = foundi->second[0];
411      curTarget = *position;                  TriangleIndex tr1 = foundi->second[1];
412                    assert(tr0 != NULL && tr1 != NULL);
413    // Etape initiale : recherche d'un triangle possédant en sommet le premier point de L                  tr0->tag |= HSide[tr0->localizeNeigh(tr1)];
414    // On peut y mettre une recherche quadtree ou par adjacente optimale (cf find_P_in_elt)                  tr1->tag |= HSide[tr1->localizeNeigh(tr0)];
415    TriangleIndex curTriangle = findVertex(curVertex);                  edgeMapping.erase(foundi);
416    BorderEdges::const_iterator i_to_rm = i++;
417    bool found = false;                  borderEdges.erase(i_to_rm);
418    while (!found) {              } else ++i;
419      assert(curTriangle->base()[0] == curVertex ||          }
420             curTriangle->base()[1] == curVertex ||  
421             curTriangle->base()[2] == curVertex);          if (!borderEdges.empty()) {
422                bool done = false;
423      unsigned localVertex=curTriangle->localizeVertex(curVertex); // position local de curVertex              while (!done) {
424    EdgeMapping::iterator e = edgeMapping.begin();
425      if (curTriangle->base()[(localVertex+1)%3] == curTarget) {                  for(unsigned n = myrand(edgeMapping.size());n>0;--n) ++e;
426        // On a trouvé curTarget dans le triangle courant, on passe au suivant.                  const PointIndex P1 = e->first[0];
427        // On marque le coté comme frontière si closed == true ie c'est pas une ligne interne                  const PointIndex P2 = e->first[1];
428        if (closed) {                  TriangleIndex tr0 = e->second[0];
429          // coté entre localVertex et LocalVertex+1 => neigh[localVertex+2]                  TriangleIndex tr1 = e->second[1];
430          curTriangle->tag |= HSide[(localVertex+2)%3];                  assert(tr0 != NULL && tr1 != NULL);
431          TriangleIndex NTri = curTriangle->neigh((localVertex+2)%3); // forcement non null                  const unsigned l0 = tr0->localizeNeigh(tr1);
432          assert(NTri != NULL); // Car tout est inclus dans la boite                  const unsigned l1 = tr1->localizeNeigh(tr0);
433          NTri->tag |= HSide[NTri->localizeNeigh(curTriangle)];                  const PointIndex Q1 = tr0->base()[l0];
434        }                  const PointIndex Q2 = tr1->base()[l1];
435        curVertex = curTarget;                  if (checkInter(*P1,*P2,*Q1,*Q2)) {
436        found = (++position == end);                      // say() << "Flip edge " << *P1 << " " << *P2 << "\n";
437        curTarget = *position;                      done = true;
438                        permutation(tr0,tr0->localizeNeigh(tr1));
439      } else if (curTriangle->base()[(localVertex+2)%3] == curTarget) {                  }
       // On a trouvé curTarget dans le triangle courant, on passe au suivant.  
       // On marque le coté comme frontière si closed == true ie c'est pas une ligne interne  
       if (closed) {  
         // coté entre localVertex et LocalVertex+1 => neigh[localVertex+2]  
         curTriangle->tag |= HSide[(localVertex+1)%3];  
         TriangleIndex NTri = curTriangle->neigh((localVertex+1)%3); // forcement non null  
         assert(NTri != NULL); // Car tout est inclus dans la boite  
         NTri->tag |= HSide[NTri->localizeNeigh(curTriangle)];  
       }  
       curVertex = curTarget;  
       found = (++position == end);  
       curTarget = *position;  
   
     } else if(checkInter(*curVertex,*curTarget,  
                          (*curTriangle)((localVertex+1)%3),  
                          (*curTriangle)((localVertex+2)%3))) {  
       // Ok : found a problem on the border  
       // gestion spécifique du cnflit arète-frontière  
       bool atteint = false;  
       while(!atteint) {  
         assert(curTriangle != NULL);  
             
         // description globale du coté opposée à localvertex dans curTriangle ie l'entrée  
         const PointIndex  
           g1 = curTriangle->base()[(localVertex+1)%3],  
           g2 = curTriangle->base()[(localVertex+2)%3];  
             
         TriangleIndex new_t = curTriangle->neigh(localVertex);  
         assert(new_t != NULL);  
             
         unsigned l=new_t->localizeNeigh(curTriangle);  
             
         const PointIndex from = curTriangle->base()[localVertex];  
         const PointIndex to   = new_t->base()[l];  
             
         // Test de convexité  
         if (checkInter(*from,*to,*g1,*g2)) { // c'est convexe ou fusinnable  
 #ifndef STEP_DETAILS  
           // Test de fusion  
           // g1 sur [from,to]  
           if (SDet(*from,*to,*g1) == 0)  
             ffout(4) << "G1 is on the border\n";  
           // g2 sur [from,to]  
           if (SDet(*from,*to,*g2) == 0)  
             ffout(4) << "G2 is on the border\n";                
           // from sur [g1,g2]  
           if (SDet(*g1,*g2,*from) == 0)  
             ffout(4) << "from is on the border\n";  
           // tp sur [g1,g2]  
           if (SDet(*g1,*g2,*to) == 0)  
             ffout(4) << "to is on the border\n";  
 #endif /* STEP_DETAILS */  
             
           permutation(curTriangle,localVertex);  
           change = true; // changement  
   
           // le cas où g3 serait curTarget  
           localVertex=curTriangle->localizeVertex(from);  
           l=new_t->localizeVertex(from);  
   
           if (to == curTarget) { // passer au point suivant.  
 #warning !  
 //          if (closed && from == curVertex) {  
 //            // on a l'arête complète  
 //            curTriangle->tag |= HSide[curTriangle->localizeNeigh(new_t)];  
 //            new_t->tag |= HSide[new_t->localizeNeigh(curTriangle)];  
 //          }  
 #warning end of !  
             curVertex = curTarget;  
             found = (++position == end);  
             curTarget = *position;  
             atteint = true;  
           } else {  
             // on pourrait y mettre les noms globaux que l'on y connait  
             if (checkInter(*curVertex,*curTarget,  
                            (*curTriangle)((localVertex+1)%3),  
                            (*curTriangle)((localVertex+2)%3))) {  
               // curTriangle : OK  
               // localVertex : OK  
             } else {  
               curTriangle = new_t;  
               localVertex = l;  
               assert(checkInter(*curVertex,*curTarget,  
                                 (*curTriangle)((l+1)%3),  
                                 (*curTriangle)((l+2)%3)));  
440              }              }
           }  
         } else {  
           // c'est pas convexe et ça coupe encore l'arète  
           curTriangle = new_t;  
           // test de continuation  
           if (checkInter(*curVertex,*curTarget,  
                          (*curTriangle)(l),  
                          (*curTriangle)((l+1)%3))) {  
             // sinon c'est l'autre  
             // c'est sur la droite [l, l+1] opp à l+2  
             localVertex = (l+2)%3;  
           } else {  
             // c'est sur la droite [l,l+2] opp à l+1  
             assert(checkInter(*curVertex,*curTarget,  
                               (*curTriangle)(l),  
                               (*curTriangle)((l+2)%3)));  
             localVertex = (l+1)%3;  
           }  
441          }          }
442        }      } while (!borderEdges.empty());
     } else {  
       // A fortment optimiser, peut coupable avec le calcul d'intersection        
       const double theta1 = Angle(*curVertex,*curTarget,  
                                   (*curTriangle)((localVertex+1)%3));  
       const double theta2 = Angle(*curVertex,*curTarget,  
                                   (*curTriangle)((localVertex+2)%3));  
       const int sens = (theta1 > theta2)?2:1;  
       //  ffout(4) << "angles : " << theta1 << " " << theta2 << " sens : " << sens << "\n";      
       curTriangle = curTriangle->neigh((localVertex+sens)%3);  
     }  
   } // end of while(!found)  
443    
444    return change;      return false;
445  }  }
446    
447  bool Triangulation::colorize(const Points & points) {  bool Triangulation::colorize(const Points & points) {

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26