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