Bundle updated version of carve. Should be no functional changes, small code cleanup
[blender-staging.git] / extern / carve / lib / intersect.cpp
1 // Begin License:
2 // Copyright (C) 2006-2011 Tobias Sargeant (tobias.sargeant@gmail.com).
3 // All rights reserved.
4 //
5 // This file is part of the Carve CSG Library (http://carve-csg.com/)
6 //
7 // This file may be used under the terms of the GNU General Public
8 // License version 2.0 as published by the Free Software Foundation
9 // and appearing in the file LICENSE.GPL2 included in the packaging of
10 // this file.
11 //
12 // This file is provided "AS IS" with NO WARRANTY OF ANY KIND,
13 // INCLUDING THE WARRANTIES OF DESIGN, MERCHANTABILITY AND FITNESS FOR
14 // A PARTICULAR PURPOSE.
15 // End:
16
17
18 #if defined(HAVE_CONFIG_H)
19 #  include <carve_config.h>
20 #endif
21
22 #include <carve/csg.hpp>
23 #include <carve/pointset.hpp>
24 #include <carve/polyline.hpp>
25
26 #include <list>
27 #include <set>
28 #include <iostream>
29
30 #include <algorithm>
31
32 #include "csg_detail.hpp"
33 #include "csg_data.hpp"
34
35 #include "intersect_debug.hpp"
36 #include "intersect_common.hpp"
37 #include "intersect_classify_common.hpp"
38
39 #include "csg_collector.hpp"
40
41 #include <carve/timing.hpp>
42 #include <carve/colour.hpp>
43
44
45
46
47 carve::csg::VertexPool::VertexPool() {
48 }
49
50 carve::csg::VertexPool::~VertexPool() {
51 }
52
53 void carve::csg::VertexPool::reset() {
54   pool.clear();
55 }
56
57 carve::csg::VertexPool::vertex_t *carve::csg::VertexPool::get(const vertex_t::vector_t &v) {
58   if (!pool.size() || pool.back().size() == blocksize) {
59     pool.push_back(std::vector<vertex_t>());
60     pool.back().reserve(blocksize);
61   }
62   pool.back().push_back(vertex_t(v));
63   return &pool.back().back();
64 }
65
66 bool carve::csg::VertexPool::inPool(vertex_t *v) const {
67   for (pool_t::const_iterator i = pool.begin(); i != pool.end(); ++i) {
68     if (v >= &(i->front()) && v <= &(i->back())) return true;
69   }
70   return false;
71 }
72
73
74
75 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
76 void writePLY(const std::string &out_file, const carve::point::PointSet *points, bool ascii);
77 void writePLY(const std::string &out_file, const carve::line::PolylineSet *lines, bool ascii);
78 void writePLY(const std::string &out_file, const carve::mesh::MeshSet<3> *poly, bool ascii);
79
80 static carve::mesh::MeshSet<3> *faceLoopsToPolyhedron(const carve::csg::FaceLoopList &fl) {
81   std::vector<carve::mesh::MeshSet<3>::face_t *> faces;
82   faces.reserve(fl.size());
83   for (carve::csg::FaceLoop *f = fl.head; f; f = f->next) {
84     faces.push_back(f->orig_face->create(f->vertices.begin(), f->vertices.end(), false));
85   }
86   carve::mesh::MeshSet<3> *poly = new carve::mesh::MeshSet<3>(faces);
87
88   return poly;
89 }
90 #endif
91
92 namespace {
93   /** 
94    * \brief Sort a range [\a beg, \a end) of vertices in order of increasing dot product of vertex - \a base on \dir.
95    * 
96    * @tparam[in] T a forward iterator type.
97    * @param[in] dir The direction in which to sort vertices.
98    * @param[in] base 
99    * @param[in] beg The start of the vertex range to sort.
100    * @param[in] end The end of the vertex range to sort.
101    * @param[out] out The sorted vertex result.
102    * @param[in] size_hint A hint regarding the size of the output
103    *            vector (to avoid needing to be able to calculate \a
104    *            end - \a beg).
105    */
106   template<typename iter_t>
107   void orderVertices(iter_t beg, const iter_t end,
108                      const carve::mesh::MeshSet<3>::vertex_t::vector_t &dir,
109                      const carve::mesh::MeshSet<3>::vertex_t::vector_t &base,
110                      std::vector<carve::mesh::MeshSet<3>::vertex_t *> &out) {
111     typedef std::vector<std::pair<double, carve::mesh::MeshSet<3>::vertex_t *> > DVVector;
112     std::vector<std::pair<double, carve::mesh::MeshSet<3>::vertex_t *> > ordered_vertices;
113
114     ordered_vertices.reserve(std::distance(beg, end));
115   
116     for (; beg != end; ++beg) {
117       carve::mesh::MeshSet<3>::vertex_t *v = (*beg);
118       ordered_vertices.push_back(std::make_pair(carve::geom::dot(v->v - base, dir), v));
119     }
120   
121     std::sort(ordered_vertices.begin(), ordered_vertices.end());
122   
123     out.clear();
124     out.reserve(ordered_vertices.size());
125     for (DVVector::const_iterator
126            i = ordered_vertices.begin(), e = ordered_vertices.end();
127          i != e;
128          ++i) {
129       out.push_back((*i).second);
130     }
131   }
132
133
134
135   /** 
136    * 
137    * 
138    * @param dir 
139    * @param base 
140    * @param beg 
141    * @param end 
142    */
143   template<typename iter_t>
144   void selectOrderingProjection(iter_t beg, const iter_t end,
145                                 carve::mesh::MeshSet<3>::vertex_t::vector_t &dir,
146                                 carve::mesh::MeshSet<3>::vertex_t::vector_t &base) {
147     double dx, dy, dz;
148     carve::mesh::MeshSet<3>::vertex_t *min_x, *min_y, *min_z, *max_x, *max_y, *max_z;
149     if (beg == end) return;
150     min_x = max_x = min_y = max_y = min_z = max_z = *beg++;
151     for (; beg != end; ++beg) {
152       if (min_x->v.x > (*beg)->v.x) min_x = *beg;
153       if (min_y->v.y > (*beg)->v.y) min_y = *beg;
154       if (min_z->v.z > (*beg)->v.z) min_z = *beg;
155       if (max_x->v.x < (*beg)->v.x) max_x = *beg;
156       if (max_y->v.y < (*beg)->v.y) max_y = *beg;
157       if (max_z->v.z < (*beg)->v.z) max_z = *beg;
158     }
159
160     dx = max_x->v.x - min_x->v.x;
161     dy = max_y->v.y - min_y->v.y;
162     dz = max_z->v.z - min_z->v.z;
163
164     if (dx > dy) {
165       if (dx > dz) {
166         dir = max_x->v - min_x->v; base = min_x->v;
167       } else {
168         dir = max_z->v - min_z->v; base = min_z->v;
169       }
170     } else {
171       if (dy > dz) {
172         dir = max_y->v - min_y->v; base = min_y->v;
173       } else {
174         dir = max_z->v - min_z->v; base = min_z->v;
175       }
176     }
177   }
178 }
179
180 namespace {
181   struct dump_data {
182     carve::mesh::MeshSet<3>::vertex_t *i_pt;
183     carve::csg::IObj i_src;
184     carve::csg::IObj i_tgt;
185     dump_data(carve::mesh::MeshSet<3>::vertex_t *_i_pt,
186               carve::csg::IObj _i_src,
187               carve::csg::IObj _i_tgt) : i_pt(_i_pt), i_src(_i_src), i_tgt(_i_tgt) {
188     }
189   };
190
191
192
193   struct dump_sort {
194     bool operator()(const dump_data &a, const dump_data &b) const {
195       if (a.i_pt->v.x < b.i_pt->v.x) return true;
196       if (a.i_pt->v.x > b.i_pt->v.x) return false;
197       if (a.i_pt->v.y < b.i_pt->v.y) return true;
198       if (a.i_pt->v.y > b.i_pt->v.y) return false;
199       if (a.i_pt->v.z < b.i_pt->v.z) return true;
200       if (a.i_pt->v.z > b.i_pt->v.z) return false;
201       return false;
202     }
203   };
204
205
206
207   void dump_intersections(std::ostream &out, carve::csg::Intersections &csg_intersections) {
208     std::vector<dump_data> temp;
209
210     for (carve::csg::Intersections::const_iterator
211         i = csg_intersections.begin(),
212         ie = csg_intersections.end();
213         i != ie;
214         ++i) {
215       const carve::csg::IObj &i_src = ((*i).first);
216
217       for (carve::csg::Intersections::mapped_type::const_iterator
218           j = (*i).second.begin(),
219           je = (*i).second.end();
220           j != je;
221           ++j) {
222         const carve::csg::IObj &i_tgt = ((*j).first);
223         carve::mesh::MeshSet<3>::vertex_t *i_pt = ((*j).second);
224         temp.push_back(dump_data(i_pt, i_src, i_tgt));
225       }
226     }
227
228     std::sort(temp.begin(), temp.end(), dump_sort());
229
230     for (size_t i = 0; i < temp.size(); ++i) {
231       const carve::csg::IObj &i_src = temp[i].i_src;
232       const carve::csg::IObj &i_tgt = temp[i].i_tgt;
233       out
234         << "INTERSECTION: " << temp[i].i_pt << " (" << temp[i].i_pt->v << ") "
235         << "is " << i_src << ".." << i_tgt << std::endl;
236     }
237
238 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
239     std::vector<carve::geom3d::Vector> vertices;
240
241     for (carve::csg::Intersections::const_iterator
242         i = csg_intersections.begin(),
243         ie = csg_intersections.end();
244         i != ie;
245         ++i) {
246       for (carve::csg::Intersections::mapped_type::const_iterator
247           j = (*i).second.begin(),
248           je = (*i).second.end();
249           j != je;
250           ++j) {
251         carve::mesh::MeshSet<3>::vertex_t *i_pt = ((*j).second);
252         vertices.push_back(i_pt->v);
253       }
254     }
255
256     carve::point::PointSet points(vertices);
257
258     std::string outf("/tmp/intersection-points.ply");
259     ::writePLY(outf, &points, true);
260 #endif
261   }
262
263
264
265   /** 
266    * \brief Populate a collection with the faces adjoining an edge.
267    * 
268    * @tparam face_set_t A collection type.
269    * @param e The edge for which to collect adjoining faces.
270    * @param faces 
271    */
272   template<typename face_set_t>
273   inline void facesForVertex(carve::mesh::MeshSet<3>::vertex_t *v,
274                              const carve::csg::detail::VEVecMap &ve,
275                              face_set_t &faces) {
276     carve::csg::detail::VEVecMap::const_iterator vi = ve.find(v);
277     if (vi != ve.end()) {
278       for (carve::csg::detail::VEVecMap::data_type::const_iterator i = (*vi).second.begin(); i != (*vi).second.end(); ++i) {
279         faces.insert((*i)->face);
280       }
281     }
282   }
283
284   /** 
285    * \brief Populate a collection with the faces adjoining an edge.
286    * 
287    * @tparam face_set_t A collection type.
288    * @param e The edge for which to collect adjoining faces.
289    * @param faces 
290    */
291   template<typename face_set_t>
292   inline void facesForEdge(carve::mesh::MeshSet<3>::edge_t *e,
293                            face_set_t &faces) {
294     faces.insert(e->face);
295   }
296
297   /** 
298    * \brief Populate a collection with the faces adjoining a face.
299    * 
300    * @tparam face_set_t A collection type.
301    * @param f The face for which to collect adjoining faces.
302    * @param faces 
303    */
304   template<typename face_set_t>
305   inline void facesForFace(carve::mesh::MeshSet<3>::face_t *f,
306                            face_set_t &faces) {
307     faces.insert(f);
308   }
309
310   /** 
311    * \brief Populate a collection with the faces adjoining an intersection object.
312    * 
313    * @tparam face_set_t A collection type holding const carve::poly::Polyhedron::face_t *.
314    * @param obj The intersection object for which to collect adjoining faces.
315    * @param faces 
316    */
317   template<typename face_set_t>
318   void facesForObject(const carve::csg::IObj &obj,
319                       const carve::csg::detail::VEVecMap &ve,
320                       face_set_t &faces) {
321     switch (obj.obtype) {
322     case carve::csg::IObj::OBTYPE_VERTEX:
323       facesForVertex(obj.vertex, ve, faces);
324       break;
325
326     case carve::csg::IObj::OBTYPE_EDGE:
327       facesForEdge(obj.edge, faces);
328       break;
329
330     case  carve::csg::IObj::OBTYPE_FACE:
331       facesForFace(obj.face, faces);
332       break;
333
334     default:
335       break;
336     }
337   }
338
339
340
341 }
342
343
344
345 bool carve::csg::CSG::Hooks::hasHook(unsigned hook_num) {
346   return hooks[hook_num].size() > 0;
347 }
348
349 void carve::csg::CSG::Hooks::intersectionVertex(const carve::mesh::MeshSet<3>::vertex_t *vertex,
350                                                 const IObjPairSet &intersections) {
351   for (std::list<Hook *>::iterator j = hooks[INTERSECTION_VERTEX_HOOK].begin();
352        j != hooks[INTERSECTION_VERTEX_HOOK].end();
353        ++j) {
354     (*j)->intersectionVertex(vertex, intersections);
355   }
356 }
357
358 void carve::csg::CSG::Hooks::processOutputFace(std::vector<carve::mesh::MeshSet<3>::face_t *> &faces,
359                                                const carve::mesh::MeshSet<3>::face_t *orig_face,
360                                                bool flipped) {
361   for (std::list<Hook *>::iterator j = hooks[PROCESS_OUTPUT_FACE_HOOK].begin();
362        j != hooks[PROCESS_OUTPUT_FACE_HOOK].end();
363        ++j) {
364     (*j)->processOutputFace(faces, orig_face, flipped);
365   }
366 }
367
368 void carve::csg::CSG::Hooks::resultFace(const carve::mesh::MeshSet<3>::face_t *new_face,
369                                         const carve::mesh::MeshSet<3>::face_t *orig_face,
370                                         bool flipped) {
371   for (std::list<Hook *>::iterator j = hooks[RESULT_FACE_HOOK].begin();
372        j != hooks[RESULT_FACE_HOOK].end();
373        ++j) {
374     (*j)->resultFace(new_face, orig_face, flipped);
375   }
376 }
377
378 void carve::csg::CSG::Hooks::registerHook(Hook *hook, unsigned hook_bits) {
379   for (unsigned i = 0; i < HOOK_MAX; ++i) {
380     if (hook_bits & (1U << i)) {
381       hooks[i].push_back(hook);
382     }
383   }
384 }
385
386 void carve::csg::CSG::Hooks::unregisterHook(Hook *hook) {
387   for (unsigned i = 0; i < HOOK_MAX; ++i) {
388     hooks[i].erase(std::remove(hooks[i].begin(), hooks[i].end(), hook), hooks[i].end());
389   }
390 }
391
392 void carve::csg::CSG::Hooks::reset() {
393   for (unsigned i = 0; i < HOOK_MAX; ++i) {
394     for (std::list<Hook *>::iterator j = hooks[i].begin(); j != hooks[i].end(); ++j) {
395       delete (*j);
396     }
397     hooks[i].clear();
398   }
399 }
400
401 carve::csg::CSG::Hooks::Hooks() : hooks() {
402   hooks.resize(HOOK_MAX);
403 }
404  
405 carve::csg::CSG::Hooks::~Hooks() {
406   reset();
407 }
408
409
410
411 void carve::csg::CSG::makeVertexIntersections() {
412   static carve::TimingName FUNC_NAME("CSG::makeVertexIntersections()");
413   carve::TimingBlock block(FUNC_NAME);
414   vertex_intersections.clear();
415   for (Intersections::const_iterator
416          i = intersections.begin(),
417          ie = intersections.end();
418        i != ie;
419        ++i) {
420     const IObj &i_src = ((*i).first);
421
422     for (Intersections::mapped_type::const_iterator
423            j = (*i).second.begin(),
424            je = (*i).second.end();
425          j != je;
426          ++j) {
427       const IObj &i_tgt = ((*j).first);
428       carve::mesh::MeshSet<3>::vertex_t *i_pt = ((*j).second);
429
430       vertex_intersections[i_pt].insert(std::make_pair(i_src, i_tgt));
431     }
432   }
433 }
434
435
436
437 static carve::mesh::MeshSet<3>::vertex_t *chooseWeldPoint(
438     const carve::csg::detail::VSet &equivalent,
439     carve::csg::VertexPool &vertex_pool) {
440   // XXX: choose a better weld point.
441   if (!equivalent.size()) return NULL;
442
443   for (carve::csg::detail::VSet::const_iterator
444          i = equivalent.begin(), e = equivalent.end();
445        i != e;
446        ++i) {
447     if (!vertex_pool.inPool((*i))) return (*i);
448   }
449   return *equivalent.begin();
450 }
451
452
453
454 static const carve::mesh::MeshSet<3>::vertex_t *weld(
455     const carve::csg::detail::VSet &equivalent,
456     carve::csg::VertexIntersections &vertex_intersections,
457     carve::csg::VertexPool &vertex_pool) {
458   carve::mesh::MeshSet<3>::vertex_t *weld_point = chooseWeldPoint(equivalent, vertex_pool);
459
460 #if defined(CARVE_DEBUG)
461   std::cerr << "weld: " << equivalent.size() << " vertices ( ";
462   for (carve::csg::detail::VSet::const_iterator
463          i = equivalent.begin(), e = equivalent.end();
464        i != e;
465        ++i) {
466     const carve::mesh::MeshSet<3>::vertex_t *v = (*i);
467     std::cerr << " " << v;
468   }
469   std::cerr << ") to " << weld_point << std::endl;
470 #endif
471
472   if (!weld_point) return NULL;
473
474   carve::csg::VertexIntersections::mapped_type &weld_tgt = (vertex_intersections[weld_point]);
475
476   for (carve::csg::detail::VSet::const_iterator
477          i = equivalent.begin(), e = equivalent.end();
478        i != e;
479        ++i) {
480     carve::mesh::MeshSet<3>::vertex_t *v = (*i);
481
482     if (v != weld_point) {
483       carve::csg::VertexIntersections::iterator j = vertex_intersections.find(v);
484
485       if (j != vertex_intersections.end()) {
486         weld_tgt.insert((*j).second.begin(), (*j).second.end());
487         vertex_intersections.erase(j);
488       }
489     }
490   }
491   return weld_point;
492 }
493
494
495
496 void carve::csg::CSG::groupIntersections() {
497 #if 0 // old code, to be removed.
498   static carve::TimingName GROUP_INTERSECTONS("groupIntersections()");
499
500   carve::TimingBlock block(GROUP_INTERSECTONS);
501   
502   std::vector<carve::mesh::MeshSet<3>::vertex_t *> vertices;
503   detail::VVSMap graph;
504 #if defined(CARVE_DEBUG)
505   std::cerr << "groupIntersections()" << ": vertex_intersections.size()==" << vertex_intersections.size() << std::endl;
506 #endif
507
508   vertices.reserve(vertex_intersections.size());
509   for (carve::csg::VertexIntersections::const_iterator
510          i = vertex_intersections.begin(),
511          e = vertex_intersections.end();
512        i != e;
513        ++i) 
514     {
515       vertices.push_back((*i).first);
516     }
517   carve::geom3d::AABB aabb;
518   aabb.fit(vertices.begin(), vertices.end(), carve::poly::vec_adapt_vertex_ptr());
519   Octree vertex_intersections_octree;
520   vertex_intersections_octree.setBounds(aabb);
521
522   vertex_intersections_octree.addVertices(vertices);
523       
524   std::vector<carve::mesh::MeshSet<3>::vertex_t *> out;
525   for (size_t i = 0, l = vertices.size(); i != l; ++i) {
526     // let's find all the vertices near this one. 
527     out.clear();
528     vertex_intersections_octree.findVerticesNearAllowDupes(vertices[i]->v, out);
529
530     for (size_t j = 0; j < out.size(); ++j) {
531       if (vertices[i] != out[j] && carve::geom::equal(vertices[i]->v, out[j]->v)) {
532 #if defined(CARVE_DEBUG)
533         std::cerr << "EQ: " << vertices[i] << "," << out[j] << " " << vertices[i]->v << "," << out[j]->v << std::endl;
534 #endif
535         graph[vertices[i]].insert(out[j]);
536         graph[out[j]].insert(vertices[i]);
537       }
538     }
539   }
540
541   detail::VSet visited, open;
542   while (graph.size()) {
543     visited.clear();
544     open.clear();
545     detail::VVSMap::iterator i = graph.begin();
546     open.insert((*i).first);
547     while (open.size()) {
548       detail::VSet::iterator t = open.begin();
549       const carve::mesh::MeshSet<3>::vertex_t *o = (*t);
550       open.erase(t);
551       i = graph.find(o);
552       CARVE_ASSERT(i != graph.end());
553       visited.insert(o);
554       for (detail::VVSMap::mapped_type::const_iterator
555              j = (*i).second.begin(),
556              je = (*i).second.end();
557            j != je;
558            ++j) {
559         if (visited.count((*j)) == 0) {
560           open.insert((*j));
561         }
562       }
563       graph.erase(i);
564     }
565     weld(visited, vertex_intersections, vertex_pool);
566   }
567 #endif
568 }
569
570
571
572 void carve::csg::CSG::intersectingFacePairs(detail::Data &data) {
573   static carve::TimingName FUNC_NAME("CSG::intersectingFacePairs()");
574   carve::TimingBlock block(FUNC_NAME);
575
576   // iterate over all intersection points.
577   for (VertexIntersections::const_iterator i = vertex_intersections.begin(), ie = vertex_intersections.end(); i != ie; ++i) {
578     carve::mesh::MeshSet<3>::vertex_t *i_pt = ((*i).first);
579     detail::VFSMap::mapped_type &face_set = (data.fmap_rev[i_pt]);
580
581     // for all pairs of intersecting objects at this point
582     for (VertexIntersections::data_type::const_iterator j = (*i).second.begin(), je = (*i).second.end(); j != je; ++j) {
583       const IObj &i_src = ((*j).first);
584       const IObj &i_tgt = ((*j).second);
585
586       // work out the faces involved. this updates fmap_rev.
587       facesForObject(i_src, data.vert_to_edges, face_set);
588       facesForObject(i_tgt, data.vert_to_edges, face_set);
589
590       // record the intersection with respect to any involved vertex.
591       if (i_src.obtype == IObj::OBTYPE_VERTEX) data.vmap[i_src.vertex] = i_pt;
592       if (i_tgt.obtype == IObj::OBTYPE_VERTEX) data.vmap[i_tgt.vertex] = i_pt;
593
594       // record the intersection with respect to any involved edge.
595       if (i_src.obtype == IObj::OBTYPE_EDGE) data.emap[i_src.edge].insert(i_pt);
596       if (i_tgt.obtype == IObj::OBTYPE_EDGE) data.emap[i_tgt.edge].insert(i_pt);
597     }
598
599     // record the intersection with respect to each face.
600     for (carve::csg::detail::VFSMap::mapped_type::const_iterator k = face_set.begin(), ke = face_set.end(); k != ke; ++k) {
601       carve::mesh::MeshSet<3>::face_t *f = (*k);
602       data.fmap[f].insert(i_pt);
603     }
604   }
605 }
606
607
608
609 void carve::csg::CSG::_generateVertexVertexIntersections(carve::mesh::MeshSet<3>::vertex_t *va,
610                                                          carve::mesh::MeshSet<3>::edge_t *eb) {
611   if (intersections.intersects(va, eb->v1())) {
612     return;
613   }
614
615   double d_v1 = carve::geom::distance2(va->v, eb->v1()->v);
616
617   if  (d_v1 < carve::EPSILON2) {
618     intersections.record(va, eb->v1(), va);
619   }
620 }
621
622
623
624 void carve::csg::CSG::generateVertexVertexIntersections(carve::mesh::MeshSet<3>::face_t *a,
625                                                         const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
626   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
627
628   ea = a->edge;
629   do {
630     for (size_t i = 0; i < b.size(); ++i) {
631       carve::mesh::MeshSet<3>::face_t *t = b[i];
632       eb = t->edge;
633       do {
634         _generateVertexVertexIntersections(ea->v1(), eb);
635         eb = eb->next;
636       } while (eb != t->edge);
637     }
638     ea = ea->next;
639   } while (ea != a->edge);
640 }
641
642
643
644 void carve::csg::CSG::_generateVertexEdgeIntersections(carve::mesh::MeshSet<3>::vertex_t *va,
645                                                        carve::mesh::MeshSet<3>::edge_t *eb) {
646   if (intersections.intersects(va, eb)) {
647     return;
648   }
649
650   carve::geom::aabb<3> eb_aabb;
651   eb_aabb.fit(eb->v1()->v, eb->v2()->v);
652   if (eb_aabb.maxAxisSeparation(va->v) > carve::EPSILON) {
653     return;
654   }
655
656   double a = cross(eb->v2()->v - eb->v1()->v, va->v - eb->v1()->v).length2();
657   double b = (eb->v2()->v - eb->v1()->v).length2();
658
659   if (a < b * carve::EPSILON2) {
660     // vertex-edge intersection
661     intersections.record(eb, va, va);
662     if (eb->rev) intersections.record(eb->rev, va, va);
663   }
664 }
665
666
667
668 void carve::csg::CSG::generateVertexEdgeIntersections(carve::mesh::MeshSet<3>::face_t *a,
669                                                       const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
670   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
671
672   ea = a->edge;
673   do {
674     for (size_t i = 0; i < b.size(); ++i) {
675       carve::mesh::MeshSet<3>::face_t *t = b[i];
676       eb = t->edge;
677       do {
678         _generateVertexEdgeIntersections(ea->v1(), eb);
679         eb = eb->next;
680       } while (eb != t->edge);
681     }
682     ea = ea->next;
683   } while (ea != a->edge);
684 }
685
686
687
688 void carve::csg::CSG::_generateEdgeEdgeIntersections(carve::mesh::MeshSet<3>::edge_t *ea,
689                                                      carve::mesh::MeshSet<3>::edge_t *eb) {
690   if (intersections.intersects(ea, eb)) {
691     return;
692   }
693
694   carve::mesh::MeshSet<3>::vertex_t *v1 = ea->v1(), *v2 = ea->v2();
695   carve::mesh::MeshSet<3>::vertex_t *v3 = eb->v1(), *v4 = eb->v2();
696
697   carve::geom::aabb<3> ea_aabb, eb_aabb;
698   ea_aabb.fit(v1->v, v2->v);
699   eb_aabb.fit(v3->v, v4->v);
700   if (ea_aabb.maxAxisSeparation(eb_aabb) > EPSILON) return;
701
702   carve::mesh::MeshSet<3>::vertex_t::vector_t p1, p2;
703   double mu1, mu2;
704
705   switch (carve::geom3d::rayRayIntersection(carve::geom3d::Ray(v2->v - v1->v, v1->v),
706                                             carve::geom3d::Ray(v4->v - v3->v, v3->v),
707                                             p1, p2, mu1, mu2)) {
708   case carve::RR_INTERSECTION: {
709     // edges intersect
710     if (mu1 >= 0.0 && mu1 <= 1.0 && mu2 >= 0.0 && mu2 <= 1.0) {
711       carve::mesh::MeshSet<3>::vertex_t *p = vertex_pool.get((p1 + p2) / 2.0);
712       intersections.record(ea, eb, p);
713       if (ea->rev) intersections.record(ea->rev, eb, p);
714       if (eb->rev) intersections.record(ea, eb->rev, p);
715       if (ea->rev && eb->rev) intersections.record(ea->rev, eb->rev, p);
716     }
717     break;
718   }
719   case carve::RR_PARALLEL: {
720     // edges parallel. any intersection of this type should have
721     // been handled by generateVertexEdgeIntersections().
722     break;
723   }
724   case carve::RR_DEGENERATE: {
725     throw carve::exception("degenerate edge");
726     break;
727   }
728   case carve::RR_NO_INTERSECTION: {
729     break;
730   }
731   }
732 }
733
734
735
736 void carve::csg::CSG::generateEdgeEdgeIntersections(carve::mesh::MeshSet<3>::face_t *a,
737                                                     const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
738   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
739
740   ea = a->edge;
741   do {
742     for (size_t i = 0; i < b.size(); ++i) {
743       carve::mesh::MeshSet<3>::face_t *t = b[i];
744       eb = t->edge;
745       do {
746         _generateEdgeEdgeIntersections(ea, eb);
747         eb = eb->next;
748       } while (eb != t->edge);
749     }
750     ea = ea->next;
751   } while (ea != a->edge);
752 }
753
754
755
756 void carve::csg::CSG::_generateVertexFaceIntersections(carve::mesh::MeshSet<3>::face_t *fa,
757                                                        carve::mesh::MeshSet<3>::edge_t *eb) {
758   if (intersections.intersects(eb->v1(), fa)) {
759     return;
760   }
761
762   double d1 = carve::geom::distance(fa->plane, eb->v1()->v);
763
764   if (fabs(d1) < carve::EPSILON &&
765       fa->containsPoint(eb->v1()->v)) {
766     intersections.record(eb->v1(), fa, eb->v1());
767   }
768 }
769
770
771
772 void carve::csg::CSG::generateVertexFaceIntersections(carve::mesh::MeshSet<3>::face_t *a,
773                                                       const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
774   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
775
776   for (size_t i = 0; i < b.size(); ++i) {
777     carve::mesh::MeshSet<3>::face_t *t = b[i];
778     eb = t->edge;
779     do {
780       _generateVertexFaceIntersections(a, eb);
781       eb = eb->next;
782     } while (eb != t->edge);
783   }
784 }
785
786
787
788 void carve::csg::CSG::_generateEdgeFaceIntersections(carve::mesh::MeshSet<3>::face_t *fa,
789                                                      carve::mesh::MeshSet<3>::edge_t *eb) {
790   if (intersections.intersects(eb, fa)) {
791     return;
792   }
793
794   carve::mesh::MeshSet<3>::vertex_t::vector_t _p;
795   if (fa->simpleLineSegmentIntersection(carve::geom3d::LineSegment(eb->v1()->v, eb->v2()->v), _p)) {
796     carve::mesh::MeshSet<3>::vertex_t *p = vertex_pool.get(_p);
797     intersections.record(eb, fa, p);
798     if (eb->rev) intersections.record(eb->rev, fa, p);
799   }
800 }
801
802
803
804 void carve::csg::CSG::generateEdgeFaceIntersections(carve::mesh::MeshSet<3>::face_t *a,
805                                                     const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
806   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
807
808   for (size_t i = 0; i < b.size(); ++i) {
809     carve::mesh::MeshSet<3>::face_t *t = b[i];
810     eb = t->edge;
811     do {
812       _generateEdgeFaceIntersections(a, eb);
813       eb = eb->next;
814     } while (eb != t->edge);
815   }
816 }
817
818
819
820 void carve::csg::CSG::generateIntersectionCandidates(carve::mesh::MeshSet<3> *a,
821                                                      const face_rtree_t *a_node,
822                                                      carve::mesh::MeshSet<3> *b,
823                                                      const face_rtree_t *b_node,
824                                                      face_pairs_t &face_pairs,
825                                                      bool descend_a) {
826   if (!a_node->bbox.intersects(b_node->bbox)) {
827     return;
828   }
829
830   if (a_node->child && (descend_a || !b_node->child)) {
831     for (face_rtree_t *node = a_node->child; node; node = node->sibling) {
832       generateIntersectionCandidates(a, node, b, b_node, face_pairs, false);
833     }
834   } else if (b_node->child) {
835     for (face_rtree_t *node = b_node->child; node; node = node->sibling) {
836       generateIntersectionCandidates(a, a_node, b, node, face_pairs, true);
837     }
838   } else {
839     for (size_t i = 0; i < a_node->data.size(); ++i) {
840       carve::mesh::MeshSet<3>::face_t *fa = a_node->data[i];
841       carve::geom::aabb<3> aabb_a = fa->getAABB();
842       if (aabb_a.maxAxisSeparation(b_node->bbox) > carve::EPSILON) continue;
843
844       for (size_t j = 0; j < b_node->data.size(); ++j) {
845         carve::mesh::MeshSet<3>::face_t *fb = b_node->data[j];
846         carve::geom::aabb<3> aabb_b = fb->getAABB();
847         if (aabb_b.maxAxisSeparation(aabb_a) > carve::EPSILON) continue;
848
849         std::pair<double, double> a_ra = fa->rangeInDirection(fa->plane.N, fa->edge->vert->v);
850         std::pair<double, double> b_ra = fb->rangeInDirection(fa->plane.N, fa->edge->vert->v);
851         if (carve::rangeSeparation(a_ra, b_ra) > carve::EPSILON) continue;
852
853         std::pair<double, double> a_rb = fa->rangeInDirection(fb->plane.N, fb->edge->vert->v);
854         std::pair<double, double> b_rb = fb->rangeInDirection(fb->plane.N, fb->edge->vert->v);
855         if (carve::rangeSeparation(a_rb, b_rb) > carve::EPSILON) continue;
856
857         if (!facesAreCoplanar(fa, fb)) {
858           face_pairs[fa].push_back(fb); 
859           face_pairs[fb].push_back(fa);
860         }
861      }
862     }
863   }
864 }
865
866
867
868
869 void carve::csg::CSG::generateIntersections(carve::mesh::MeshSet<3> *a,
870                                             const face_rtree_t *a_rtree,
871                                             carve::mesh::MeshSet<3> *b,
872                                             const face_rtree_t *b_rtree,
873                                             detail::Data &data) {
874   face_pairs_t face_pairs;
875   generateIntersectionCandidates(a, a_rtree, b, b_rtree, face_pairs);
876
877   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
878     carve::mesh::MeshSet<3>::face_t *f = (*i).first;
879     carve::mesh::MeshSet<3>::edge_t *e = f->edge;
880     do {
881       data.vert_to_edges[e->v1()].push_back(e);
882       e = e->next;
883     } while (e != f->edge);
884   }
885
886   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
887     generateVertexVertexIntersections((*i).first, (*i).second);
888   }
889
890   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
891     generateVertexEdgeIntersections((*i).first, (*i).second);
892   }
893
894   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
895     generateEdgeEdgeIntersections((*i).first, (*i).second);
896   }
897
898   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
899     generateVertexFaceIntersections((*i).first, (*i).second);
900   }
901
902   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
903     generateEdgeFaceIntersections((*i).first, (*i).second);
904   }
905
906
907 #if defined(CARVE_DEBUG)
908   std::cerr << "makeVertexIntersections" << std::endl;
909 #endif
910   makeVertexIntersections();
911
912 #if defined(CARVE_DEBUG)
913   std::cerr << "  intersections.size() " << intersections.size() << std::endl;
914   map_histogram(std::cerr, intersections);
915   std::cerr << "  vertex_intersections.size() " << vertex_intersections.size() << std::endl;
916   map_histogram(std::cerr, vertex_intersections);
917 #endif
918
919 #if defined(CARVE_DEBUG) && defined(DEBUG_DRAW_INTERSECTIONS)
920   HOOK(drawIntersections(vertex_intersections););
921 #endif
922
923 #if defined(CARVE_DEBUG)
924   std::cerr << "  intersections.size() " << intersections.size() << std::endl;
925   std::cerr << "  vertex_intersections.size() " << vertex_intersections.size() << std::endl;
926 #endif
927
928   // notify about intersections.
929   if (hooks.hasHook(Hooks::INTERSECTION_VERTEX_HOOK)) {
930     for (VertexIntersections::const_iterator i = vertex_intersections.begin();
931          i != vertex_intersections.end();
932          ++i) {
933       hooks.intersectionVertex((*i).first, (*i).second);
934     }
935   }
936
937   // from here on, only vertex_intersections is used for intersection
938   // information.
939
940   // intersections still contains the vertex_to_face map. maybe that
941   // should be moved out into another class.
942   static_cast<Intersections::super>(intersections).clear();
943 }
944
945
946
947 carve::csg::CSG::CSG() {
948 }
949
950
951
952 /** 
953  * \brief For each intersected edge, decompose into a set of vertex pairs representing an ordered set of edge fragments.
954  * 
955  * @tparam[in,out] data Internal intersection data. data.emap is used to produce data.divided_edges.
956  */
957 void carve::csg::CSG::divideIntersectedEdges(detail::Data &data) {
958   static carve::TimingName FUNC_NAME("CSG::divideIntersectedEdges()");
959   carve::TimingBlock block(FUNC_NAME);
960
961   for (detail::EVSMap::const_iterator i = data.emap.begin(), ei = data.emap.end(); i != ei; ++i) {
962     carve::mesh::MeshSet<3>::edge_t *edge = (*i).first;
963     const detail::EVSMap::mapped_type &vertices = (*i).second;
964     std::vector<carve::mesh::MeshSet<3>::vertex_t *> &verts = data.divided_edges[edge];
965     orderVertices(vertices.begin(), vertices.end(),
966                   edge->v2()->v - edge->v1()->v, edge->v1()->v,
967                   verts);
968   }
969 }
970
971
972
973 carve::csg::CSG::~CSG() {
974 }
975
976
977
978 void carve::csg::CSG::makeFaceEdges(carve::csg::EdgeClassification &eclass,
979                                     detail::Data &data) {
980   detail::FSet face_b_set;
981   for (detail::FVSMap::const_iterator
982          i = data.fmap.begin(), ie = data.fmap.end();
983        i != ie;
984        ++i) {
985     carve::mesh::MeshSet<3>::face_t *face_a = (*i).first;
986     const detail::FVSMap::mapped_type &face_a_intersections = ((*i).second);
987     face_b_set.clear();
988
989     // work out the set of faces from the opposing polyhedron that intersect face_a.
990     for (detail::FVSMap::mapped_type::const_iterator
991            j = face_a_intersections.begin(), je = face_a_intersections.end();
992          j != je;
993          ++j) {
994       for (detail::VFSMap::mapped_type::const_iterator
995              k = data.fmap_rev[*j].begin(), ke = data.fmap_rev[*j].end();
996            k != ke;
997            ++k) {
998         carve::mesh::MeshSet<3>::face_t *face_b = (*k);
999         if (face_a != face_b && face_b->mesh->meshset != face_a->mesh->meshset) {
1000           face_b_set.insert(face_b);
1001         }
1002       }
1003     }
1004
1005     // run through each intersecting face.
1006     for (detail::FSet::const_iterator
1007            j = face_b_set.begin(), je = face_b_set.end();
1008          j != je;
1009          ++j) {
1010       carve::mesh::MeshSet<3>::face_t *face_b = (*j);
1011       const detail::FVSMap::mapped_type &face_b_intersections = (data.fmap[face_b]);
1012
1013       std::vector<carve::mesh::MeshSet<3>::vertex_t *> vertices;
1014       vertices.reserve(std::min(face_a_intersections.size(), face_b_intersections.size()));
1015
1016       // record the points of intersection between face_a and face_b
1017       std::set_intersection(face_a_intersections.begin(),
1018                             face_a_intersections.end(),
1019                             face_b_intersections.begin(),
1020                             face_b_intersections.end(),
1021                             std::back_inserter(vertices));
1022
1023 #if defined(CARVE_DEBUG)
1024       std::cerr << "face pair: "
1025                 << face_a << ":" << face_b
1026                 << " N(verts) " << vertices.size() << std::endl;
1027       for (std::vector<carve::mesh::MeshSet<3>::vertex_t *>::const_iterator i = vertices.begin(), e = vertices.end(); i != e; ++i) {
1028         std::cerr << (*i) << " " << (*i)->v << " ("
1029                   << carve::geom::distance(face_a->plane, (*i)->v) << ","
1030                   << carve::geom::distance(face_b->plane, (*i)->v) << ")"
1031                   << std::endl;
1032         //CARVE_ASSERT(carve::geom3d::distance(face_a->plane_eqn, *(*i)) < EPSILON);
1033         //CARVE_ASSERT(carve::geom3d::distance(face_b->plane_eqn, *(*i)) < EPSILON);
1034       }
1035 #endif
1036
1037       // if there are two points of intersection, then the added edge is simple to determine.
1038       if (vertices.size() == 2) {
1039         carve::mesh::MeshSet<3>::vertex_t *v1 = vertices[0];
1040         carve::mesh::MeshSet<3>::vertex_t *v2 = vertices[1];
1041         carve::geom3d::Vector c = (v1->v + v2->v) / 2;
1042
1043         // determine whether the midpoint of the implied edge is contained in face_a and face_b
1044
1045 #if defined(CARVE_DEBUG)
1046         std::cerr << "face_a->nVertices() = " << face_a->nVertices() << " face_a->containsPointInProjection(c) = " << face_a->containsPointInProjection(c) << std::endl;
1047         std::cerr << "face_b->nVertices() = " << face_b->nVertices() << " face_b->containsPointInProjection(c) = " << face_b->containsPointInProjection(c) << std::endl;
1048 #endif
1049
1050         if (face_a->containsPointInProjection(c) && face_b->containsPointInProjection(c)) {
1051 #if defined(CARVE_DEBUG)
1052           std::cerr << "adding edge: " << v1 << "-" << v2 << std::endl;
1053 #if defined(DEBUG_DRAW_FACE_EDGES)
1054           HOOK(drawEdge(v1, v2, 1, 1, 1, 1, 1, 1, 1, 1, 2.0););
1055 #endif
1056 #endif
1057           // record the edge, with class information.
1058           if (v1 > v2) std::swap(v1, v2);
1059           eclass[ordered_edge(v1, v2)] = carve::csg::EC2(carve::csg::EDGE_ON, carve::csg::EDGE_ON);
1060           data.face_split_edges[face_a].insert(std::make_pair(v1, v2));
1061           data.face_split_edges[face_b].insert(std::make_pair(v1, v2));
1062         }
1063         continue;
1064       }
1065
1066       // otherwise, it's more complex.
1067       carve::geom3d::Vector base, dir;
1068       std::vector<carve::mesh::MeshSet<3>::vertex_t *> ordered;
1069
1070       // skip coplanar edges. this simplifies the resulting
1071       // mesh. eventually all coplanar face regions of two polyhedra
1072       // must reach a point where they are no longer coplanar (or the
1073       // polyhedra are identical).
1074       if (!facesAreCoplanar(face_a, face_b)) {
1075         // order the intersection vertices (they must lie along a
1076         // vector, as the faces aren't coplanar).
1077         selectOrderingProjection(vertices.begin(), vertices.end(), dir, base);
1078         orderVertices(vertices.begin(), vertices.end(), dir, base, ordered);
1079
1080         // for each possible edge in the ordering, test the midpoint,
1081         // and record if it's contained in face_a and face_b.
1082         for (int k = 0, ke = (int)ordered.size() - 1; k < ke; ++k) {
1083           carve::mesh::MeshSet<3>::vertex_t *v1 = ordered[k];
1084           carve::mesh::MeshSet<3>::vertex_t *v2 = ordered[k + 1];
1085           carve::geom3d::Vector c = (v1->v + v2->v) / 2;
1086
1087 #if defined(CARVE_DEBUG)
1088           std::cerr << "testing edge: " << v1 << "-" << v2 << " at " << c << std::endl;
1089           std::cerr << "a: " << face_a->containsPointInProjection(c) << " b: " << face_b->containsPointInProjection(c) << std::endl;
1090           std::cerr << "face_a->containsPointInProjection(c): " << face_a->containsPointInProjection(c) << std::endl;
1091           std::cerr << "face_b->containsPointInProjection(c): " << face_b->containsPointInProjection(c) << std::endl;
1092 #endif
1093
1094           if (face_a->containsPointInProjection(c) && face_b->containsPointInProjection(c)) {
1095 #if defined(CARVE_DEBUG)
1096             std::cerr << "adding edge: " << v1 << "-" << v2 << std::endl;
1097 #if defined(DEBUG_DRAW_FACE_EDGES)
1098             HOOK(drawEdge(v1, v2, .5, .5, .5, 1, .5, .5, .5, 1, 2.0););
1099 #endif
1100 #endif
1101             // record the edge, with class information.
1102             if (v1 > v2) std::swap(v1, v2);
1103             eclass[ordered_edge(v1, v2)] = carve::csg::EC2(carve::csg::EDGE_ON, carve::csg::EDGE_ON);
1104             data.face_split_edges[face_a].insert(std::make_pair(v1, v2));
1105             data.face_split_edges[face_b].insert(std::make_pair(v1, v2));
1106           }
1107         }
1108       }
1109     }
1110   }
1111
1112
1113 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
1114   {
1115     V2Set edges;
1116     for (detail::FV2SMap::const_iterator i = data.face_split_edges.begin(); i != data.face_split_edges.end(); ++i) {
1117       edges.insert((*i).second.begin(), (*i).second.end());
1118     }
1119
1120     detail::VSet vertices;
1121     for (V2Set::const_iterator i = edges.begin(); i != edges.end(); ++i) {
1122       vertices.insert((*i).first);
1123       vertices.insert((*i).second);
1124     }
1125
1126     carve::line::PolylineSet intersection_graph;
1127     intersection_graph.vertices.resize(vertices.size());
1128     std::map<const carve::mesh::MeshSet<3>::vertex_t *, size_t> vmap;
1129
1130     size_t j = 0;
1131     for (detail::VSet::const_iterator i = vertices.begin(); i != vertices.end(); ++i) {
1132       intersection_graph.vertices[j].v = (*i)->v;
1133       vmap[(*i)] = j++;
1134     }
1135
1136     for (V2Set::const_iterator i = edges.begin(); i != edges.end(); ++i) {
1137       size_t line[2];
1138       line[0] = vmap[(*i).first];
1139       line[1] = vmap[(*i).second];
1140       intersection_graph.addPolyline(false, line, line + 2);
1141     }
1142
1143     std::string out("/tmp/intersection-edges.ply");
1144     ::writePLY(out, &intersection_graph, true);
1145   }
1146 #endif
1147 }
1148
1149
1150
1151 /** 
1152  * 
1153  * 
1154  * @param fll 
1155  */
1156 static void checkFaceLoopIntegrity(carve::csg::FaceLoopList &fll) {
1157   static carve::TimingName FUNC_NAME("CSG::checkFaceLoopIntegrity()");
1158   carve::TimingBlock block(FUNC_NAME);
1159
1160   std::unordered_map<carve::csg::V2, int> counts;
1161   for (carve::csg::FaceLoop *fl = fll.head; fl; fl = fl->next) {
1162     std::vector<carve::mesh::MeshSet<3>::vertex_t *> &loop = (fl->vertices);
1163     carve::mesh::MeshSet<3>::vertex_t *v1, *v2;
1164     v1 = loop[loop.size() - 1];
1165     for (unsigned i = 0; i < loop.size(); ++i) {
1166       v2 = loop[i];
1167       if (v1 < v2) {
1168         counts[std::make_pair(v1, v2)]++;
1169       } else {
1170         counts[std::make_pair(v2, v1)]--;
1171       }
1172       v1 = v2;
1173     }
1174   }
1175   for (std::unordered_map<carve::csg::V2, int>::const_iterator
1176          x = counts.begin(), xe = counts.end(); x != xe; ++x) {
1177     if ((*x).second) {
1178       std::cerr << "FACE LOOP ERROR: " << (*x).first.first << "-" << (*x).first.second << " : " << (*x).second << std::endl;
1179     }
1180   }
1181 }
1182
1183
1184
1185 /** 
1186  * 
1187  * 
1188  * @param a 
1189  * @param b 
1190  * @param vclass 
1191  * @param eclass 
1192  * @param a_face_loops 
1193  * @param b_face_loops 
1194  * @param a_edge_count 
1195  * @param b_edge_count 
1196  * @param hooks 
1197  */
1198 void carve::csg::CSG::calc(carve::mesh::MeshSet<3> *a,
1199                            const face_rtree_t *a_rtree,
1200                            carve::mesh::MeshSet<3> *b,
1201                            const face_rtree_t *b_rtree,
1202                            carve::csg::VertexClassification &vclass,
1203                            carve::csg::EdgeClassification &eclass,
1204                            carve::csg::FaceLoopList &a_face_loops,
1205                            carve::csg::FaceLoopList &b_face_loops,
1206                            size_t &a_edge_count,
1207                            size_t &b_edge_count) {
1208   detail::Data data;
1209
1210 #if defined(CARVE_DEBUG)
1211   std::cerr << "init" << std::endl;
1212 #endif
1213   init();
1214
1215   generateIntersections(a, a_rtree, b, b_rtree, data);
1216
1217 #if defined(CARVE_DEBUG)
1218   std::cerr << "intersectingFacePairs" << std::endl;
1219 #endif
1220   intersectingFacePairs(data);
1221
1222 #if defined(CARVE_DEBUG)
1223   std::cerr << "emap:" << std::endl;
1224   map_histogram(std::cerr, data.emap);
1225   std::cerr << "fmap:" << std::endl;
1226   map_histogram(std::cerr, data.fmap);
1227   std::cerr << "fmap_rev:" << std::endl;
1228   map_histogram(std::cerr, data.fmap_rev);
1229 #endif
1230
1231   // std::cerr << "removeCoplanarFaces" << std::endl;
1232   // fp_intersections.removeCoplanarFaces();
1233
1234 #if defined(CARVE_DEBUG) && defined(DEBUG_DRAW_OCTREE)
1235   HOOK(drawOctree(a->octree););
1236   HOOK(drawOctree(b->octree););
1237 #endif
1238
1239 #if defined(CARVE_DEBUG)
1240   std::cerr << "divideIntersectedEdges" << std::endl;
1241 #endif
1242   divideIntersectedEdges(data);
1243
1244 #if defined(CARVE_DEBUG)
1245   std::cerr << "makeFaceEdges" << std::endl;
1246 #endif
1247   // makeFaceEdges(data.face_split_edges, eclass, data.fmap, data.fmap_rev);
1248   makeFaceEdges(eclass, data);
1249
1250 #if defined(CARVE_DEBUG)
1251   std::cerr << "generateFaceLoops" << std::endl;
1252 #endif
1253   a_edge_count = generateFaceLoops(a, data, a_face_loops);
1254   b_edge_count = generateFaceLoops(b, data, b_face_loops);
1255
1256 #if defined(CARVE_DEBUG)
1257   std::cerr << "generated " << a_edge_count << " edges for poly a" << std::endl;
1258   std::cerr << "generated " << b_edge_count << " edges for poly b" << std::endl;
1259 #endif
1260
1261 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
1262   {
1263     std::string out("/tmp/a_split.ply");
1264     writePLY(out, faceLoopsToPolyhedron(a_face_loops), false);
1265   }
1266   {
1267     std::string out("/tmp/b_split.ply");
1268     writePLY(out, faceLoopsToPolyhedron(b_face_loops), false);
1269   }
1270 #endif
1271
1272   checkFaceLoopIntegrity(a_face_loops);
1273   checkFaceLoopIntegrity(b_face_loops);
1274
1275 #if defined(CARVE_DEBUG)
1276   std::cerr << "classify" << std::endl;
1277 #endif
1278   // initialize some classification information.
1279   for (std::vector<carve::mesh::MeshSet<3>::vertex_t>::iterator
1280          i = a->vertex_storage.begin(), e = a->vertex_storage.end(); i != e; ++i) {
1281     vclass[map_vertex(data.vmap, &(*i))].cls[0] = POINT_ON;
1282   }
1283   for (std::vector<carve::mesh::MeshSet<3>::vertex_t>::iterator
1284          i = b->vertex_storage.begin(), e = b->vertex_storage.end(); i != e; ++i) {
1285     vclass[map_vertex(data.vmap, &(*i))].cls[1] = POINT_ON;
1286   }
1287   for (VertexIntersections::const_iterator
1288          i = vertex_intersections.begin(), e = vertex_intersections.end(); i != e; ++i) {
1289     vclass[(*i).first] = PC2(POINT_ON, POINT_ON);
1290   }
1291
1292 #if defined(CARVE_DEBUG)
1293   std::cerr << data.divided_edges.size() << " edges are split" << std::endl;
1294   std::cerr << data.face_split_edges.size() << " faces are split" << std::endl;
1295
1296   std::cerr << "poly a: " << a_face_loops.size() << " face loops" << std::endl;
1297   std::cerr << "poly b: " << b_face_loops.size() << " face loops" << std::endl;
1298 #endif
1299
1300   // std::cerr << "OCTREE A:" << std::endl;
1301   // dump_octree_stats(a->octree.root, 0);
1302   // std::cerr << "OCTREE B:" << std::endl;
1303   // dump_octree_stats(b->octree.root, 0);
1304 }
1305
1306
1307
1308 /** 
1309  * 
1310  * 
1311  * @param shared_edges 
1312  * @param result_list 
1313  * @param shared_edge_ptr 
1314  */
1315 void returnSharedEdges(carve::csg::V2Set &shared_edges, 
1316                        std::list<carve::mesh::MeshSet<3> *> &result_list,
1317                        carve::csg::V2Set *shared_edge_ptr) {
1318   // need to convert shared edges to point into result
1319   typedef std::map<carve::geom3d::Vector, carve::mesh::MeshSet<3>::vertex_t *> remap_type;
1320   remap_type remap;
1321   for (std::list<carve::mesh::MeshSet<3> *>::iterator list_it =
1322          result_list.begin(); list_it != result_list.end(); list_it++) {
1323     carve::mesh::MeshSet<3> *result = *list_it;
1324     if (result) {
1325       for (std::vector<carve::mesh::MeshSet<3>::vertex_t>::iterator it =
1326              result->vertex_storage.begin(); it != result->vertex_storage.end(); it++) {
1327         remap.insert(std::make_pair((*it).v, &(*it)));
1328       }
1329     }
1330   }
1331   for (carve::csg::V2Set::iterator it = shared_edges.begin(); 
1332        it != shared_edges.end(); it++) {
1333     remap_type::iterator first_it = remap.find(((*it).first)->v);
1334     remap_type::iterator second_it = remap.find(((*it).second)->v);
1335     CARVE_ASSERT(first_it != remap.end() && second_it != remap.end());
1336     shared_edge_ptr->insert(std::make_pair(first_it->second, second_it->second));
1337   }
1338 }
1339
1340
1341
1342 /** 
1343  * 
1344  * 
1345  * @param a 
1346  * @param b 
1347  * @param collector 
1348  * @param hooks 
1349  * @param shared_edges_ptr 
1350  * @param classify_type 
1351  * 
1352  * @return 
1353  */
1354 carve::mesh::MeshSet<3> *carve::csg::CSG::compute(carve::mesh::MeshSet<3> *a,
1355                                                   carve::mesh::MeshSet<3> *b,
1356                                                   carve::csg::CSG::Collector &collector,
1357                                                   carve::csg::V2Set *shared_edges_ptr,
1358                                                   CLASSIFY_TYPE classify_type) {
1359   static carve::TimingName FUNC_NAME("CSG::compute");
1360   carve::TimingBlock block(FUNC_NAME);
1361
1362   VertexClassification vclass;
1363   EdgeClassification eclass;
1364
1365   FLGroupList a_loops_grouped;
1366   FLGroupList b_loops_grouped;
1367
1368   FaceLoopList a_face_loops;
1369   FaceLoopList b_face_loops;
1370
1371   size_t a_edge_count;
1372   size_t b_edge_count;
1373
1374   face_rtree_t *a_rtree = face_rtree_t::construct_STR(a->faceBegin(), a->faceEnd(), 4, 4);
1375   face_rtree_t *b_rtree = face_rtree_t::construct_STR(b->faceBegin(), b->faceEnd(), 4, 4);
1376
1377   {
1378     static carve::TimingName FUNC_NAME("CSG::compute - calc()");
1379     carve::TimingBlock block(FUNC_NAME);
1380     calc(a, a_rtree, b, b_rtree, vclass, eclass,a_face_loops, b_face_loops, a_edge_count, b_edge_count);
1381   }
1382
1383   detail::LoopEdges a_edge_map;
1384   detail::LoopEdges b_edge_map;
1385
1386   {
1387     static carve::TimingName FUNC_NAME("CSG::compute - makeEdgeMap()");
1388     carve::TimingBlock block(FUNC_NAME);
1389     makeEdgeMap(a_face_loops, a_edge_count, a_edge_map);
1390     makeEdgeMap(b_face_loops, b_edge_count, b_edge_map);
1391
1392   }
1393   
1394   {
1395     static carve::TimingName FUNC_NAME("CSG::compute - sortFaceLoopLists()");
1396     carve::TimingBlock block(FUNC_NAME);
1397     a_edge_map.sortFaceLoopLists();
1398     b_edge_map.sortFaceLoopLists();
1399   }
1400
1401   V2Set shared_edges;
1402   
1403   {
1404     static carve::TimingName FUNC_NAME("CSG::compute - findSharedEdges()");
1405     carve::TimingBlock block(FUNC_NAME);
1406     findSharedEdges(a_edge_map, b_edge_map, shared_edges);
1407   }
1408
1409   {
1410     static carve::TimingName FUNC_NAME("CSG::compute - groupFaceLoops()");
1411     carve::TimingBlock block(FUNC_NAME);
1412     groupFaceLoops(a, a_face_loops, a_edge_map, shared_edges, a_loops_grouped);
1413     groupFaceLoops(b, b_face_loops, b_edge_map, shared_edges, b_loops_grouped);
1414 #if defined(CARVE_DEBUG)
1415     std::cerr << "*** a_loops_grouped.size(): " << a_loops_grouped.size() << std::endl;
1416     std::cerr << "*** b_loops_grouped.size(): " << b_loops_grouped.size() << std::endl;
1417 #endif
1418   }
1419
1420 #if defined(CARVE_DEBUG) && defined(DEBUG_DRAW_GROUPS)
1421   {
1422     float n = 1.0 / (a_loops_grouped.size() + b_loops_grouped.size() + 1);
1423     float H = 0.0, S = 1.0, V = 1.0;
1424     float r, g, b;
1425     for (FLGroupList::const_iterator i = a_loops_grouped.begin(); i != a_loops_grouped.end(); ++i) {
1426       carve::colour::HSV2RGB(H, S, V, r, g, b); H += n;
1427       drawFaceLoopList((*i).face_loops, r, g, b, 1.0, r * .5, g * .5, b * .5, 1.0, true);
1428     }
1429     for (FLGroupList::const_iterator i = b_loops_grouped.begin(); i != b_loops_grouped.end(); ++i) {
1430       carve::colour::HSV2RGB(H, S, V, r, g, b); H += n;
1431       drawFaceLoopList((*i).face_loops, r, g, b, 1.0, r * .5, g * .5, b * .5, 1.0, true);
1432     }
1433
1434     for (FLGroupList::const_iterator i = a_loops_grouped.begin(); i != a_loops_grouped.end(); ++i) {
1435       drawFaceLoopListWireframe((*i).face_loops);
1436     }
1437     for (FLGroupList::const_iterator i = b_loops_grouped.begin(); i != b_loops_grouped.end(); ++i) {
1438       drawFaceLoopListWireframe((*i).face_loops);
1439     }
1440   }
1441 #endif
1442
1443   switch (classify_type) {
1444   case CLASSIFY_EDGE:
1445     classifyFaceGroupsEdge(shared_edges,
1446                            vclass,
1447                            a,
1448                            a_rtree,
1449                            a_loops_grouped,
1450                            a_edge_map,
1451                            b,
1452                            b_rtree,
1453                            b_loops_grouped,
1454                            b_edge_map,
1455                            collector);
1456     break;
1457   case CLASSIFY_NORMAL:
1458     classifyFaceGroups(shared_edges,
1459                        vclass,
1460                        a,
1461                        a_rtree,
1462                        a_loops_grouped,
1463                        a_edge_map,
1464                        b,
1465                        b_rtree,
1466                        b_loops_grouped,
1467                        b_edge_map,
1468                        collector);
1469     break;
1470   }
1471
1472   carve::mesh::MeshSet<3> *result = collector.done(hooks);
1473   if (result != NULL && shared_edges_ptr != NULL) {
1474     std::list<carve::mesh::MeshSet<3> *> result_list;
1475     result_list.push_back(result);
1476     returnSharedEdges(shared_edges, result_list, shared_edges_ptr);
1477   }
1478   return result;
1479 }
1480
1481
1482
1483 /** 
1484  * 
1485  * 
1486  * @param a 
1487  * @param b 
1488  * @param op 
1489  * @param hooks 
1490  * @param shared_edges 
1491  * @param classify_type 
1492  * 
1493  * @return 
1494  */
1495 carve::mesh::MeshSet<3> *carve::csg::CSG::compute(carve::mesh::MeshSet<3> *a,
1496                                                   carve::mesh::MeshSet<3> *b,
1497                                                   carve::csg::CSG::OP op,
1498                                                   carve::csg::V2Set *shared_edges,
1499                                                   CLASSIFY_TYPE classify_type) {
1500   Collector *coll = makeCollector(op, a, b);
1501   if (!coll) return NULL;
1502
1503   carve::mesh::MeshSet<3> *result = compute(a, b, *coll, shared_edges, classify_type);
1504      
1505   delete coll;
1506
1507   return result;
1508 }
1509
1510
1511
1512 /** 
1513  * 
1514  * 
1515  * @param closed 
1516  * @param open 
1517  * @param FaceClass 
1518  * @param result 
1519  * @param hooks 
1520  * @param shared_edges_ptr 
1521  * 
1522  * @return 
1523  */
1524 bool carve::csg::CSG::sliceAndClassify(carve::mesh::MeshSet<3> *closed,
1525                                        carve::mesh::MeshSet<3> *open,
1526                                        std::list<std::pair<FaceClass, carve::mesh::MeshSet<3> *> > &result,
1527                                        carve::csg::V2Set *shared_edges_ptr) {
1528   if (!closed->isClosed()) return false;
1529   carve::csg::VertexClassification vclass;
1530   carve::csg::EdgeClassification eclass;
1531
1532   carve::csg::FLGroupList a_loops_grouped;
1533   carve::csg::FLGroupList b_loops_grouped;
1534
1535   carve::csg::FaceLoopList a_face_loops;
1536   carve::csg::FaceLoopList b_face_loops;
1537
1538   size_t a_edge_count;
1539   size_t b_edge_count;
1540
1541   face_rtree_t *closed_rtree = face_rtree_t::construct_STR(closed->faceBegin(), closed->faceEnd(), 4, 4);
1542   face_rtree_t *open_rtree = face_rtree_t::construct_STR(open->faceBegin(), open->faceEnd(), 4, 4);
1543
1544   calc(closed, closed_rtree, open, open_rtree, vclass, eclass,a_face_loops, b_face_loops, a_edge_count, b_edge_count);
1545
1546   detail::LoopEdges a_edge_map;
1547   detail::LoopEdges b_edge_map;
1548
1549   makeEdgeMap(a_face_loops, a_edge_count, a_edge_map);
1550   makeEdgeMap(b_face_loops, b_edge_count, b_edge_map);
1551   
1552   carve::csg::V2Set shared_edges;
1553   
1554   findSharedEdges(a_edge_map, b_edge_map, shared_edges);
1555   
1556   groupFaceLoops(closed, a_face_loops, a_edge_map, shared_edges, a_loops_grouped);
1557   groupFaceLoops(open,   b_face_loops, b_edge_map, shared_edges, b_loops_grouped);
1558
1559   halfClassifyFaceGroups(shared_edges,
1560                          vclass,
1561                          closed,
1562                          closed_rtree,
1563                          a_loops_grouped,
1564                          a_edge_map,
1565                          open,
1566                          open_rtree,
1567                          b_loops_grouped,
1568                          b_edge_map,
1569                          result);
1570
1571   if (shared_edges_ptr != NULL) {
1572     std::list<carve::mesh::MeshSet<3> *> result_list;
1573     for (std::list<std::pair<FaceClass, carve::mesh::MeshSet<3> *> >::iterator it = result.begin(); it != result.end(); it++) {
1574       result_list.push_back(it->second);
1575     }
1576     returnSharedEdges(shared_edges, result_list, shared_edges_ptr);
1577   }
1578   return true;
1579 }
1580
1581
1582
1583 /** 
1584  * 
1585  * 
1586  * @param a 
1587  * @param b 
1588  * @param a_sliced 
1589  * @param b_sliced 
1590  * @param hooks 
1591  * @param shared_edges_ptr 
1592  */
1593 void carve::csg::CSG::slice(carve::mesh::MeshSet<3> *a,
1594                             carve::mesh::MeshSet<3> *b,
1595                             std::list<carve::mesh::MeshSet<3> *> &a_sliced,
1596                             std::list<carve::mesh::MeshSet<3> *> &b_sliced,
1597                             carve::csg::V2Set *shared_edges_ptr) {
1598   carve::csg::VertexClassification vclass;
1599   carve::csg::EdgeClassification eclass;
1600
1601   carve::csg::FLGroupList a_loops_grouped;
1602   carve::csg::FLGroupList b_loops_grouped;
1603
1604   carve::csg::FaceLoopList a_face_loops;
1605   carve::csg::FaceLoopList b_face_loops;
1606
1607   size_t a_edge_count;
1608   size_t b_edge_count;
1609
1610   face_rtree_t *a_rtree = face_rtree_t::construct_STR(a->faceBegin(), a->faceEnd(), 4, 4);
1611   face_rtree_t *b_rtree = face_rtree_t::construct_STR(b->faceBegin(), b->faceEnd(), 4, 4);
1612
1613   calc(a, a_rtree, b, b_rtree, vclass, eclass,a_face_loops, b_face_loops, a_edge_count, b_edge_count);
1614
1615   detail::LoopEdges a_edge_map;
1616   detail::LoopEdges b_edge_map;
1617       
1618   makeEdgeMap(a_face_loops, a_edge_count, a_edge_map);
1619   makeEdgeMap(b_face_loops, b_edge_count, b_edge_map);
1620   
1621   carve::csg::V2Set shared_edges;
1622   
1623   findSharedEdges(a_edge_map, b_edge_map, shared_edges);
1624   
1625   groupFaceLoops(a, a_face_loops, a_edge_map, shared_edges, a_loops_grouped);
1626   groupFaceLoops(b, b_face_loops, b_edge_map, shared_edges, b_loops_grouped);
1627
1628   for (carve::csg::FLGroupList::iterator
1629          i = a_loops_grouped.begin(), e = a_loops_grouped.end();
1630        i != e; ++i) {
1631     Collector *all = makeCollector(ALL, a, b);
1632     all->collect(&*i, hooks);
1633     a_sliced.push_back(all->done(hooks));
1634
1635     delete all;
1636   }
1637
1638   for (carve::csg::FLGroupList::iterator
1639          i = b_loops_grouped.begin(), e = b_loops_grouped.end();
1640        i != e; ++i) {
1641     Collector *all = makeCollector(ALL, a, b);
1642     all->collect(&*i, hooks);
1643     b_sliced.push_back(all->done(hooks));
1644
1645     delete all;
1646   }
1647   if (shared_edges_ptr != NULL) {
1648     std::list<carve::mesh::MeshSet<3> *> result_list;
1649     result_list.insert(result_list.end(), a_sliced.begin(), a_sliced.end());
1650     result_list.insert(result_list.end(), b_sliced.begin(), b_sliced.end());
1651     returnSharedEdges(shared_edges, result_list, shared_edges_ptr);
1652   }
1653 }
1654
1655
1656
1657 /** 
1658  * 
1659  * 
1660  */
1661 void carve::csg::CSG::init() {
1662   intersections.clear();
1663   vertex_intersections.clear();
1664   vertex_pool.reset();
1665 }