self_intersect.h

  1 // compute self-intersection of a CGAL triangle polyhedron mesh
  2 // original code from Lutz Kettner
  3 #ifndef _SELF_INTERSECT_H_
  4 #define _SELF_INTERSECT_H_
  5 
  6 #include <CGAL/box_intersection_d.h>
  7 #include <CGAL/intersections.h>
  8 #include <CGAL/Bbox_3.h>
  9 
 10 template <class Polyhedron, class Kernel, class OutputIterator>
 11 struct Intersect_facets
 12 {
 13   typedef typename Kernel::Point_3      Point;
 14   typedef typename Kernel::Vector_3     Vector;
 15   typedef typename Kernel::Segment_3    Segment;
 16   typedef typename Kernel::Triangle_3   Triangle;
 17   typedef typename Polyhedron::Facet_const_handle       Facet_const_handle;
 18   typedef typename Polyhedron::Facet_const_iterator     Facet_const_iterator;
 19   typedef typename Polyhedron::Halfedge_const_handle    Halfedge_const_handle;
 20   typedef typename CGAL::Box_intersection_d::Box_with_handle_d<double, 3, Facet_const_handle> Box;
 21   mutable OutputIterator m_iterator;
 22 
 23 public:
 24 
 25   Intersect_facets(OutputIterator it)
 26     : m_iterator(it)
 27   {
 28   }
 29 
 30   void operator()(const Box* b,
 31     const Box* c) const
 32   {
 33     Halfedge_const_handle h = b->handle()->halfedge();
 34 
 35     // check for shared egde --> no intersection
 36     if(h->opposite()->facet() == c->handle() ||
 37       h->next()->opposite()->facet() == c->handle() ||
 38       h->next()->next()->opposite()->facet() == c->handle())
 39       return;
 40 
 41     // check for shared vertex --> maybe intersection, maybe not
 42     Halfedge_const_handle g = c->handle()->halfedge();
 43     Halfedge_const_handle v;
 44 
 45     if(h->vertex() == g->vertex())
 46       v = g;
 47     if(h->vertex() == g->next()->vertex())
 48       v = g->next();
 49     if(h->vertex() == g->next()->next()->vertex())
 50       v = g->next()->next();
 51 
 52     if(v == Halfedge_const_handle())
 53     {
 54       h = h->next();
 55       if(h->vertex() == g->vertex())
 56     v = g;
 57       if(h->vertex() == g->next()->vertex())
 58     v = g->next();
 59       if(h->vertex() == g->next()->next()->vertex())
 60     v = g->next()->next();
 61       if(v == Halfedge_const_handle())
 62       {
 63     h = h->next();
 64     if(h->vertex() == g->vertex())
 65       v = g;
 66     if(h->vertex() == g->next()->vertex())
 67       v = g->next();
 68     if(h->vertex() == g->next()->next()->vertex())
 69       v = g->next()->next();
 70       }
 71     }
 72 
 73     if(v != Halfedge_const_handle())
 74     {
 75       // found shared vertex: 
 76       CGAL_assertion(h->vertex() == v->vertex());
 77       // geometric check if the opposite segments intersect the triangles
 78       Triangle t1( h->vertex()->point(),
 79     h->next()->vertex()->point(),
 80     h->next()->next()->vertex()->point());
 81       Triangle t2( v->vertex()->point(),
 82     v->next()->vertex()->point(),
 83     v->next()->next()->vertex()->point());
 84       Segment s1( h->next()->vertex()->point(),
 85     h->next()->next()->vertex()->point());
 86       Segment s2( v->next()->vertex()->point(),
 87     v->next()->next()->vertex()->point());
 88 
 89       if(CGAL::do_intersect(t1,s2))
 90       {
 91     *m_iterator++ = t1;
 92     *m_iterator++ = t2;
 93       }
 94       else
 95     if(CGAL::do_intersect(t2,s1))
 96     {
 97       *m_iterator++ = t1;
 98       *m_iterator++ = t2;
 99     }
100     return;
101     }
102 
103     // check for geometric intersection
104     Triangle t1( h->vertex()->point(),
105       h->next()->vertex()->point(),
106       h->next()->next()->vertex()->point());
107     Triangle t2( g->vertex()->point(),
108       g->next()->vertex()->point(),
109       g->next()->next()->vertex()->point());
110     if(CGAL::do_intersect(t1, t2))
111     {
112       *m_iterator++ = t1;
113       *m_iterator++ = t2;
114     }
115   } // end operator ()
116 }; // end struct Intersect_facets
117 
118 template <class Polyhedron, class Kernel, class OutputIterator>
119 void self_intersect(const Polyhedron& polyhedron,
120             OutputIterator out)
121 {
122   typedef CGAL::Bbox_3                  Bbox; // always double
123   typedef typename Kernel::Point_3      Point;
124   typedef typename Kernel::Vector_3     Vector;
125   typedef typename Kernel::Triangle_3   Triangle;
126   typedef typename Kernel::Segment_3    Segment;
127   typedef typename Polyhedron::Halfedge_const_handle    Halfedge_const_handle;
128   typedef typename Polyhedron::Facet_const_iterator     Facet_const_iterator;
129   typedef typename Polyhedron::Facet_const_handle       Facet_const_handle;
130   typedef typename CGAL::Box_intersection_d::Box_with_handle_d<double, 3, Facet_const_handle> Box;
131 
132   // make one box per facet
133   std::vector<Box> boxes;
134   boxes.reserve(polyhedron.size_of_facets());
135 
136   Facet_const_iterator f;
137   for(f = polyhedron.facets_begin();
138     f != polyhedron.facets_end();
139     f++)
140     boxes.push_back(Box( f->halfedge()->vertex()->point().bbox() +
141     f->halfedge()->next()->vertex()->point().bbox() +
142     f->halfedge()->next()->next()->vertex()->point().bbox(),
143     f));
144 
145   // generate box pointers
146   std::vector<const Box*> box_ptr;
147   box_ptr.reserve(polyhedron.size_of_facets());
148   typename std::vector<Box>::iterator b;
149   for(b = boxes.begin();
150     b != boxes.end();
151     b++)
152     box_ptr.push_back(&*b);
153 
154   // compute self-intersections filtered out by boxes
155   Intersect_facets<Polyhedron,Kernel,OutputIterator> intersect_facets(out);
156   std::ptrdiff_t cutoff = 2000;
157   CGAL::box_self_intersection_d(box_ptr.begin(), box_ptr.end(),intersect_facets,cutoff);
158 
159 } // end self_intersect
160 
161 #endif // _SELF_INTERSECT_H_