boolops_triangulation.hpp

Go to the documentation of this file.

// Copyright (c) 2012-2019 University of Lyon and CNRS (France).
// All rights reserved.
//
// This file is part of MEPP2; 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 3 of the License,
// or (at your option) any later version.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 
#pragma once
 
#include <CGAL/Triangulation_2.h>
#include <CGAL/Constrained_Delaunay_triangulation_2.h>
 
template <class Gt, class Vb = CGAL::Triangulation_vertex_base_2<Gt> >
class Enriched_vertex_base : public Vb
{
public:
  template < typename TDS2 >
  struct Rebind_TDS {
    typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2;
    typedef Enriched_vertex_base<Gt,Vb2> Other;
  };
 
private:
  unsigned long m_Label;
 
public:
  void set_Label(unsigned long Label) {m_Label = Label;}
  unsigned long get_Label() {return m_Label;}
};
 
template <class Gt, class Fb = CGAL::Constrained_triangulation_face_base_2<Gt> >
class Enriched_face_base : public Fb
{
public:
  typedef typename Fb::Triangulation_data_structure::Vertex_handle          Vertex_handle;
 
  typedef typename Fb::Triangulation_data_structure::Face_handle            Face_handle;
 
  template < typename TDS2 >
  struct Rebind_TDS {
    typedef typename Fb::template Rebind_TDS<TDS2>::Other  Fb2;
    typedef Enriched_face_base<Gt,Fb2>            Other;
  };
 
private:
  bool m_OK;
  bool m_Ext;
 
public:
  Enriched_face_base() : Fb() {}
  Enriched_face_base(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2) : Fb(v0,v1,v2) {}
  Enriched_face_base(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2,
    Face_handle n0, Face_handle n1, Face_handle n2) : Fb(v0,v1,v2,n0,n1,n2) {}
  Enriched_face_base(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Face_handle n0, Face_handle n1, Face_handle n2,
          bool c0, bool c1, bool c2 ) : Fb(v0,v1,v2,n0,n1,n2) {}
  void set_Ext(bool Ext) {m_Ext = Ext;}
  bool get_Ext() {return m_Ext;}
  void set_OK(bool OK) {m_OK = OK;}
  bool get_OK() {return m_OK;}
};
 
template <class K>
class Triangulation
{
  typedef typename K::Point_3                            Point_3;
 
        typedef /*typename*/ Enriched_vertex_base<K>                    Tri_vb;
 
        typedef /*typename*/ Enriched_face_base<K>                      Tri_fb;
 
  typedef typename CGAL::Triangulation_data_structure_2<Tri_vb,Tri_fb>      Tri_DS;
 
//DBG #define STR(x) #x
//DBG #define XSTR(x) STR(x)
//DBG #pragma message("CGAL_VERSION_NR=" XSTR(CGAL_VERSION_NR))
//DBG #pragma message("CGAL_VERSION_STR=" XSTR(CGAL_VERSION_STR))
 
// Note: the following test should be
//          #if( CGAL_VERSION_NR >= CGAL_VERSION_NUMBER(5,1,0) )
//       but CGAL_VERSION_NUMBER(5,1,0)=1050101000 whereas the version of CGAL
//       packaged with vcpkg has CGAL_VERSION_NR=1050100900 (it seems that it
//       is not a public release).
#if( CGAL_VERSION_NR >= 1050100000 )
  typedef typename CGAL::No_constraint_intersection_requiring_constructions_tag
      Itag;
#else
  typedef typename CGAL::No_intersection_tag                    Itag;
#endif
 
  typedef typename CGAL::Constrained_Delaunay_triangulation_2<K, Tri_DS, Itag>  Constrained_Delaunay_tri;
 
  typedef typename Constrained_Delaunay_tri::Vertex_handle            Vertex_handle_tri;
 
  typedef typename Constrained_Delaunay_tri::Face_handle              Face_handle_tri;
 
  typedef typename Constrained_Delaunay_tri::Point                Point_tri;
 
  typedef typename Constrained_Delaunay_tri::Face_iterator            Face_iterator_tri;
 
public:
  Triangulation(Halfedge_handle &he, Vector_exact &norm_dir)
  {
    //find the longest coordinate of the normal vector, and its sign
    double x = to_double(norm_dir.x());
    double y = to_double(norm_dir.y());
    double z = to_double(norm_dir.z());
    double absx = std::abs(x);
    double absy = std::abs(y);
    double absz = std::abs(z);
 
    //this information is stored using a code :
    //0 : The coordinate X is the longest, and positive
    //1 : The coordinate Y is the longest, and positive
    //2 : The coordinate Z is the longest, and positive
    //3 : The coordinate X is the longest, and negative
    //4 : The coordinate Y is the longest, and negative
    //5 : The coordinate Z is the longest, and negative
    if (absx >= absy && absx >= absz) max_coordinate = (x>0)?0:3;
    else if (absy >= absx && absy >= absz) max_coordinate = (y>0)?1:4;
    else if (absz >= absx && absz >= absy) max_coordinate = (z>0)?2:5;
 
    //we add the three vertices of the facet to the triangulation
    //The Label of these vertices is set for the corresponding point in the triangulation
    v1 = add_new_pt(point_to_exact(he->vertex()->point()), he->vertex()->Label);
    v2 = add_new_pt(point_to_exact(he->next()->vertex()->point()), he->next()->vertex()->Label);
    v3 = add_new_pt(point_to_exact(he->next()->next()->vertex()->point()), he->next()->next()->vertex()->Label);
 
    //if the vertices does not have an Id (Label = OxFFFFFFFF), the labels
    //of the points in the triangulation is set as follows :
    //0xFFFFFFFF for the first point
    //0xFFFFFFFE for the second point
    //0xFFFFFFFD for the third point
    if(v2->get_Label() == 0xFFFFFFFF) v2->set_Label(0xFFFFFFFE);
    if(v3->get_Label() == 0xFFFFFFFF) v3->set_Label(0xFFFFFFFD);
  }
 
  Point_tri get_minvar_point_2(Point_3 &p)
  {
    switch(max_coordinate)
    {
    case 0:
      return Point_tri(p.y(),p.z());
      break;
    case 1:
      return Point_tri(p.z(),p.x());
      break;
    case 2:
      return Point_tri(p.x(),p.y());
      break;
    case 3:
      return Point_tri(p.z(),p.y());
      break;
    case 4:
      return Point_tri(p.x(),p.z());
      break;
    case 5:
      return Point_tri(p.y(),p.x());
      break;
    default:
      return Point_tri(p.y(),p.z());
    }
  }
 
        Vertex_handle_tri add_new_pt(Point_3 p, unsigned long &Label)   // MT: suppression r�f�rence
  {
    //if the point is not a new one, we verify that the point has not already been added
    if(Label != 0xFFFFFFFF)
      for(unsigned int i = 0;i != pts_point.size();++i)
        if(Label == pts_point[i])
          //if the point is already in the triangulation, we return its handle
          return pts_vertex[i];
    Vertex_handle_tri v;
    v = ct.insert(get_minvar_point_2(p));
    v->set_Label(Label);
    pts_point.push_back(Label);
    pts_vertex.push_back(v);
    return v;
  }
 
  void add_segment(Point_3 &p1, Point_3 &p2, unsigned long &Label1, unsigned long &Label2)
  {
    // we add the two points in the triangulation and store their handles in c1 and c2
    c1 = add_new_pt(p1, Label1);
    c2 = add_new_pt(p2, Label2);
    // we set a constrained segment between these two points
    ct.insert_constraint(c1, c2);
    // if an other segment is added, we will overwrite c1 and c2
    // what is important is to memorize the handles of the last segment added
  }
 
  std::vector< std::vector<unsigned long> > get_triangles(bool inv_triangles, bool *IsExt)
  {
    //init
    IsExt[0] = false;
    IsExt[1] = false;
    IsExt[2] = false;
    std::vector<std::vector<unsigned long> > tris;
    for(Face_iterator_tri fi = ct.faces_begin();fi != ct.faces_end();fi++)
      fi->set_OK(false);
 
    //the constrained segments are oriented, we search the triangle (c1, c2, X), (X, c1, c2) or (c2, X, c1)
    //where c1 and c2 are the two points related to the last constrained segment added, and X another point
    //this triangle belongs to the result (thanks to the orientation of the segments)
    Face_handle_tri f, f2 = c1->face();
                int i=0; // MT
    do {
      f = f2;
      f->has_vertex(c1,i);
      f2 = f->neighbor(f->ccw(i));
    } while( ! ( f->has_vertex(c2) && f2->has_vertex(c2) ) );
 
    //dans le cas particulier ou la frontiere se trouve exactement sur un bord de la triangulation,
    //et que ce triangle n'appartient pas a la triangulation, on d�marrera avec l'autre triangle
    //incluant le segment c1, c2 (et donc, n'appartenant pas au r�sultat
 
    //if the segment is exactly on the border of the triangulation, the triangle could be outside the triangulation
    //in that case, we will search the other triangle including the points c1 and c2
    //this triangle does not belong to the result
    if(f->has_vertex(ct.infinite_vertex()))
    {
      f = f2;
      f->set_Ext(false);
    }
    else
    {
      f->set_Ext(true);
    }
 
    std::stack<Face_handle_tri> sfh;
    f->set_OK(true);
    sfh.push(f);
 
    //we decide for all the triangles, if they belongs to the result, starting from the first triangle f,
    //by moving on the triangulation using the connectivity between the triangles.
    //If a constrained segment is crossed, the value of the tag "isext" is inverted
    while(!sfh.empty())
    {
      f = sfh.top();
      sfh.pop();
 
      if(f->get_Ext())
      {
        std::vector<unsigned long> tri;
        int i;
        tri.push_back(f->vertex(0)->get_Label());
 
        //verify if the neighboring facets belongs to the result or not
        if(f->has_vertex(v1,i) && f->neighbor(f->ccw(i))->has_vertex(ct.infinite_vertex())) IsExt[0] = true;
        if(f->has_vertex(v2,i) && f->neighbor(f->ccw(i))->has_vertex(ct.infinite_vertex())) IsExt[1] = true;
        if(f->has_vertex(v3,i) && f->neighbor(f->ccw(i))->has_vertex(ct.infinite_vertex())) IsExt[2] = true;
 
        if(inv_triangles)
        {
          tri.push_back(f->vertex(2)->get_Label());
          tri.push_back(f->vertex(1)->get_Label());
        }
        else
        {
          tri.push_back(f->vertex(1)->get_Label());
          tri.push_back(f->vertex(2)->get_Label());
        }
        tris.push_back(tri);
      }
      for(i = 0;i!=3;i++)
      {
        if(!(f->neighbor(i)->get_OK() || f->neighbor(i)->has_vertex(ct.infinite_vertex())))
        {
          f->neighbor(i)->set_OK(true);
          f->neighbor(i)->set_Ext((f->is_constrained(i))?!f->get_Ext():f->get_Ext());
          sfh.push(f->neighbor(i));
        }
      }
    }
    return tris;
  }
 
  std::vector< std::vector<unsigned long> > get_all_triangles(bool inv_triangles)
  {
    std::vector< std::vector<unsigned long> > tris;
    for(Face_iterator_tri f = ct.faces_begin();f != ct.faces_end();f++)
    {
      std::vector<unsigned long> tri;
      tri.push_back(f->vertex(0)->get_Label());
      if(inv_triangles)
      {
        tri.push_back(f->vertex(2)->get_Label());
        tri.push_back(f->vertex(1)->get_Label());
      }
      else
      {
        tri.push_back(f->vertex(1)->get_Label());
        tri.push_back(f->vertex(2)->get_Label());
      }
      tris.push_back(tri);
    }
    return tris;
  }
 
private:
  Constrained_Delaunay_tri ct;
  std::vector<InterId> pts_point;
  std::vector<Vertex_handle_tri> pts_vertex;
  Vertex_handle_tri v1;
  Vertex_handle_tri v2;
  Vertex_handle_tri v3;
  Vertex_handle_tri c1;
  Vertex_handle_tri c2;
  int max_coordinate;
};