Carve booleans library integration
[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   if (std::min(eb->v1()->v.x, eb->v2()->v.x) - carve::EPSILON > va->v.x ||
651       std::max(eb->v1()->v.x, eb->v2()->v.x) + carve::EPSILON < va->v.x ||
652       std::min(eb->v1()->v.y, eb->v2()->v.y) - carve::EPSILON > va->v.y ||
653       std::max(eb->v1()->v.y, eb->v2()->v.y) + carve::EPSILON < va->v.y ||
654       std::min(eb->v1()->v.z, eb->v2()->v.z) - carve::EPSILON > va->v.z ||
655       std::max(eb->v1()->v.z, eb->v2()->v.z) + carve::EPSILON < va->v.z) {
656     return;
657   }
658
659   double a = cross(eb->v2()->v - eb->v1()->v, va->v - eb->v1()->v).length2();
660   double b = (eb->v2()->v - eb->v1()->v).length2();
661
662   if (a < b * carve::EPSILON2) {
663     // vertex-edge intersection
664     intersections.record(eb, va, va);
665     if (eb->rev) intersections.record(eb->rev, va, va);
666   }
667 }
668
669
670
671 void carve::csg::CSG::generateVertexEdgeIntersections(carve::mesh::MeshSet<3>::face_t *a,
672                                                       const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
673   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
674
675   ea = a->edge;
676   do {
677     for (size_t i = 0; i < b.size(); ++i) {
678       carve::mesh::MeshSet<3>::face_t *t = b[i];
679       eb = t->edge;
680       do {
681         _generateVertexEdgeIntersections(ea->v1(), eb);
682         eb = eb->next;
683       } while (eb != t->edge);
684     }
685     ea = ea->next;
686   } while (ea != a->edge);
687 }
688
689
690
691 void carve::csg::CSG::_generateEdgeEdgeIntersections(carve::mesh::MeshSet<3>::edge_t *ea,
692                                                      carve::mesh::MeshSet<3>::edge_t *eb) {
693   if (intersections.intersects(ea, eb)) {
694     return;
695   }
696
697   carve::mesh::MeshSet<3>::vertex_t *v1 = ea->v1(), *v2 = ea->v2();
698   carve::mesh::MeshSet<3>::vertex_t *v3 = eb->v1(), *v4 = eb->v2();
699
700   carve::geom::aabb<3> ea_aabb, eb_aabb;
701   ea_aabb.fit(v1->v, v2->v);
702   eb_aabb.fit(v3->v, v4->v);
703   if (ea_aabb.maxAxisSeparation(eb_aabb) > EPSILON) return;
704
705   carve::mesh::MeshSet<3>::vertex_t::vector_t p1, p2;
706   double mu1, mu2;
707
708   switch (carve::geom3d::rayRayIntersection(carve::geom3d::Ray(v2->v - v1->v, v1->v),
709                                             carve::geom3d::Ray(v4->v - v3->v, v3->v),
710                                             p1, p2, mu1, mu2)) {
711   case carve::RR_INTERSECTION: {
712     // edges intersect
713     if (mu1 >= 0.0 && mu1 <= 1.0 && mu2 >= 0.0 && mu2 <= 1.0) {
714       carve::mesh::MeshSet<3>::vertex_t *p = vertex_pool.get((p1 + p2) / 2.0);
715       intersections.record(ea, eb, p);
716       if (ea->rev) intersections.record(ea->rev, eb, p);
717       if (eb->rev) intersections.record(ea, eb->rev, p);
718       if (ea->rev && eb->rev) intersections.record(ea->rev, eb->rev, p);
719     }
720     break;
721   }
722   case carve::RR_PARALLEL: {
723     // edges parallel. any intersection of this type should have
724     // been handled by generateVertexEdgeIntersections().
725     break;
726   }
727   case carve::RR_DEGENERATE: {
728     throw carve::exception("degenerate edge");
729     break;
730   }
731   case carve::RR_NO_INTERSECTION: {
732     break;
733   }
734   }
735 }
736
737
738
739 void carve::csg::CSG::generateEdgeEdgeIntersections(carve::mesh::MeshSet<3>::face_t *a,
740                                                     const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
741   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
742
743   ea = a->edge;
744   do {
745     for (size_t i = 0; i < b.size(); ++i) {
746       carve::mesh::MeshSet<3>::face_t *t = b[i];
747       eb = t->edge;
748       do {
749         _generateEdgeEdgeIntersections(ea, eb);
750         eb = eb->next;
751       } while (eb != t->edge);
752     }
753     ea = ea->next;
754   } while (ea != a->edge);
755 }
756
757
758
759 void carve::csg::CSG::_generateVertexFaceIntersections(carve::mesh::MeshSet<3>::face_t *fa,
760                                                        carve::mesh::MeshSet<3>::edge_t *eb) {
761   if (intersections.intersects(eb->v1(), fa)) {
762     return;
763   }
764
765   double d1 = carve::geom::distance(fa->plane, eb->v1()->v);
766
767   if (fabs(d1) < carve::EPSILON &&
768       fa->containsPoint(eb->v1()->v)) {
769     intersections.record(eb->v1(), fa, eb->v1());
770   }
771 }
772
773
774
775 void carve::csg::CSG::generateVertexFaceIntersections(carve::mesh::MeshSet<3>::face_t *a,
776                                                       const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
777   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
778
779   for (size_t i = 0; i < b.size(); ++i) {
780     carve::mesh::MeshSet<3>::face_t *t = b[i];
781     eb = t->edge;
782     do {
783       _generateVertexFaceIntersections(a, eb);
784       eb = eb->next;
785     } while (eb != t->edge);
786   }
787 }
788
789
790
791 void carve::csg::CSG::_generateEdgeFaceIntersections(carve::mesh::MeshSet<3>::face_t *fa,
792                                                      carve::mesh::MeshSet<3>::edge_t *eb) {
793   if (intersections.intersects(eb, fa)) {
794     return;
795   }
796
797   carve::mesh::MeshSet<3>::vertex_t::vector_t _p;
798   if (fa->simpleLineSegmentIntersection(carve::geom3d::LineSegment(eb->v1()->v, eb->v2()->v), _p)) {
799     carve::mesh::MeshSet<3>::vertex_t *p = vertex_pool.get(_p);
800     intersections.record(eb, fa, p);
801     if (eb->rev) intersections.record(eb->rev, fa, p);
802   }
803 }
804
805
806
807 void carve::csg::CSG::generateEdgeFaceIntersections(carve::mesh::MeshSet<3>::face_t *a,
808                                                     const std::vector<carve::mesh::MeshSet<3>::face_t *> &b) {
809   carve::mesh::MeshSet<3>::edge_t *ea, *eb;
810
811   for (size_t i = 0; i < b.size(); ++i) {
812     carve::mesh::MeshSet<3>::face_t *t = b[i];
813     eb = t->edge;
814     do {
815       _generateEdgeFaceIntersections(a, eb);
816       eb = eb->next;
817     } while (eb != t->edge);
818   }
819 }
820
821
822
823 void carve::csg::CSG::generateIntersectionCandidates(carve::mesh::MeshSet<3> *a,
824                                                      const face_rtree_t *a_node,
825                                                      carve::mesh::MeshSet<3> *b,
826                                                      const face_rtree_t *b_node,
827                                                      face_pairs_t &face_pairs,
828                                                      bool descend_a) {
829   if (!a_node->bbox.intersects(b_node->bbox)) {
830     return;
831   }
832
833   if (a_node->child && (descend_a || !b_node->child)) {
834     for (face_rtree_t *node = a_node->child; node; node = node->sibling) {
835       generateIntersectionCandidates(a, node, b, b_node, face_pairs, false);
836     }
837   } else if (b_node->child) {
838     for (face_rtree_t *node = b_node->child; node; node = node->sibling) {
839       generateIntersectionCandidates(a, a_node, b, node, face_pairs, true);
840     }
841   } else {
842     for (size_t i = 0; i < a_node->data.size(); ++i) {
843       carve::mesh::MeshSet<3>::face_t *fa = a_node->data[i];
844       carve::geom::aabb<3> aabb_a = fa->getAABB();
845       if (aabb_a.maxAxisSeparation(b_node->bbox) > carve::EPSILON) continue;
846
847       for (size_t j = 0; j < b_node->data.size(); ++j) {
848         carve::mesh::MeshSet<3>::face_t *fb = b_node->data[j];
849         carve::geom::aabb<3> aabb_b = fb->getAABB();
850         if (aabb_b.maxAxisSeparation(aabb_a) > carve::EPSILON) continue;
851
852         std::pair<double, double> a_ra = fa->rangeInDirection(fa->plane.N, fa->edge->vert->v);
853         std::pair<double, double> b_ra = fb->rangeInDirection(fa->plane.N, fa->edge->vert->v);
854         if (carve::rangeSeparation(a_ra, b_ra) > carve::EPSILON) continue;
855
856         std::pair<double, double> a_rb = fa->rangeInDirection(fb->plane.N, fb->edge->vert->v);
857         std::pair<double, double> b_rb = fb->rangeInDirection(fb->plane.N, fb->edge->vert->v);
858         if (carve::rangeSeparation(a_rb, b_rb) > carve::EPSILON) continue;
859
860         if (!facesAreCoplanar(fa, fb)) {
861           face_pairs[fa].push_back(fb); 
862           face_pairs[fb].push_back(fa);
863         }
864      }
865     }
866   }
867 }
868
869
870
871
872 void carve::csg::CSG::generateIntersections(carve::mesh::MeshSet<3> *a,
873                                             const face_rtree_t *a_rtree,
874                                             carve::mesh::MeshSet<3> *b,
875                                             const face_rtree_t *b_rtree,
876                                             detail::Data &data) {
877   face_pairs_t face_pairs;
878   generateIntersectionCandidates(a, a_rtree, b, b_rtree, face_pairs);
879
880   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
881     carve::mesh::MeshSet<3>::face_t *f = (*i).first;
882     carve::mesh::MeshSet<3>::edge_t *e = f->edge;
883     do {
884       data.vert_to_edges[e->v1()].push_back(e);
885       e = e->next;
886     } while (e != f->edge);
887   }
888
889   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
890     generateVertexVertexIntersections((*i).first, (*i).second);
891   }
892
893   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
894     generateVertexEdgeIntersections((*i).first, (*i).second);
895   }
896
897   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
898     generateEdgeEdgeIntersections((*i).first, (*i).second);
899   }
900
901   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
902     generateVertexFaceIntersections((*i).first, (*i).second);
903   }
904
905   for (face_pairs_t::const_iterator i = face_pairs.begin(); i != face_pairs.end(); ++i) {
906     generateEdgeFaceIntersections((*i).first, (*i).second);
907   }
908
909
910 #if defined(CARVE_DEBUG)
911   std::cerr << "makeVertexIntersections" << std::endl;
912 #endif
913   makeVertexIntersections();
914
915 #if defined(CARVE_DEBUG)
916   std::cerr << "  intersections.size() " << intersections.size() << std::endl;
917   map_histogram(std::cerr, intersections);
918   std::cerr << "  vertex_intersections.size() " << vertex_intersections.size() << std::endl;
919   map_histogram(std::cerr, vertex_intersections);
920 #endif
921
922 #if defined(CARVE_DEBUG) && defined(DEBUG_DRAW_INTERSECTIONS)
923   HOOK(drawIntersections(vertex_intersections););
924 #endif
925
926 #if defined(CARVE_DEBUG)
927   std::cerr << "  intersections.size() " << intersections.size() << std::endl;
928   std::cerr << "  vertex_intersections.size() " << vertex_intersections.size() << std::endl;
929 #endif
930
931   // notify about intersections.
932   if (hooks.hasHook(Hooks::INTERSECTION_VERTEX_HOOK)) {
933     for (VertexIntersections::const_iterator i = vertex_intersections.begin();
934          i != vertex_intersections.end();
935          ++i) {
936       hooks.intersectionVertex((*i).first, (*i).second);
937     }
938   }
939
940   // from here on, only vertex_intersections is used for intersection
941   // information.
942
943   // intersections still contains the vertex_to_face map. maybe that
944   // should be moved out into another class.
945   static_cast<Intersections::super>(intersections).clear();
946 }
947
948
949
950 carve::csg::CSG::CSG() {
951 }
952
953
954
955 /** 
956  * \brief For each intersected edge, decompose into a set of vertex pairs representing an ordered set of edge fragments.
957  * 
958  * @tparam[in,out] data Internal intersection data. data.emap is used to produce data.divided_edges.
959  */
960 void carve::csg::CSG::divideIntersectedEdges(detail::Data &data) {
961   static carve::TimingName FUNC_NAME("CSG::divideIntersectedEdges()");
962   carve::TimingBlock block(FUNC_NAME);
963
964   for (detail::EVSMap::const_iterator i = data.emap.begin(), ei = data.emap.end(); i != ei; ++i) {
965     carve::mesh::MeshSet<3>::edge_t *edge = (*i).first;
966     const detail::EVSMap::mapped_type &vertices = (*i).second;
967     std::vector<carve::mesh::MeshSet<3>::vertex_t *> &verts = data.divided_edges[edge];
968     orderVertices(vertices.begin(), vertices.end(),
969                   edge->v2()->v - edge->v1()->v, edge->v1()->v,
970                   verts);
971   }
972 }
973
974
975
976 carve::csg::CSG::~CSG() {
977 }
978
979
980
981 void carve::csg::CSG::makeFaceEdges(carve::csg::EdgeClassification &eclass,
982                                     detail::Data &data) {
983   detail::FSet face_b_set;
984   for (detail::FVSMap::const_iterator
985          i = data.fmap.begin(), ie = data.fmap.end();
986        i != ie;
987        ++i) {
988     carve::mesh::MeshSet<3>::face_t *face_a = (*i).first;
989     const detail::FVSMap::mapped_type &face_a_intersections = ((*i).second);
990     face_b_set.clear();
991
992     // work out the set of faces from the opposing polyhedron that intersect face_a.
993     for (detail::FVSMap::mapped_type::const_iterator
994            j = face_a_intersections.begin(), je = face_a_intersections.end();
995          j != je;
996          ++j) {
997       for (detail::VFSMap::mapped_type::const_iterator
998              k = data.fmap_rev[*j].begin(), ke = data.fmap_rev[*j].end();
999            k != ke;
1000            ++k) {
1001         carve::mesh::MeshSet<3>::face_t *face_b = (*k);
1002         if (face_a != face_b && face_b->mesh->meshset != face_a->mesh->meshset) {
1003           face_b_set.insert(face_b);
1004         }
1005       }
1006     }
1007
1008     // run through each intersecting face.
1009     for (detail::FSet::const_iterator
1010            j = face_b_set.begin(), je = face_b_set.end();
1011          j != je;
1012          ++j) {
1013       carve::mesh::MeshSet<3>::face_t *face_b = (*j);
1014       const detail::FVSMap::mapped_type &face_b_intersections = (data.fmap[face_b]);
1015
1016       std::vector<carve::mesh::MeshSet<3>::vertex_t *> vertices;
1017       vertices.reserve(std::min(face_a_intersections.size(), face_b_intersections.size()));
1018
1019       // record the points of intersection between face_a and face_b
1020       std::set_intersection(face_a_intersections.begin(),
1021                             face_a_intersections.end(),
1022                             face_b_intersections.begin(),
1023                             face_b_intersections.end(),
1024                             std::back_inserter(vertices));
1025
1026 #if defined(CARVE_DEBUG)
1027       std::cerr << "face pair: "
1028                 << face_a << ":" << face_b
1029                 << " N(verts) " << vertices.size() << std::endl;
1030       for (std::vector<carve::mesh::MeshSet<3>::vertex_t *>::const_iterator i = vertices.begin(), e = vertices.end(); i != e; ++i) {
1031         std::cerr << (*i) << " " << (*i)->v << " ("
1032                   << carve::geom::distance(face_a->plane, (*i)->v) << ","
1033                   << carve::geom::distance(face_b->plane, (*i)->v) << ")"
1034                   << std::endl;
1035         //CARVE_ASSERT(carve::geom3d::distance(face_a->plane_eqn, *(*i)) < EPSILON);
1036         //CARVE_ASSERT(carve::geom3d::distance(face_b->plane_eqn, *(*i)) < EPSILON);
1037       }
1038 #endif
1039
1040       // if there are two points of intersection, then the added edge is simple to determine.
1041       if (vertices.size() == 2) {
1042         carve::mesh::MeshSet<3>::vertex_t *v1 = vertices[0];
1043         carve::mesh::MeshSet<3>::vertex_t *v2 = vertices[1];
1044         carve::geom3d::Vector c = (v1->v + v2->v) / 2;
1045
1046         // determine whether the midpoint of the implied edge is contained in face_a and face_b
1047
1048 #if defined(CARVE_DEBUG)
1049         std::cerr << "face_a->nVertices() = " << face_a->nVertices() << " face_a->containsPointInProjection(c) = " << face_a->containsPointInProjection(c) << std::endl;
1050         std::cerr << "face_b->nVertices() = " << face_b->nVertices() << " face_b->containsPointInProjection(c) = " << face_b->containsPointInProjection(c) << std::endl;
1051 #endif
1052
1053         if (face_a->containsPointInProjection(c) && face_b->containsPointInProjection(c)) {
1054 #if defined(CARVE_DEBUG)
1055           std::cerr << "adding edge: " << v1 << "-" << v2 << std::endl;
1056 #if defined(DEBUG_DRAW_FACE_EDGES)
1057           HOOK(drawEdge(v1, v2, 1, 1, 1, 1, 1, 1, 1, 1, 2.0););
1058 #endif
1059 #endif
1060           // record the edge, with class information.
1061           if (v1 > v2) std::swap(v1, v2);
1062           eclass[ordered_edge(v1, v2)] = carve::csg::EC2(carve::csg::EDGE_ON, carve::csg::EDGE_ON);
1063           data.face_split_edges[face_a].insert(std::make_pair(v1, v2));
1064           data.face_split_edges[face_b].insert(std::make_pair(v1, v2));
1065         }
1066         continue;
1067       }
1068
1069       // otherwise, it's more complex.
1070       carve::geom3d::Vector base, dir;
1071       std::vector<carve::mesh::MeshSet<3>::vertex_t *> ordered;
1072
1073       // skip coplanar edges. this simplifies the resulting
1074       // mesh. eventually all coplanar face regions of two polyhedra
1075       // must reach a point where they are no longer coplanar (or the
1076       // polyhedra are identical).
1077       if (!facesAreCoplanar(face_a, face_b)) {
1078         // order the intersection vertices (they must lie along a
1079         // vector, as the faces aren't coplanar).
1080         selectOrderingProjection(vertices.begin(), vertices.end(), dir, base);
1081         orderVertices(vertices.begin(), vertices.end(), dir, base, ordered);
1082
1083         // for each possible edge in the ordering, test the midpoint,
1084         // and record if it's contained in face_a and face_b.
1085         for (int k = 0, ke = (int)ordered.size() - 1; k < ke; ++k) {
1086           carve::mesh::MeshSet<3>::vertex_t *v1 = ordered[k];
1087           carve::mesh::MeshSet<3>::vertex_t *v2 = ordered[k + 1];
1088           carve::geom3d::Vector c = (v1->v + v2->v) / 2;
1089
1090 #if defined(CARVE_DEBUG)
1091           std::cerr << "testing edge: " << v1 << "-" << v2 << " at " << c << std::endl;
1092           std::cerr << "a: " << face_a->containsPointInProjection(c) << " b: " << face_b->containsPointInProjection(c) << std::endl;
1093           std::cerr << "face_a->containsPointInProjection(c): " << face_a->containsPointInProjection(c) << std::endl;
1094           std::cerr << "face_b->containsPointInProjection(c): " << face_b->containsPointInProjection(c) << std::endl;
1095 #endif
1096
1097           if (face_a->containsPointInProjection(c) && face_b->containsPointInProjection(c)) {
1098 #if defined(CARVE_DEBUG)
1099             std::cerr << "adding edge: " << v1 << "-" << v2 << std::endl;
1100 #if defined(DEBUG_DRAW_FACE_EDGES)
1101             HOOK(drawEdge(v1, v2, .5, .5, .5, 1, .5, .5, .5, 1, 2.0););
1102 #endif
1103 #endif
1104             // record the edge, with class information.
1105             if (v1 > v2) std::swap(v1, v2);
1106             eclass[ordered_edge(v1, v2)] = carve::csg::EC2(carve::csg::EDGE_ON, carve::csg::EDGE_ON);
1107             data.face_split_edges[face_a].insert(std::make_pair(v1, v2));
1108             data.face_split_edges[face_b].insert(std::make_pair(v1, v2));
1109           }
1110         }
1111       }
1112     }
1113   }
1114
1115
1116 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
1117   {
1118     V2Set edges;
1119     for (detail::FV2SMap::const_iterator i = data.face_split_edges.begin(); i != data.face_split_edges.end(); ++i) {
1120       edges.insert((*i).second.begin(), (*i).second.end());
1121     }
1122
1123     detail::VSet vertices;
1124     for (V2Set::const_iterator i = edges.begin(); i != edges.end(); ++i) {
1125       vertices.insert((*i).first);
1126       vertices.insert((*i).second);
1127     }
1128
1129     carve::line::PolylineSet intersection_graph;
1130     intersection_graph.vertices.resize(vertices.size());
1131     std::map<const carve::mesh::MeshSet<3>::vertex_t *, size_t> vmap;
1132
1133     size_t j = 0;
1134     for (detail::VSet::const_iterator i = vertices.begin(); i != vertices.end(); ++i) {
1135       intersection_graph.vertices[j].v = (*i)->v;
1136       vmap[(*i)] = j++;
1137     }
1138
1139     for (V2Set::const_iterator i = edges.begin(); i != edges.end(); ++i) {
1140       size_t line[2];
1141       line[0] = vmap[(*i).first];
1142       line[1] = vmap[(*i).second];
1143       intersection_graph.addPolyline(false, line, line + 2);
1144     }
1145
1146     std::string out("/tmp/intersection-edges.ply");
1147     ::writePLY(out, &intersection_graph, true);
1148   }
1149 #endif
1150 }
1151
1152
1153
1154 /** 
1155  * 
1156  * 
1157  * @param fll 
1158  */
1159 static void checkFaceLoopIntegrity(carve::csg::FaceLoopList &fll) {
1160   static carve::TimingName FUNC_NAME("CSG::checkFaceLoopIntegrity()");
1161   carve::TimingBlock block(FUNC_NAME);
1162
1163   std::unordered_map<carve::csg::V2, int> counts;
1164   for (carve::csg::FaceLoop *fl = fll.head; fl; fl = fl->next) {
1165     std::vector<carve::mesh::MeshSet<3>::vertex_t *> &loop = (fl->vertices);
1166     carve::mesh::MeshSet<3>::vertex_t *v1, *v2;
1167     v1 = loop[loop.size() - 1];
1168     for (unsigned i = 0; i < loop.size(); ++i) {
1169       v2 = loop[i];
1170       if (v1 < v2) {
1171         counts[std::make_pair(v1, v2)]++;
1172       } else {
1173         counts[std::make_pair(v2, v1)]--;
1174       }
1175       v1 = v2;
1176     }
1177   }
1178   for (std::unordered_map<carve::csg::V2, int>::const_iterator
1179          x = counts.begin(), xe = counts.end(); x != xe; ++x) {
1180     if ((*x).second) {
1181       std::cerr << "FACE LOOP ERROR: " << (*x).first.first << "-" << (*x).first.second << " : " << (*x).second << std::endl;
1182     }
1183   }
1184 }
1185
1186
1187
1188 /** 
1189  * 
1190  * 
1191  * @param a 
1192  * @param b 
1193  * @param vclass 
1194  * @param eclass 
1195  * @param a_face_loops 
1196  * @param b_face_loops 
1197  * @param a_edge_count 
1198  * @param b_edge_count 
1199  * @param hooks 
1200  */
1201 void carve::csg::CSG::calc(carve::mesh::MeshSet<3> *a,
1202                            const face_rtree_t *a_rtree,
1203                            carve::mesh::MeshSet<3> *b,
1204                            const face_rtree_t *b_rtree,
1205                            carve::csg::VertexClassification &vclass,
1206                            carve::csg::EdgeClassification &eclass,
1207                            carve::csg::FaceLoopList &a_face_loops,
1208                            carve::csg::FaceLoopList &b_face_loops,
1209                            size_t &a_edge_count,
1210                            size_t &b_edge_count) {
1211   detail::Data data;
1212
1213 #if defined(CARVE_DEBUG)
1214   std::cerr << "init" << std::endl;
1215 #endif
1216   init();
1217
1218   generateIntersections(a, a_rtree, b, b_rtree, data);
1219
1220 #if defined(CARVE_DEBUG)
1221   std::cerr << "intersectingFacePairs" << std::endl;
1222 #endif
1223   intersectingFacePairs(data);
1224
1225 #if defined(CARVE_DEBUG)
1226   std::cerr << "emap:" << std::endl;
1227   map_histogram(std::cerr, data.emap);
1228   std::cerr << "fmap:" << std::endl;
1229   map_histogram(std::cerr, data.fmap);
1230   std::cerr << "fmap_rev:" << std::endl;
1231   map_histogram(std::cerr, data.fmap_rev);
1232 #endif
1233
1234   // std::cerr << "removeCoplanarFaces" << std::endl;
1235   // fp_intersections.removeCoplanarFaces();
1236
1237 #if defined(CARVE_DEBUG) && defined(DEBUG_DRAW_OCTREE)
1238   HOOK(drawOctree(a->octree););
1239   HOOK(drawOctree(b->octree););
1240 #endif
1241
1242 #if defined(CARVE_DEBUG)
1243   std::cerr << "divideIntersectedEdges" << std::endl;
1244 #endif
1245   divideIntersectedEdges(data);
1246
1247 #if defined(CARVE_DEBUG)
1248   std::cerr << "makeFaceEdges" << std::endl;
1249 #endif
1250   // makeFaceEdges(data.face_split_edges, eclass, data.fmap, data.fmap_rev);
1251   makeFaceEdges(eclass, data);
1252
1253 #if defined(CARVE_DEBUG)
1254   std::cerr << "generateFaceLoops" << std::endl;
1255 #endif
1256   a_edge_count = generateFaceLoops(a, data, a_face_loops);
1257   b_edge_count = generateFaceLoops(b, data, b_face_loops);
1258
1259 #if defined(CARVE_DEBUG)
1260   std::cerr << "generated " << a_edge_count << " edges for poly a" << std::endl;
1261   std::cerr << "generated " << b_edge_count << " edges for poly b" << std::endl;
1262 #endif
1263
1264 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
1265   {
1266     std::string out("/tmp/a_split.ply");
1267     writePLY(out, faceLoopsToPolyhedron(a_face_loops), false);
1268   }
1269   {
1270     std::string out("/tmp/b_split.ply");
1271     writePLY(out, faceLoopsToPolyhedron(b_face_loops), false);
1272   }
1273 #endif
1274
1275   checkFaceLoopIntegrity(a_face_loops);
1276   checkFaceLoopIntegrity(b_face_loops);
1277
1278 #if defined(CARVE_DEBUG)
1279   std::cerr << "classify" << std::endl;
1280 #endif
1281   // initialize some classification information.
1282   for (std::vector<carve::mesh::MeshSet<3>::vertex_t>::iterator
1283          i = a->vertex_storage.begin(), e = a->vertex_storage.end(); i != e; ++i) {
1284     vclass[map_vertex(data.vmap, &(*i))].cls[0] = POINT_ON;
1285   }
1286   for (std::vector<carve::mesh::MeshSet<3>::vertex_t>::iterator
1287          i = b->vertex_storage.begin(), e = b->vertex_storage.end(); i != e; ++i) {
1288     vclass[map_vertex(data.vmap, &(*i))].cls[1] = POINT_ON;
1289   }
1290   for (VertexIntersections::const_iterator
1291          i = vertex_intersections.begin(), e = vertex_intersections.end(); i != e; ++i) {
1292     vclass[(*i).first] = PC2(POINT_ON, POINT_ON);
1293   }
1294
1295 #if defined(CARVE_DEBUG)
1296   std::cerr << data.divided_edges.size() << " edges are split" << std::endl;
1297   std::cerr << data.face_split_edges.size() << " faces are split" << std::endl;
1298
1299   std::cerr << "poly a: " << a_face_loops.size() << " face loops" << std::endl;
1300   std::cerr << "poly b: " << b_face_loops.size() << " face loops" << std::endl;
1301 #endif
1302
1303   // std::cerr << "OCTREE A:" << std::endl;
1304   // dump_octree_stats(a->octree.root, 0);
1305   // std::cerr << "OCTREE B:" << std::endl;
1306   // dump_octree_stats(b->octree.root, 0);
1307 }
1308
1309
1310
1311 /** 
1312  * 
1313  * 
1314  * @param shared_edges 
1315  * @param result_list 
1316  * @param shared_edge_ptr 
1317  */
1318 void returnSharedEdges(carve::csg::V2Set &shared_edges, 
1319                        std::list<carve::mesh::MeshSet<3> *> &result_list,
1320                        carve::csg::V2Set *shared_edge_ptr) {
1321   // need to convert shared edges to point into result
1322   typedef std::map<carve::geom3d::Vector, carve::mesh::MeshSet<3>::vertex_t *> remap_type;
1323   remap_type remap;
1324   for (std::list<carve::mesh::MeshSet<3> *>::iterator list_it =
1325          result_list.begin(); list_it != result_list.end(); list_it++) {
1326     carve::mesh::MeshSet<3> *result = *list_it;
1327     if (result) {
1328       for (std::vector<carve::mesh::MeshSet<3>::vertex_t>::iterator it =
1329              result->vertex_storage.begin(); it != result->vertex_storage.end(); it++) {
1330         remap.insert(std::make_pair((*it).v, &(*it)));
1331       }
1332     }
1333   }
1334   for (carve::csg::V2Set::iterator it = shared_edges.begin(); 
1335        it != shared_edges.end(); it++) {
1336     remap_type::iterator first_it = remap.find(((*it).first)->v);
1337     remap_type::iterator second_it = remap.find(((*it).second)->v);
1338     CARVE_ASSERT(first_it != remap.end() && second_it != remap.end());
1339     shared_edge_ptr->insert(std::make_pair(first_it->second, second_it->second));
1340   }
1341 }
1342
1343
1344
1345 /** 
1346  * 
1347  * 
1348  * @param a 
1349  * @param b 
1350  * @param collector 
1351  * @param hooks 
1352  * @param shared_edges_ptr 
1353  * @param classify_type 
1354  * 
1355  * @return 
1356  */
1357 carve::mesh::MeshSet<3> *carve::csg::CSG::compute(carve::mesh::MeshSet<3> *a,
1358                                                   carve::mesh::MeshSet<3> *b,
1359                                                   carve::csg::CSG::Collector &collector,
1360                                                   carve::csg::V2Set *shared_edges_ptr,
1361                                                   CLASSIFY_TYPE classify_type) {
1362   static carve::TimingName FUNC_NAME("CSG::compute");
1363   carve::TimingBlock block(FUNC_NAME);
1364
1365   VertexClassification vclass;
1366   EdgeClassification eclass;
1367
1368   FLGroupList a_loops_grouped;
1369   FLGroupList b_loops_grouped;
1370
1371   FaceLoopList a_face_loops;
1372   FaceLoopList b_face_loops;
1373
1374   size_t a_edge_count;
1375   size_t b_edge_count;
1376
1377   face_rtree_t *a_rtree = face_rtree_t::construct_STR(a->faceBegin(), a->faceEnd(), 4, 4);
1378   face_rtree_t *b_rtree = face_rtree_t::construct_STR(b->faceBegin(), b->faceEnd(), 4, 4);
1379
1380   {
1381     static carve::TimingName FUNC_NAME("CSG::compute - calc()");
1382     carve::TimingBlock block(FUNC_NAME);
1383     calc(a, a_rtree, b, b_rtree, vclass, eclass,a_face_loops, b_face_loops, a_edge_count, b_edge_count);
1384   }
1385
1386   detail::LoopEdges a_edge_map;
1387   detail::LoopEdges b_edge_map;
1388
1389   {
1390     static carve::TimingName FUNC_NAME("CSG::compute - makeEdgeMap()");
1391     carve::TimingBlock block(FUNC_NAME);
1392     makeEdgeMap(a_face_loops, a_edge_count, a_edge_map);
1393     makeEdgeMap(b_face_loops, b_edge_count, b_edge_map);
1394
1395   }
1396   
1397   {
1398     static carve::TimingName FUNC_NAME("CSG::compute - sortFaceLoopLists()");
1399     carve::TimingBlock block(FUNC_NAME);
1400     a_edge_map.sortFaceLoopLists();
1401     b_edge_map.sortFaceLoopLists();
1402   }
1403
1404   V2Set shared_edges;
1405   
1406   {
1407     static carve::TimingName FUNC_NAME("CSG::compute - findSharedEdges()");
1408     carve::TimingBlock block(FUNC_NAME);
1409     findSharedEdges(a_edge_map, b_edge_map, shared_edges);
1410   }
1411
1412   {
1413     static carve::TimingName FUNC_NAME("CSG::compute - groupFaceLoops()");
1414     carve::TimingBlock block(FUNC_NAME);
1415     groupFaceLoops(a, a_face_loops, a_edge_map, shared_edges, a_loops_grouped);
1416     groupFaceLoops(b, b_face_loops, b_edge_map, shared_edges, b_loops_grouped);
1417 #if defined(CARVE_DEBUG)
1418     std::cerr << "*** a_loops_grouped.size(): " << a_loops_grouped.size() << std::endl;
1419     std::cerr << "*** b_loops_grouped.size(): " << b_loops_grouped.size() << std::endl;
1420 #endif
1421   }
1422
1423 #if defined(CARVE_DEBUG) && defined(DEBUG_DRAW_GROUPS)
1424   {
1425     float n = 1.0 / (a_loops_grouped.size() + b_loops_grouped.size() + 1);
1426     float H = 0.0, S = 1.0, V = 1.0;
1427     float r, g, b;
1428     for (FLGroupList::const_iterator i = a_loops_grouped.begin(); i != a_loops_grouped.end(); ++i) {
1429       carve::colour::HSV2RGB(H, S, V, r, g, b); H += n;
1430       drawFaceLoopList((*i).face_loops, r, g, b, 1.0, r * .5, g * .5, b * .5, 1.0, true);
1431     }
1432     for (FLGroupList::const_iterator i = b_loops_grouped.begin(); i != b_loops_grouped.end(); ++i) {
1433       carve::colour::HSV2RGB(H, S, V, r, g, b); H += n;
1434       drawFaceLoopList((*i).face_loops, r, g, b, 1.0, r * .5, g * .5, b * .5, 1.0, true);
1435     }
1436
1437     for (FLGroupList::const_iterator i = a_loops_grouped.begin(); i != a_loops_grouped.end(); ++i) {
1438       drawFaceLoopListWireframe((*i).face_loops);
1439     }
1440     for (FLGroupList::const_iterator i = b_loops_grouped.begin(); i != b_loops_grouped.end(); ++i) {
1441       drawFaceLoopListWireframe((*i).face_loops);
1442     }
1443   }
1444 #endif
1445
1446   switch (classify_type) {
1447   case CLASSIFY_EDGE:
1448     classifyFaceGroupsEdge(shared_edges,
1449                            vclass,
1450                            a,
1451                            a_rtree,
1452                            a_loops_grouped,
1453                            a_edge_map,
1454                            b,
1455                            b_rtree,
1456                            b_loops_grouped,
1457                            b_edge_map,
1458                            collector);
1459     break;
1460   case CLASSIFY_NORMAL:
1461     classifyFaceGroups(shared_edges,
1462                        vclass,
1463                        a,
1464                        a_rtree,
1465                        a_loops_grouped,
1466                        a_edge_map,
1467                        b,
1468                        b_rtree,
1469                        b_loops_grouped,
1470                        b_edge_map,
1471                        collector);
1472     break;
1473   }
1474
1475   carve::mesh::MeshSet<3> *result = collector.done(hooks);
1476   if (result != NULL && shared_edges_ptr != NULL) {
1477     std::list<carve::mesh::MeshSet<3> *> result_list;
1478     result_list.push_back(result);
1479     returnSharedEdges(shared_edges, result_list, shared_edges_ptr);
1480   }
1481   return result;
1482 }
1483
1484
1485
1486 /** 
1487  * 
1488  * 
1489  * @param a 
1490  * @param b 
1491  * @param op 
1492  * @param hooks 
1493  * @param shared_edges 
1494  * @param classify_type 
1495  * 
1496  * @return 
1497  */
1498 carve::mesh::MeshSet<3> *carve::csg::CSG::compute(carve::mesh::MeshSet<3> *a,
1499                                                   carve::mesh::MeshSet<3> *b,
1500                                                   carve::csg::CSG::OP op,
1501                                                   carve::csg::V2Set *shared_edges,
1502                                                   CLASSIFY_TYPE classify_type) {
1503   Collector *coll = makeCollector(op, a, b);
1504   if (!coll) return NULL;
1505
1506   carve::mesh::MeshSet<3> *result = compute(a, b, *coll, shared_edges, classify_type);
1507      
1508   delete coll;
1509
1510   return result;
1511 }
1512
1513
1514
1515 /** 
1516  * 
1517  * 
1518  * @param closed 
1519  * @param open 
1520  * @param FaceClass 
1521  * @param result 
1522  * @param hooks 
1523  * @param shared_edges_ptr 
1524  * 
1525  * @return 
1526  */
1527 bool carve::csg::CSG::sliceAndClassify(carve::mesh::MeshSet<3> *closed,
1528                                        carve::mesh::MeshSet<3> *open,
1529                                        std::list<std::pair<FaceClass, carve::mesh::MeshSet<3> *> > &result,
1530                                        carve::csg::V2Set *shared_edges_ptr) {
1531   if (!closed->isClosed()) return false;
1532   carve::csg::VertexClassification vclass;
1533   carve::csg::EdgeClassification eclass;
1534
1535   carve::csg::FLGroupList a_loops_grouped;
1536   carve::csg::FLGroupList b_loops_grouped;
1537
1538   carve::csg::FaceLoopList a_face_loops;
1539   carve::csg::FaceLoopList b_face_loops;
1540
1541   size_t a_edge_count;
1542   size_t b_edge_count;
1543
1544   face_rtree_t *closed_rtree = face_rtree_t::construct_STR(closed->faceBegin(), closed->faceEnd(), 4, 4);
1545   face_rtree_t *open_rtree = face_rtree_t::construct_STR(open->faceBegin(), open->faceEnd(), 4, 4);
1546
1547   calc(closed, closed_rtree, open, open_rtree, vclass, eclass,a_face_loops, b_face_loops, a_edge_count, b_edge_count);
1548
1549   detail::LoopEdges a_edge_map;
1550   detail::LoopEdges b_edge_map;
1551
1552   makeEdgeMap(a_face_loops, a_edge_count, a_edge_map);
1553   makeEdgeMap(b_face_loops, b_edge_count, b_edge_map);
1554   
1555   carve::csg::V2Set shared_edges;
1556   
1557   findSharedEdges(a_edge_map, b_edge_map, shared_edges);
1558   
1559   groupFaceLoops(closed, a_face_loops, a_edge_map, shared_edges, a_loops_grouped);
1560   groupFaceLoops(open,   b_face_loops, b_edge_map, shared_edges, b_loops_grouped);
1561
1562   halfClassifyFaceGroups(shared_edges,
1563                          vclass,
1564                          closed,
1565                          closed_rtree,
1566                          a_loops_grouped,
1567                          a_edge_map,
1568                          open,
1569                          open_rtree,
1570                          b_loops_grouped,
1571                          b_edge_map,
1572                          result);
1573
1574   if (shared_edges_ptr != NULL) {
1575     std::list<carve::mesh::MeshSet<3> *> result_list;
1576     for (std::list<std::pair<FaceClass, carve::mesh::MeshSet<3> *> >::iterator it = result.begin(); it != result.end(); it++) {
1577       result_list.push_back(it->second);
1578     }
1579     returnSharedEdges(shared_edges, result_list, shared_edges_ptr);
1580   }
1581   return true;
1582 }
1583
1584
1585
1586 /** 
1587  * 
1588  * 
1589  * @param a 
1590  * @param b 
1591  * @param a_sliced 
1592  * @param b_sliced 
1593  * @param hooks 
1594  * @param shared_edges_ptr 
1595  */
1596 void carve::csg::CSG::slice(carve::mesh::MeshSet<3> *a,
1597                             carve::mesh::MeshSet<3> *b,
1598                             std::list<carve::mesh::MeshSet<3> *> &a_sliced,
1599                             std::list<carve::mesh::MeshSet<3> *> &b_sliced,
1600                             carve::csg::V2Set *shared_edges_ptr) {
1601   carve::csg::VertexClassification vclass;
1602   carve::csg::EdgeClassification eclass;
1603
1604   carve::csg::FLGroupList a_loops_grouped;
1605   carve::csg::FLGroupList b_loops_grouped;
1606
1607   carve::csg::FaceLoopList a_face_loops;
1608   carve::csg::FaceLoopList b_face_loops;
1609
1610   size_t a_edge_count;
1611   size_t b_edge_count;
1612
1613   face_rtree_t *a_rtree = face_rtree_t::construct_STR(a->faceBegin(), a->faceEnd(), 4, 4);
1614   face_rtree_t *b_rtree = face_rtree_t::construct_STR(b->faceBegin(), b->faceEnd(), 4, 4);
1615
1616   calc(a, a_rtree, b, b_rtree, vclass, eclass,a_face_loops, b_face_loops, a_edge_count, b_edge_count);
1617
1618   detail::LoopEdges a_edge_map;
1619   detail::LoopEdges b_edge_map;
1620       
1621   makeEdgeMap(a_face_loops, a_edge_count, a_edge_map);
1622   makeEdgeMap(b_face_loops, b_edge_count, b_edge_map);
1623   
1624   carve::csg::V2Set shared_edges;
1625   
1626   findSharedEdges(a_edge_map, b_edge_map, shared_edges);
1627   
1628   groupFaceLoops(a, a_face_loops, a_edge_map, shared_edges, a_loops_grouped);
1629   groupFaceLoops(b, b_face_loops, b_edge_map, shared_edges, b_loops_grouped);
1630
1631   for (carve::csg::FLGroupList::iterator
1632          i = a_loops_grouped.begin(), e = a_loops_grouped.end();
1633        i != e; ++i) {
1634     Collector *all = makeCollector(ALL, a, b);
1635     all->collect(&*i, hooks);
1636     a_sliced.push_back(all->done(hooks));
1637
1638     delete all;
1639   }
1640
1641   for (carve::csg::FLGroupList::iterator
1642          i = b_loops_grouped.begin(), e = b_loops_grouped.end();
1643        i != e; ++i) {
1644     Collector *all = makeCollector(ALL, a, b);
1645     all->collect(&*i, hooks);
1646     b_sliced.push_back(all->done(hooks));
1647
1648     delete all;
1649   }
1650   if (shared_edges_ptr != NULL) {
1651     std::list<carve::mesh::MeshSet<3> *> result_list;
1652     result_list.insert(result_list.end(), a_sliced.begin(), a_sliced.end());
1653     result_list.insert(result_list.end(), b_sliced.begin(), b_sliced.end());
1654     returnSharedEdges(shared_edges, result_list, shared_edges_ptr);
1655   }
1656 }
1657
1658
1659
1660 /** 
1661  * 
1662  * 
1663  */
1664 void carve::csg::CSG::init() {
1665   intersections.clear();
1666   vertex_intersections.clear();
1667   vertex_pool.reset();
1668 }