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_