Silence some annoying warnings when doing full build with strict flags
[blender.git] / extern / carve / lib / intersect_face_division.cpp
1 // Begin License:
2 // Copyright (C) 2006-2014 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 either the GNU General
8 // Public License version 2 or 3 (at your option) as published by the
9 // Free Software Foundation and appearing in the files LICENSE.GPL2
10 // and LICENSE.GPL3 included in the packaging of 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/polyline.hpp>
24 #include <carve/debug_hooks.hpp>
25 #include <carve/timing.hpp>
26 #include <carve/triangulator.hpp>
27
28 #include <list>
29 #include <set>
30 #include <iostream>
31
32 #include <algorithm>
33
34 #include "csg_detail.hpp"
35 #include "csg_data.hpp"
36
37 #include "intersect_common.hpp"
38
39
40
41 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
42 void writePLY(const std::string &out_file, const carve::line::PolylineSet *lines, bool ascii);
43 #endif
44
45
46
47 namespace {
48
49
50
51 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
52   template<typename iter_t>
53   void dumpFacesAndHoles(iter_t f_begin, iter_t f_end,
54                          iter_t h_begin, iter_t h_end,
55                          const std::string &fname) {
56     std::cerr << "dumping " << std::distance(f_begin, f_end) << " faces, " << std::distance(h_begin, h_end) << " holes." << std::endl;
57     std::map<carve::mesh::MeshSet<3>::vertex_t *, size_t> v_included;
58
59     for (iter_t i = f_begin; i != f_end; ++i) {
60       for (size_t j = 0; j < (*i).size(); ++j) {
61         if (v_included.find((*i)[j]) == v_included.end()) {
62           size_t &p = v_included[(*i)[j]];
63           p = v_included.size() - 1;
64         }
65       }
66     }
67
68     for (iter_t i = h_begin; i != h_end; ++i) {
69       for (size_t j = 0; j < (*i).size(); ++j) {
70         if (v_included.find((*i)[j]) == v_included.end()) {
71           size_t &p = v_included[(*i)[j]];
72           p = v_included.size() - 1;
73         }
74       }
75     }
76
77     carve::line::PolylineSet fh;
78     fh.vertices.resize(v_included.size());
79     for (std::map<carve::mesh::MeshSet<3>::vertex_t *, size_t>::const_iterator
80            i = v_included.begin(); i != v_included.end(); ++i) {
81       fh.vertices[(*i).second].v = (*i).first->v;
82     }
83
84     {
85       std::vector<size_t> connected;
86       for (iter_t i = f_begin; i != f_end; ++i) {
87         connected.clear();
88         for (size_t j = 0; j < (*i).size(); ++j) {
89           connected.push_back(v_included[(*i)[j]]);
90         }
91         fh.addPolyline(true, connected.begin(), connected.end());
92       }
93       for (iter_t i = h_begin; i != h_end; ++i) {
94         connected.clear();
95         for (size_t j = 0; j < (*i).size(); ++j) {
96           connected.push_back(v_included[(*i)[j]]);
97         }
98         fh.addPolyline(true, connected.begin(), connected.end());
99       }
100     }
101
102     ::writePLY(fname, &fh, true);
103   }
104 #endif
105
106
107
108   template<typename T>
109   void populateVectorFromList(std::list<T> &l, std::vector<T> &v) {
110     v.clear();
111     v.reserve(l.size());
112     for (typename std::list<T>::iterator i = l.begin(); i != l.end(); ++i) {
113       v.push_back(T());
114       std::swap(*i, v.back());
115     }
116     l.clear();
117   }
118
119   template<typename T>
120   void populateListFromVector(std::vector<T> &v, std::list<T> &l) {
121     l.clear();
122     for (size_t i = 0; i < v.size(); ++i) {
123       l.push_back(T());
124       std::swap(v[i], l.back());
125     }
126     v.clear();
127   }
128
129
130
131   struct GraphEdge {
132     GraphEdge *next;
133     GraphEdge *prev;
134     GraphEdge *loop_next;
135     carve::mesh::MeshSet<3>::vertex_t *src;
136     carve::mesh::MeshSet<3>::vertex_t *tgt;
137     double ang;
138     int visited;
139
140     GraphEdge(carve::mesh::MeshSet<3>::vertex_t *_src, carve::mesh::MeshSet<3>::vertex_t *_tgt) :
141       next(NULL), prev(NULL), loop_next(NULL),
142       src(_src), tgt(_tgt),
143       ang(0.0), visited(-1) {
144     }
145   };
146
147
148
149   struct GraphEdges {
150     GraphEdge *edges;
151     carve::geom2d::P2 proj;
152
153     GraphEdges() : edges(NULL), proj() {
154     }
155   };
156
157
158
159   struct Graph {
160     typedef std::unordered_map<carve::mesh::MeshSet<3>::vertex_t *, GraphEdges> graph_t;
161
162     graph_t graph;
163
164     Graph() : graph() {
165     }
166
167     ~Graph() {
168       int c = 0;
169
170       GraphEdge *edge;
171       for (graph_t::iterator i = graph.begin(), e =  graph.end(); i != e; ++i) {
172         edge = (*i).second.edges;
173         while (edge) {
174           GraphEdge *temp = edge;
175           ++c;
176           edge = edge->next;
177           delete temp;
178         }
179       }
180
181       if (c) {
182         std::cerr << "warning: "
183                   << c
184                   << " edges should have already been removed at graph destruction time"
185                   << std::endl;
186       }
187     }
188
189     const carve::geom2d::P2 &projection(carve::mesh::MeshSet<3>::vertex_t *v) const {
190       graph_t::const_iterator i = graph.find(v);
191       CARVE_ASSERT(i != graph.end());
192       return (*i).second.proj;
193     }
194
195     void computeProjection(carve::mesh::MeshSet<3>::face_t *face) {
196       for (graph_t::iterator i = graph.begin(), e =  graph.end(); i != e; ++i) {
197         (*i).second.proj = face->project((*i).first->v);
198       }
199       for (graph_t::iterator i = graph.begin(), e =  graph.end(); i != e; ++i) {
200         for (GraphEdge *e = (*i).second.edges; e; e = e->next) {
201           e->ang = carve::math::ANG(carve::geom2d::atan2(projection(e->tgt) - projection(e->src)));
202         }
203       }
204     }
205
206     void print(std::ostream &out, const carve::csg::VertexIntersections *vi) const {
207       for (graph_t::const_iterator i = graph.begin(), e =  graph.end(); i != e; ++i) {
208         out << (*i).first << (*i).first->v << '(' << projection((*i).first).x << ',' << projection((*i).first).y << ") :";
209         for (const GraphEdge *e = (*i).second.edges; e; e = e->next) {
210           out << ' ' << e->tgt << e->tgt->v << '(' << projection(e->tgt).x << ',' << projection(e->tgt).y << ')';
211         }
212         out << std::endl;
213         if (vi) {
214           carve::csg::VertexIntersections::const_iterator j = vi->find((*i).first);
215           if (j != vi->end()) {
216             out << "   (int) ";
217             for (carve::csg::IObjPairSet::const_iterator
218                    k = (*j).second.begin(), ke = (*j).second.end(); k != ke; ++k) {
219               if ((*k).first < (*k).second) {
220                 out << (*k).first << ".." << (*k).second << "; ";
221               }
222             }
223             out << std::endl;
224           }
225         }
226       }
227     }
228
229     void addEdge(carve::mesh::MeshSet<3>::vertex_t *v1, carve::mesh::MeshSet<3>::vertex_t *v2) {
230       GraphEdges &edges = graph[v1];
231       GraphEdge *edge = new GraphEdge(v1, v2);
232       if (edges.edges) edges.edges->prev = edge;
233       edge->next = edges.edges;
234       edges.edges = edge;
235     }
236
237     void removeEdge(GraphEdge *edge) {
238       if (edge->prev != NULL) {
239         edge->prev->next = edge->next;
240       } else {
241         if (edge->next != NULL) {
242           GraphEdges &edges = (graph[edge->src]);
243           edges.edges = edge->next;
244         } else {
245           graph.erase(edge->src);
246         }
247       }
248       if (edge->next != NULL) {
249         edge->next->prev = edge->prev;
250       }
251       delete edge;
252     }
253
254     bool empty() const {
255       return graph.size() == 0;
256     }
257
258     GraphEdge *pickStartEdge() {
259       // Try and find a vertex from which there is only one outbound edge. Won't always succeed.
260       for (graph_t::iterator i = graph.begin(); i != graph.end(); ++i) {
261         GraphEdges &ge = i->second;
262         if (ge.edges->next == NULL) {
263           return ge.edges;
264         }
265       }
266       return (*graph.begin()).second.edges;
267     }
268
269     GraphEdge *outboundEdges(carve::mesh::MeshSet<3>::vertex_t *v) {
270       return graph[v].edges;
271     }
272   };
273
274
275
276   /** 
277    * \brief Take a set of new edges and split a face based upon those edges.
278    * 
279    * @param[in] face The face to be split.
280    * @param[in] edges 
281    * @param[out] face_loops Output list of face loops
282    * @param[out] hole_loops Output list of hole loops
283    * @param vi 
284    */
285   static void splitFace(carve::mesh::MeshSet<3>::face_t *face,
286                         const carve::csg::V2Set &edges,
287                         std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &face_loops,
288                         std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &hole_loops,
289                         const carve::csg::VertexIntersections & /* vi */) {
290     Graph graph;
291
292     for (carve::csg::V2Set::const_iterator
293            i = edges.begin(), e = edges.end();
294          i != e;
295          ++i) {
296       carve::mesh::MeshSet<3>::vertex_t *v1 = ((*i).first), *v2 = ((*i).second);
297       if (carve::geom::equal(v1->v, v2->v)) std::cerr << "WARNING! " << v1->v << "==" << v2->v << std::endl;
298       graph.addEdge(v1, v2);
299     }
300
301     graph.computeProjection(face);
302
303     while (!graph.empty()) {
304       GraphEdge *edge;
305       GraphEdge *start;
306       start = edge = graph.pickStartEdge();
307
308       edge->visited = 0;
309
310       int len = 0;
311
312       for (;;) {
313         double in_ang = M_PI + edge->ang;
314         if (in_ang > M_TWOPI) in_ang -= M_TWOPI;
315
316         GraphEdge *opts;
317         GraphEdge *out = NULL;
318         double best = M_TWOPI + 1.0;
319
320         for (opts = graph.outboundEdges(edge->tgt); opts; opts = opts->next) {
321           if (opts->tgt == edge->src) {
322             if (out == NULL && opts->next == NULL) out = opts;
323           } else {
324             double out_ang = carve::math::ANG(in_ang - opts->ang);
325
326             if (out == NULL || out_ang < best) {
327               out = opts;
328               best = out_ang;
329             }
330           }
331         }
332
333         CARVE_ASSERT(out != NULL);
334
335         edge->loop_next = out;
336
337         if (out->visited >= 0) {
338           while (start != out) {
339             GraphEdge *e = start;
340             start = start->loop_next;
341             e->loop_next = NULL;
342             e->visited = -1;
343           }
344           len = edge->visited - out->visited + 1;
345           break;
346         }
347
348         out->visited = edge->visited + 1;
349         edge = out;
350       }
351
352       std::vector<carve::mesh::MeshSet<3>::vertex_t *> loop(len);
353       std::vector<carve::geom2d::P2> projected(len);
354
355       edge = start;
356       for (int i = 0; i < len; ++i) {
357         GraphEdge *next = edge->loop_next;
358         loop[i] = edge->src;
359         projected[i] = graph.projection(edge->src);
360         graph.removeEdge(edge);
361         edge = next;
362       }
363
364       CARVE_ASSERT(edge == start);
365
366       if (carve::geom2d::signedArea(projected) < 0) {
367         face_loops.push_back(std::vector<carve::mesh::MeshSet<3>::vertex_t *>());
368         face_loops.back().swap(loop);
369       } else {
370         hole_loops.push_back(std::vector<carve::mesh::MeshSet<3>::vertex_t *>());
371         hole_loops.back().swap(loop);
372       }
373     }
374   }
375
376
377
378   /** 
379    * \brief Determine the relationship between a face loop and a hole loop.
380    * 
381    * Determine whether a face and hole share an edge, or a vertex,
382    * or do not touch. Find a hole vertex that is not part of the
383    * face, and a hole,face vertex pair that are coincident, if such
384    * a pair exists.
385    *
386    * @param[in] f A face loop.
387    * @param[in] f_sort A vector indexing \a f in address order
388    * @param[in] h A hole loop.
389    * @param[in] h_sort A vector indexing \a h in address order
390    * @param[out] f_idx Index of a face vertex that is shared with the hole.
391    * @param[out] h_idx Index of the hole vertex corresponding to \a f_idx.
392    * @param[out] unmatched_h_idx Index of a hole vertex that is not part of the face.
393    * @param[out] shares_vertex Boolean indicating that the face and the hole share a vertex.
394    * @param[out] shares_edge Boolean indicating that the face and the hole share an edge.
395    */
396   static void compareFaceLoopAndHoleLoop(const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &f,
397                                          const std::vector<unsigned> &f_sort,
398                                          const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &h,
399                                          const std::vector<unsigned> &h_sort,
400                                          unsigned &f_idx,
401                                          unsigned &h_idx,
402                                          int &unmatched_h_idx,
403                                          bool &shares_vertex,
404                                          bool &shares_edge) {
405     const size_t F = f.size();
406     const size_t H = h.size();
407
408     shares_vertex = shares_edge = false;
409     unmatched_h_idx = -1;
410
411     unsigned I, J;
412     for (I = J = 0; I < F && J < H;) {
413       unsigned i = f_sort[I], j = h_sort[J];
414       if (f[i] == h[j]) {
415         shares_vertex = true;
416         f_idx = i;
417         h_idx = j;
418         if (f[(i + F - 1) % F] == h[(j + 1) % H]) {
419           shares_edge = true;
420         }
421         carve::mesh::MeshSet<3>::vertex_t *t = f[i];
422         do { ++I; } while (I < F && f[f_sort[I]] == t);
423         do { ++J; } while (J < H && h[h_sort[J]] == t);
424       } else if (f[i] < h[j]) {
425         ++I;
426       } else {
427         unmatched_h_idx = j;
428         ++J;
429       }
430     }
431     if (J < H) {
432       unmatched_h_idx = h_sort[J];
433     }
434   }
435
436
437
438   /** 
439    * \brief Compute an embedding for a set of face loops and hole loops.
440    *
441    * Because face and hole loops may be contained within each other,
442    * it must be determined which hole loops are directly contained
443    * within a face loop.
444    * 
445    * @param[in] face The face from which these face and hole loops derive.
446    * @param[in] face_loops 
447    * @param[in] hole_loops 
448    * @param[out] containing_faces     A vector which for each hole loop
449    *                                  lists the indices of the face
450    *                                  loops it is containined in.
451    * @param[out] hole_shared_vertices A map from a face,hole pair to
452    *                                  a shared vertex pair.
453    */
454   static void computeContainment(carve::mesh::MeshSet<3>::face_t *face,
455                                  std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &face_loops,
456                                  std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &hole_loops,
457                                  std::vector<std::vector<int> > &containing_faces,
458                                  std::map<int, std::map<int, std::pair<unsigned, unsigned> > > &hole_shared_vertices) {
459 #if defined(CARVE_DEBUG)
460     std::cerr << "input: "
461               << face_loops.size() << "faces, "
462               << hole_loops.size() << "holes."
463               << std::endl;
464 #endif
465
466     std::vector<std::vector<carve::geom2d::P2> > face_loops_projected, hole_loops_projected;
467     std::vector<carve::geom::aabb<2> > face_loop_aabb, hole_loop_aabb;
468     std::vector<std::vector<unsigned> > face_loops_sorted, hole_loops_sorted;
469
470     std::vector<double> face_loop_areas, hole_loop_areas;
471
472     face_loops_projected.resize(face_loops.size());
473     face_loops_sorted.resize(face_loops.size());
474     face_loop_aabb.resize(face_loops.size());
475     face_loop_areas.resize(face_loops.size());
476
477     hole_loops_projected.resize(hole_loops.size());
478     hole_loops_sorted.resize(hole_loops.size());
479     hole_loop_aabb.resize(hole_loops.size());
480     hole_loop_areas.resize(hole_loops.size());
481
482     // produce a projection of each face loop onto a 2D plane, and an
483     // index vector which sorts vertices by address.
484     for (size_t m = 0; m < face_loops.size(); ++m) {
485       const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &f_loop = (face_loops[m]);
486       face_loops_projected[m].reserve(f_loop.size());
487       face_loops_sorted[m].reserve(f_loop.size());
488       for (size_t n = 0; n < f_loop.size(); ++n) {
489         face_loops_projected[m].push_back(face->project(f_loop[n]->v));
490         face_loops_sorted[m].push_back(n);
491       }
492       face_loop_areas.push_back(carve::geom2d::signedArea(face_loops_projected[m]));
493
494       std::sort(face_loops_sorted[m].begin(), face_loops_sorted[m].end(), 
495                 carve::make_index_sort(face_loops[m].begin()));
496       face_loop_aabb[m].fit(face_loops_projected[m].begin(), face_loops_projected[m].end());
497     }
498
499     // produce a projection of each hole loop onto a 2D plane, and an
500     // index vector which sorts vertices by address.
501     for (size_t m = 0; m < hole_loops.size(); ++m) {
502       const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &h_loop = (hole_loops[m]);
503       hole_loops_projected[m].reserve(h_loop.size());
504       hole_loops_projected[m].reserve(h_loop.size());
505       for (size_t n = 0; n < h_loop.size(); ++n) {
506         hole_loops_projected[m].push_back(face->project(h_loop[n]->v));
507         hole_loops_sorted[m].push_back(n);
508       }
509       hole_loop_areas.push_back(carve::geom2d::signedArea(hole_loops_projected[m]));
510
511       std::sort(hole_loops_sorted[m].begin(), hole_loops_sorted[m].end(), 
512                 carve::make_index_sort(hole_loops[m].begin()));
513       hole_loop_aabb[m].fit(hole_loops_projected[m].begin(), hole_loops_projected[m].end());
514     }
515
516     containing_faces.resize(hole_loops.size());
517
518     for (unsigned i = 0; i < hole_loops.size(); ++i) {
519
520       for (unsigned j = 0; j < face_loops.size(); ++j) {
521         if (!face_loop_aabb[j].completelyContains(hole_loop_aabb[i])) {
522 #if defined(CARVE_DEBUG)
523           std::cerr << "face: " << j
524                     << " hole: " << i
525                     << " skipped test (aabb fail)"
526                     << std::endl;
527 #endif
528           continue;
529         }
530
531         unsigned f_idx, h_idx;
532         int unmatched_h_idx;
533         bool shares_vertex, shares_edge;
534         compareFaceLoopAndHoleLoop(face_loops[j],
535                                    face_loops_sorted[j],
536                                    hole_loops[i],
537                                    hole_loops_sorted[i],
538                                    f_idx, h_idx,
539                                    unmatched_h_idx,
540                                    shares_vertex,
541                                    shares_edge);
542
543 #if defined(CARVE_DEBUG)
544         std::cerr << "face: " << j
545                   << " hole: " << i
546                   << " shares_vertex: " << shares_vertex
547                   << " shares_edge: " << shares_edge
548                   << std::endl;
549 #endif
550
551         carve::geom3d::Vector test = hole_loops[i][0]->v;
552         carve::geom2d::P2 test_p = face->project(test);
553
554         if (shares_vertex) {
555           hole_shared_vertices[i][j] = std::make_pair(h_idx, f_idx);
556           // Hole touches face. Should be able to connect it up
557           // trivially. Still need to record its containment, so that
558           // the assignment below works.
559           if (unmatched_h_idx != -1) {
560 #if defined(CARVE_DEBUG)
561             std::cerr << "using unmatched vertex: " << unmatched_h_idx << std::endl;
562 #endif
563             test = hole_loops[i][unmatched_h_idx]->v;
564             test_p = face->project(test);
565           } else {
566             // XXX: hole shares ALL vertices with face. Pick a point
567             // internal to the projected poly.
568             if (shares_edge) {
569               // Hole shares edge with face => face can't contain hole.
570               continue;
571             }
572
573             // XXX: how is this possible? Doesn't share an edge, but
574             // also doesn't have any vertices that are not in
575             // common. Degenerate hole?
576
577             // XXX: come up with a test case for this.
578             CARVE_FAIL("implement me");
579           }
580         }
581
582
583         // XXX: use loop area to avoid some point-in-poly tests? Loop
584         // area is faster, but not sure which is more robust.
585         if (carve::geom2d::pointInPolySimple(face_loops_projected[j], test_p)) {
586 #if defined(CARVE_DEBUG)
587           std::cerr << "contains: " << i << " - " << j << std::endl;
588 #endif
589           containing_faces[i].push_back(j);
590         } else {
591 #if defined(CARVE_DEBUG)
592           std::cerr << "does not contain: " << i << " - " << j << std::endl;
593 #endif
594         }
595       }
596
597 #if defined(CARVE_DEBUG)
598       if (containing_faces[i].size() == 0) {
599         //HOOK(drawFaceLoopWireframe(hole_loops[i], face->normal, 1.0, 0.0, 0.0, 1.0););
600         std::cerr << "hole loop: ";
601         for (unsigned j = 0; j < hole_loops[i].size(); ++j) {
602           std::cerr << " " << hole_loops[i][j] << ":" << hole_loops[i][j]->v;
603         }
604         std::cerr << std::endl;
605         for (unsigned j = 0; j < face_loops.size(); ++j) {
606           //HOOK(drawFaceLoopWireframe(face_loops[j], face->normal, 0.0, 1.0, 0.0, 1.0););
607         }
608       }
609 #endif
610
611       // CARVE_ASSERT(containing_faces[i].size() >= 1);
612     }
613   }
614
615
616
617   /** 
618    * \brief Merge face loops and hole loops to produce a set of face loops without holes.
619    * 
620    * @param[in] face The face from which these face loops derive.
621    * @param[in,out] f_loops A list of face loops.
622    * @param[in] h_loops A list of hole loops to be incorporated into face loops.
623    */
624   static void mergeFacesAndHoles(carve::mesh::MeshSet<3>::face_t *face,
625                                  std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &f_loops,
626                                  std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &h_loops,
627                                  carve::csg::CSG::Hooks & /* hooks */) {
628     std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > face_loops;
629     std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > hole_loops;
630
631     std::vector<std::vector<int> > containing_faces;
632     std::map<int, std::map<int, std::pair<unsigned, unsigned> > > hole_shared_vertices;
633
634 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
635     dumpFacesAndHoles(f_loops.begin(), f_loops.end(), h_loops.begin(), h_loops.end(), "/tmp/pre_merge.ply");
636 #endif
637
638     {
639       // move input face and hole loops to temp vectors.
640       size_t m;
641       face_loops.resize(f_loops.size());
642       m = 0;
643       for (std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> >::iterator
644              i = f_loops.begin(), ie = f_loops.end();
645            i != ie;
646            ++i, ++m) {
647         face_loops[m].swap((*i));
648       }
649
650       hole_loops.resize(h_loops.size());
651       m = 0;
652       for (std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> >::iterator
653              i = h_loops.begin(), ie = h_loops.end();
654            i != ie;
655            ++i, ++m) {
656         hole_loops[m].swap((*i));
657       }
658       f_loops.clear();
659       h_loops.clear();
660     }
661
662     // work out the embedding of holes and faces.
663     computeContainment(face, face_loops, hole_loops, containing_faces, hole_shared_vertices);
664
665     int unassigned = (int)hole_loops.size();
666
667     std::vector<std::vector<int> > face_holes;
668     face_holes.resize(face_loops.size());
669
670     for (unsigned i = 0; i < containing_faces.size(); ++i) {
671       if (containing_faces[i].size() == 0) {
672         std::map<int, std::map<int, std::pair<unsigned, unsigned> > >::iterator it = hole_shared_vertices.find(i);
673         if (it != hole_shared_vertices.end()) {
674           std::map<int, std::pair<unsigned, unsigned> >::iterator it2 = (*it).second.begin();
675           int f = (*it2).first;
676           unsigned h_idx = (*it2).second.first;
677           unsigned f_idx = (*it2).second.second;
678
679           // patch the hole into the face directly. because
680           // f_loop[f_idx] == h_loop[h_idx], we don't need to
681           // duplicate the f_loop vertex.
682
683           std::vector<carve::mesh::MeshSet<3>::vertex_t *> &f_loop = face_loops[f];
684           std::vector<carve::mesh::MeshSet<3>::vertex_t *> &h_loop = hole_loops[i];
685
686           f_loop.insert(f_loop.begin() + f_idx + 1, h_loop.size(), NULL);
687
688           unsigned p = f_idx + 1;
689           for (unsigned a = h_idx + 1; a < h_loop.size(); ++a, ++p) {
690             f_loop[p] = h_loop[a];
691           }
692           for (unsigned a = 0; a <= h_idx; ++a, ++p) {
693             f_loop[p] = h_loop[a];
694           }
695
696 #if defined(CARVE_DEBUG)
697           std::cerr << "hook face " << f << " to hole " << i << "(vertex)" << std::endl;
698 #endif
699         } else {
700           std::cerr << "uncontained hole loop does not share vertices with any face loop!" << std::endl;
701         }
702         unassigned--;
703       }
704     }
705
706
707     // work out which holes are directly contained within which faces.
708     while (unassigned) {
709       std::set<int> removed;
710
711       for (unsigned i = 0; i < containing_faces.size(); ++i) {
712         if (containing_faces[i].size() == 1) {
713           int f = containing_faces[i][0];
714           face_holes[f].push_back(i);
715 #if defined(CARVE_DEBUG)
716           std::cerr << "hook face " << f << " to hole " << i << std::endl;
717 #endif
718           removed.insert(f);
719           unassigned--;
720         }
721       }
722
723       if (!removed.size())
724         throw carve::exception("Failed to merge holes");
725
726       for (std::set<int>::iterator f = removed.begin(); f != removed.end(); ++f) {
727         for (unsigned i = 0; i < containing_faces.size(); ++i) {
728           containing_faces[i].erase(std::remove(containing_faces[i].begin(),
729                                                 containing_faces[i].end(),
730                                                 *f),
731                                     containing_faces[i].end());
732         }
733       }
734     }
735
736 #if 0
737     // use old templated projection code to patch holes into faces.
738     for (unsigned i = 0; i < face_loops.size(); ++i) {
739       std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > face_hole_loops;
740       face_hole_loops.resize(face_holes[i].size());
741       for (unsigned j = 0; j < face_holes[i].size(); ++j) {
742         face_hole_loops[j].swap(hole_loops[face_holes[i][j]]);
743       }
744       if (face_hole_loops.size()) {
745
746         f_loops.push_back(carve::triangulate::incorporateHolesIntoPolygon(
747           carve::mesh::MeshSet<3>::face_t::projection_mapping(face->project),
748           face_loops[i],
749           face_hole_loops));
750       } else {
751         f_loops.push_back(face_loops[i]);
752       }
753     }
754
755 #else
756     // use new 2d-only hole patching code.
757     for (size_t i = 0; i < face_loops.size(); ++i) {
758       if (!face_holes[i].size()) {
759         f_loops.push_back(face_loops[i]);
760         continue;
761       }
762
763       std::vector<std::vector<carve::geom2d::P2> > projected_poly;
764       projected_poly.resize(face_holes[i].size() + 1);
765       projected_poly[0].reserve(face_loops[i].size());
766       for (size_t j = 0; j < face_loops[i].size(); ++j) {
767         projected_poly[0].push_back(face->project(face_loops[i][j]->v));
768       }
769       for (size_t j = 0; j < face_holes[i].size(); ++j) {
770         projected_poly[j+1].reserve(hole_loops[face_holes[i][j]].size());
771         for (size_t k = 0; k < hole_loops[face_holes[i][j]].size(); ++k) {
772           projected_poly[j+1].push_back(face->project(hole_loops[face_holes[i][j]][k]->v));
773         }
774       }
775
776       std::vector<std::pair<size_t, size_t> > result = carve::triangulate::incorporateHolesIntoPolygon(projected_poly);
777
778       f_loops.push_back(std::vector<carve::mesh::MeshSet<3>::vertex_t *>());
779       std::vector<carve::mesh::MeshSet<3>::vertex_t *> &out = f_loops.back();
780       out.reserve(result.size());
781       for (size_t j = 0; j < result.size(); ++j) {
782         if (result[j].first == 0) {
783           out.push_back(face_loops[i][result[j].second]);
784         } else {
785           out.push_back(hole_loops[face_holes[i][result[j].first-1]][result[j].second]);
786         }
787       }
788     }
789 #endif
790 #if defined(CARVE_DEBUG_WRITE_PLY_DATA)
791     dumpFacesAndHoles(f_loops.begin(), f_loops.end(), h_loops.begin(), h_loops.end(), "/tmp/post_merge.ply");
792 #endif
793
794   }
795
796
797
798   /** 
799    * \brief Assemble the base loop for a face.
800    *
801    * The base loop is the original face loop, including vertices
802    * created by intersections crossing any of its edges.
803    * 
804    * @param[in] face The face to process.
805    * @param[in] vmap 
806    * @param[in] face_split_edges 
807    * @param[in] divided_edges A mapping from edge pointer to sets of
808    *            ordered vertices corrsponding to the intersection points
809    *            on that edge.
810    * @param[out] base_loop A vector of the vertices of the base loop.
811    */
812   static bool assembleBaseLoop(carve::mesh::MeshSet<3>::face_t *face,
813                                const carve::csg::detail::Data &data,
814                                std::vector<carve::mesh::MeshSet<3>::vertex_t *> &base_loop,
815                                carve::csg::CSG::Hooks &hooks) {
816     base_loop.clear();
817
818     // XXX: assumes that face->edges is in the same order as
819     // face->vertices. (Which it is)
820     carve::mesh::MeshSet<3>::edge_t *e = face->edge;
821     size_t e_idx = 0;
822     bool face_edge_intersected = false;
823     do {
824       base_loop.push_back(carve::csg::map_vertex(data.vmap, e->vert));
825
826       carve::csg::detail::EVVMap::const_iterator ev = data.divided_edges.find(e);
827
828       if (ev != data.divided_edges.end()) {
829         const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &ev_vec = ((*ev).second);
830
831         for (size_t k = 0, ke = ev_vec.size(); k < ke;) {
832           base_loop.push_back(ev_vec[k++]);
833         }
834
835         if (ev_vec.size() && hooks.hasHook(carve::csg::CSG::Hooks::EDGE_DIVISION_HOOK)) {
836           carve::mesh::MeshSet<3>::vertex_t *v1 = e->vert;
837           carve::mesh::MeshSet<3>::vertex_t *v2;
838           for (size_t k = 0, ke = ev_vec.size(); k < ke;) {
839             v2 = ev_vec[k++];
840             hooks.edgeDivision(e, e_idx, v1, v2);
841             v1 = v2;
842           }
843           v2 = e->v2();
844           hooks.edgeDivision(e, e_idx, v1, v2);
845         }
846
847         face_edge_intersected = true;
848       }
849       e = e->next;
850       ++e_idx;
851     } while (e != face->edge);
852
853     return face_edge_intersected;
854   }
855
856
857
858   // the crossing_data structure holds temporary information regarding
859   // paths, and their relationship to the loop of edges that forms the
860   // face perimeter.
861   struct crossing_data {
862     std::vector<carve::mesh::MeshSet<3>::vertex_t *> *path;
863     size_t edge_idx[2];
864
865     crossing_data(std::vector<carve::mesh::MeshSet<3>::vertex_t *> *p, size_t e1, size_t e2) : path(p) {
866       edge_idx[0] = e1; edge_idx[1] = e2;
867     }
868
869     bool operator<(const crossing_data &c) const {
870       // the sort order for paths is in order of increasing initial
871       // position on the edge loop, but decreasing final position.
872       return edge_idx[0] < c.edge_idx[0] || (edge_idx[0] == c.edge_idx[0] && edge_idx[1] > c.edge_idx[1]);
873     }
874   };
875
876
877
878   bool processCrossingEdges(carve::mesh::MeshSet<3>::face_t *face,
879                             const carve::csg::VertexIntersections &vertex_intersections,
880                             carve::csg::CSG::Hooks &hooks,
881                             std::vector<carve::mesh::MeshSet<3>::vertex_t *> &base_loop,
882                             std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &paths,
883                             std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &face_loops_out) {
884     const size_t N = base_loop.size();
885     std::vector<crossing_data> endpoint_indices;
886
887     endpoint_indices.reserve(paths.size());
888
889     for (size_t i = 0; i < paths.size(); ++i) {
890       endpoint_indices.push_back(crossing_data(&paths[i], N, N));
891     }
892
893     // Step 1:
894     // locate endpoints of paths on the base loop.
895     for (size_t i = 0; i < N; ++i) {
896       for (size_t j = 0; j < paths.size(); ++j) {
897         // test beginning of path.
898         if (paths[j].front() == base_loop[i]) {
899           if (endpoint_indices[j].edge_idx[0] == N) {
900             endpoint_indices[j].edge_idx[0] = i;
901           } else {
902             // there is a duplicated vertex in the face perimeter. The
903             // path might attach to either of the duplicate instances
904             // so we have to work out which is the right one to attach
905             // to. We assume it's the index currently being examined,
906             // if the path heads in a direction that's internal to the
907             // angle made by the prior and next edges of the face
908             // perimeter. Otherwise, leave it as the currently
909             // selected index (until another duplicate is found, if it
910             // exists, and is tested).
911             const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &p = *endpoint_indices[j].path;
912             const size_t pN = p.size();
913
914             carve::mesh::MeshSet<3>::vertex_t *a, *b, *c;
915             a = base_loop[(i+N-1)%N];
916             b = base_loop[i];
917             c = base_loop[(i+1)%N];
918
919             carve::mesh::MeshSet<3>::vertex_t *adj = (p[0] == base_loop[i]) ? p[1] : p[pN-2];
920
921             if (carve::geom2d::internalToAngle(face->project(c->v),
922                                                face->project(b->v),
923                                                face->project(a->v),
924                                                face->project(adj->v))) {
925               endpoint_indices[j].edge_idx[0] = i;
926             }
927           }
928         }
929
930         // test end of path.
931         if (paths[j].back() == base_loop[i]) {
932           if (endpoint_indices[j].edge_idx[1] == N) {
933             endpoint_indices[j].edge_idx[1] = i;
934           } else {
935             // Work out which of the duplicated vertices is the right
936             // one to attach to, as above.
937             const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &p = *endpoint_indices[j].path;
938             const size_t pN = p.size();
939
940             carve::mesh::MeshSet<3>::vertex_t *a, *b, *c;
941             a = base_loop[(i+N-1)%N];
942             b = base_loop[i];
943             c = base_loop[(i+1)%N];
944
945             carve::mesh::MeshSet<3>::vertex_t *adj = (p[0] == base_loop[i]) ? p[1] : p[pN-2];
946
947             if (carve::geom2d::internalToAngle(face->project(c->v),
948                                                face->project(b->v),
949                                                face->project(a->v),
950                                                face->project(adj->v))) {
951               endpoint_indices[j].edge_idx[1] = i;
952             }
953           }
954         }
955       }
956     }
957
958 #if defined(CARVE_DEBUG)
959     std::cerr << "### N: " << N << std::endl;
960     for (size_t i = 0; i < paths.size(); ++i) {
961       std::cerr << "### path: " << i << " endpoints: " << endpoint_indices[i].edge_idx[0] << " - " << endpoint_indices[i].edge_idx[1] << std::endl;
962     }
963 #endif
964
965
966     // Step 2:
967     // divide paths up into those that connect to the base loop in two
968     // places (cross), and those that do not (noncross).
969     std::vector<crossing_data> cross, noncross;
970     cross.reserve(endpoint_indices.size() + 1);
971     noncross.reserve(endpoint_indices.size());
972
973     for (size_t i = 0; i < endpoint_indices.size(); ++i) {
974 #if defined(CARVE_DEBUG)
975       std::cerr << "### orienting path: " << i << " endpoints: " << endpoint_indices[i].edge_idx[0] << " - " << endpoint_indices[i].edge_idx[1] << std::endl;
976 #endif
977       if (endpoint_indices[i].edge_idx[0] != N && endpoint_indices[i].edge_idx[1] != N) {
978         // Orient each path correctly. Paths should progress from
979         // smaller perimeter index to larger, but if the path starts
980         // and ends at the same perimeter index, then the decision
981         // needs to be made based upon area.
982         if (endpoint_indices[i].edge_idx[0] == endpoint_indices[i].edge_idx[1]) {
983           // The path forms a loop that starts and ends at the same
984           // vertex of the perimeter. In this case, we need to orient
985           // the path so that the constructed loop has the right
986           // signed area.
987           double area = carve::geom2d::signedArea(endpoint_indices[i].path->begin() + 1,
988                                                   endpoint_indices[i].path->end(),
989                                                   carve::mesh::MeshSet<3>::face_t::projection_mapping(face->project));
990           if (area < 0) {
991             // XXX: Create test case to check that this is the correct sign for the area.
992             std::reverse(endpoint_indices[i].path->begin(), endpoint_indices[i].path->end());
993           }
994         } else {
995           if (endpoint_indices[i].edge_idx[0] > endpoint_indices[i].edge_idx[1]) {
996             std::swap(endpoint_indices[i].edge_idx[0], endpoint_indices[i].edge_idx[1]);
997             std::reverse(endpoint_indices[i].path->begin(), endpoint_indices[i].path->end());
998           }
999         }
1000       }
1001
1002       if (endpoint_indices[i].edge_idx[0] != N &&
1003           endpoint_indices[i].edge_idx[1] != N &&
1004           endpoint_indices[i].edge_idx[0] != endpoint_indices[i].edge_idx[1]) {
1005         cross.push_back(endpoint_indices[i]);
1006       } else {
1007         noncross.push_back(endpoint_indices[i]);
1008       }
1009     }
1010
1011     // Step 3:
1012     // add a temporary crossing path that connects the beginning and the
1013     // end of the base loop. this stops us from needing special case
1014     // code to handle the left over loop after all the other crossing
1015     // paths are considered.
1016     std::vector<carve::mesh::MeshSet<3>::vertex_t *> base_loop_temp_path;
1017     base_loop_temp_path.reserve(2);
1018     base_loop_temp_path.push_back(base_loop.front());
1019     base_loop_temp_path.push_back(base_loop.back());
1020
1021     cross.push_back(crossing_data(&base_loop_temp_path, 0, base_loop.size() - 1));
1022 #if defined(CARVE_DEBUG)
1023     std::cerr << "### crossing edge count (with sentinel): " << cross.size() << std::endl;
1024 #endif
1025
1026     // Step 4:
1027     // sort paths by increasing beginning point and decreasing ending point.
1028     std::sort(cross.begin(), cross.end());
1029     std::sort(noncross.begin(), noncross.end());
1030
1031     // Step 5:
1032     // divide up the base loop based upon crossing paths.
1033     std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > divided_base_loop;
1034     divided_base_loop.reserve(cross.size());
1035     std::vector<carve::mesh::MeshSet<3>::vertex_t *> out;
1036
1037     for (size_t i = 0; i < cross.size(); ++i) {
1038       size_t j;
1039       for (j = i + 1;
1040            j < cross.size() &&
1041              cross[i].edge_idx[0] == cross[j].edge_idx[0] && 
1042              cross[i].edge_idx[1] == cross[j].edge_idx[1];
1043            ++j) {}
1044       if (j - i >= 2) {
1045         // when there are multiple paths that begin and end at the
1046         // same point, they need to be ordered so that the constructed
1047         // loops have the right orientation. this means that the loop
1048         // made by taking path(i+1) forward, then path(i) backward
1049         // needs to have negative area. this combined area is equal to
1050         // the area of path(i+1) minus the area of path(i). in turn
1051         // this means that the loop made by path path(i+1) alone has
1052         // to have smaller signed area than loop made by path(i).
1053         // thus, we sort paths in order of decreasing area.
1054
1055         std::vector<std::pair<double, std::vector<carve::mesh::MeshSet<3>::vertex_t *> *> > order;
1056         order.reserve(j - i);
1057         for (size_t k = i; k < j; ++k) {
1058           double area = carve::geom2d::signedArea(cross[k].path->begin(),
1059                                                   cross[k].path->end(),
1060                                                   carve::mesh::MeshSet<3>::face_t::projection_mapping(face->project));
1061 #if defined(CARVE_DEBUG)
1062           std::cerr << "### k=" << k << " area=" << area << std::endl;
1063 #endif
1064           order.push_back(std::make_pair(-area, cross[k].path));
1065         }
1066         std::sort(order.begin(), order.end());
1067         for (size_t k = i; k < j; ++k) {
1068           cross[k].path = order[k-i].second;
1069 #if defined(CARVE_DEBUG)
1070           std::cerr << "### post-sort k=" << k << " cross[k].path->size()=" << cross[k].path->size() << std::endl;
1071 #endif
1072         }
1073       }
1074     }
1075
1076     // Step 6:
1077     for (size_t i = 0; i < cross.size(); ++i) {
1078 #if defined(CARVE_DEBUG)
1079       std::cerr << "### i=" << i << " working on edge: " << cross[i].edge_idx[0] << " - " << cross[i].edge_idx[1] << std::endl;
1080 #endif
1081       size_t e1_0 = cross[i].edge_idx[0];
1082       size_t e1_1 = cross[i].edge_idx[1];
1083       std::vector<carve::mesh::MeshSet<3>::vertex_t *> &p1 = *cross[i].path;
1084 #if defined(CARVE_DEBUG)
1085       std::cerr << "###     path size = " << p1.size() << std::endl;
1086 #endif
1087
1088       out.clear();
1089
1090       if (i < cross.size() - 1 &&
1091           cross[i+1].edge_idx[1] <= cross[i].edge_idx[1]) {
1092 #if defined(CARVE_DEBUG)
1093         std::cerr << "###     complex case" << std::endl;
1094 #endif
1095         // complex case. crossing path with other crossing paths embedded within.
1096         size_t pos = e1_0;
1097
1098         size_t skip = i+1;
1099
1100         while (pos != e1_1) {
1101
1102           std::vector<carve::mesh::MeshSet<3>::vertex_t *> &p2 = *cross[skip].path;
1103           size_t e2_0 = cross[skip].edge_idx[0];
1104           size_t e2_1 = cross[skip].edge_idx[1];
1105
1106           // copy up to the beginning of the next path.
1107           std::copy(base_loop.begin() + pos, base_loop.begin() + e2_0, std::back_inserter(out));
1108
1109           CARVE_ASSERT(base_loop[e2_0] == p2[0]);
1110           // copy the next path in the right direction.
1111           std::copy(p2.begin(), p2.end() - 1, std::back_inserter(out));
1112
1113           // move to the position of the end of the path.
1114           pos = e2_1;
1115
1116           // advance to the next hit path.
1117           do {
1118             ++skip;
1119           } while(skip != cross.size() && cross[skip].edge_idx[0] < e2_1);
1120
1121           if (skip == cross.size()) break;
1122
1123           // if the next hit path is past the start point of the current path, we're done.
1124           if (cross[skip].edge_idx[0] >= e1_1) break;
1125         }
1126
1127         // copy up to the end of the path.
1128         if (pos < e1_1) {
1129             std::copy(base_loop.begin() + pos, base_loop.begin() + e1_1, std::back_inserter(out));
1130         }
1131
1132         CARVE_ASSERT(base_loop[e1_1] == p1.back());
1133         std::copy(p1.rbegin(), p1.rend() - 1, std::back_inserter(out));
1134       } else {
1135         size_t loop_size = (e1_1 - e1_0) + (p1.size() - 1);
1136         out.reserve(loop_size);
1137
1138         std::copy(base_loop.begin() + e1_0, base_loop.begin() + e1_1, std::back_inserter(out));
1139         std::copy(p1.rbegin(), p1.rend() - 1, std::back_inserter(out));
1140
1141         CARVE_ASSERT(out.size() == loop_size);
1142       }
1143       divided_base_loop.push_back(out);
1144
1145 #if defined(CARVE_DEBUG)
1146       {
1147         std::vector<carve::geom2d::P2> projected;
1148         projected.reserve(out.size());
1149         for (size_t n = 0; n < out.size(); ++n) {
1150           projected.push_back(face->project(out[n]->v));
1151         }
1152
1153         double A = carve::geom2d::signedArea(projected);
1154         std::cerr << "### out area=" << A << std::endl;
1155         CARVE_ASSERT(A <= 0);
1156       }
1157 #endif
1158     }
1159
1160     if (!noncross.size()) {
1161       // If there are no non-crossing paths then we're done.
1162       populateListFromVector(divided_base_loop, face_loops_out);
1163       return true;
1164     }
1165
1166     // for each divided base loop, work out which noncrossing paths and
1167     // loops are part of it. use the old algorithm to combine these into
1168     // the divided base loop. if none, the divided base loop is just
1169     // output.
1170     std::vector<std::vector<carve::geom2d::P2> > proj;
1171     std::vector<carve::geom::aabb<2> > proj_aabb;
1172     proj.resize(divided_base_loop.size());
1173     proj_aabb.resize(divided_base_loop.size());
1174
1175     // calculate an aabb for each divided base loop, to avoid expensive
1176     // point-in-poly tests.
1177     for (size_t i = 0; i < divided_base_loop.size(); ++i) {
1178       proj[i].reserve(divided_base_loop[i].size());
1179       for (size_t j = 0; j < divided_base_loop[i].size(); ++j) {
1180         proj[i].push_back(face->project(divided_base_loop[i][j]->v));
1181       }
1182       proj_aabb[i].fit(proj[i].begin(), proj[i].end());
1183     }
1184
1185     for (size_t i = 0; i < divided_base_loop.size(); ++i) {
1186       std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> *> inc;
1187       carve::geom2d::P2 test;
1188
1189       // for each noncrossing path, choose an endpoint that isn't on the
1190       // base loop as a test point.
1191       for (size_t j = 0; j < noncross.size(); ++j) {
1192         if (noncross[j].edge_idx[0] < N) {
1193           if (noncross[j].path->front() == base_loop[noncross[j].edge_idx[0]]) {
1194             // noncrossing paths may be loops that run from the edge, back to the same vertex.
1195             if (noncross[j].path->front() == noncross[j].path->back()) {
1196               CARVE_ASSERT(noncross[j].path->size() > 2);
1197               test = face->project((*noncross[j].path)[1]->v);
1198             } else {
1199               test = face->project(noncross[j].path->back()->v);
1200             }
1201           } else {
1202             test = face->project(noncross[j].path->front()->v);
1203           }
1204         } else {
1205           test = face->project(noncross[j].path->front()->v);
1206         }
1207
1208         if (proj_aabb[i].intersects(test) &&
1209             carve::geom2d::pointInPoly(proj[i], test).iclass != carve::POINT_OUT) {
1210           inc.push_back(noncross[j].path);
1211         }
1212       }
1213
1214 #if defined(CARVE_DEBUG)
1215       std::cerr << "### divided base loop:" << i << " inc.size()=" << inc.size() << std::endl;
1216       std::cerr << "### inc = [";
1217       for (size_t j = 0; j < inc.size(); ++j) {
1218         std::cerr << " " << inc[j];
1219       }
1220       std::cerr << " ]" << std::endl;
1221 #endif
1222
1223       if (inc.size()) {
1224         carve::csg::V2Set face_edges;
1225
1226         for (size_t j = 0; j < divided_base_loop[i].size() - 1; ++j) {
1227           face_edges.insert(std::make_pair(divided_base_loop[i][j],
1228                                            divided_base_loop[i][j+1]));
1229         }
1230
1231         face_edges.insert(std::make_pair(divided_base_loop[i].back(),
1232                                          divided_base_loop[i].front()));
1233
1234         for (size_t j = 0; j < inc.size(); ++j) {
1235           std::vector<carve::mesh::MeshSet<3>::vertex_t *> &path = *inc[j];
1236           for (size_t k = 0; k < path.size() - 1; ++k) {
1237             face_edges.insert(std::make_pair(path[k], path[k+1]));
1238             face_edges.insert(std::make_pair(path[k+1], path[k]));
1239           }
1240         }
1241
1242         std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > face_loops;
1243         std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > hole_loops;
1244
1245         splitFace(face, face_edges, face_loops, hole_loops, vertex_intersections);
1246
1247         if (hole_loops.size()) {
1248           mergeFacesAndHoles(face, face_loops, hole_loops, hooks);
1249         }
1250         std::copy(face_loops.begin(), face_loops.end(), std::back_inserter(face_loops_out));
1251       } else {
1252         face_loops_out.push_back(divided_base_loop[i]);
1253       }
1254     }
1255     return true;
1256   }
1257
1258
1259
1260   void composeEdgesIntoPaths(const carve::csg::V2Set &edges,
1261                              const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &extra_endpoints,
1262                              std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &paths,
1263                              std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &cuts,
1264                              std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &loops) {
1265     using namespace carve::csg;
1266
1267     detail::VVSMap vertex_graph;
1268     detail::VSet endpoints;
1269     detail::VSet cut_endpoints;
1270
1271     typedef std::vector<carve::mesh::MeshSet<3>::vertex_t *> vvec_t;
1272     vvec_t path;
1273
1274     std::list<vvec_t> path_list, cut_list, loop_list;
1275
1276     // build graph from edges.
1277     for (V2Set::const_iterator i = edges.begin(); i != edges.end(); ++i) {
1278 #if defined(CARVE_DEBUG)
1279       std::cerr << "###    edge: " << (*i).first << " - " << (*i).second << std::endl;
1280 #endif
1281       vertex_graph[(*i).first].insert((*i).second);
1282       vertex_graph[(*i).second].insert((*i).first);
1283     }
1284
1285     // find the endpoints in the graph.
1286     // every vertex with number of incident edges != 2 is an endpoint.
1287     for (detail::VVSMap::const_iterator i = vertex_graph.begin(); i != vertex_graph.end(); ++i) {
1288       if ((*i).second.size() != 2) {
1289 #if defined(CARVE_DEBUG)
1290         std::cerr << "###    endpoint: " << (*i).first << std::endl;
1291 #endif
1292         endpoints.insert((*i).first);
1293         if ((*i).second.size() == 1) {
1294           cut_endpoints.insert((*i).first);
1295         }
1296       }
1297     }
1298
1299     // every vertex on the perimeter of the face is also an endpoint.
1300     for (size_t i = 0; i < extra_endpoints.size(); ++i) {
1301       if (vertex_graph.find(extra_endpoints[i]) != vertex_graph.end()) {
1302 #if defined(CARVE_DEBUG)
1303         std::cerr << "###    extra endpoint: " << extra_endpoints[i] << std::endl;
1304 #endif
1305         endpoints.insert(extra_endpoints[i]);
1306         cut_endpoints.erase(extra_endpoints[i]);
1307       }
1308     }
1309
1310     while (endpoints.size()) {
1311       carve::mesh::MeshSet<3>::vertex_t *v = *endpoints.begin();
1312       detail::VVSMap::iterator p = vertex_graph.find(v);
1313       if (p == vertex_graph.end()) {
1314         endpoints.erase(endpoints.begin());
1315         continue;
1316       }
1317
1318       path.clear();
1319       path.push_back(v);
1320
1321       for (;;) {
1322         CARVE_ASSERT(p != vertex_graph.end());
1323
1324         // pick a connected vertex to move to.
1325         if ((*p).second.size() == 0) break;
1326
1327         carve::mesh::MeshSet<3>::vertex_t *n = *((*p).second.begin());
1328         detail::VVSMap::iterator q = vertex_graph.find(n);
1329
1330         // remove the link.
1331         (*p).second.erase(n);
1332         (*q).second.erase(v);
1333
1334         // move on.
1335         v = n;
1336         path.push_back(v);
1337
1338         if ((*p).second.size() == 0) vertex_graph.erase(p);
1339         if ((*q).second.size() == 0) {
1340           vertex_graph.erase(q);
1341           q = vertex_graph.end();
1342         }
1343
1344         p = q;
1345
1346         if (v == path[0] || p == vertex_graph.end() || endpoints.find(v) != endpoints.end()) break;
1347       }
1348       CARVE_ASSERT(endpoints.find(path.back()) != endpoints.end());
1349
1350       bool is_cut =
1351         cut_endpoints.find(path.front()) != cut_endpoints.end() &&
1352         cut_endpoints.find(path.back()) != cut_endpoints.end();
1353
1354       if (is_cut) {
1355         cut_list.push_back(vvec_t()); path.swap(cut_list.back());
1356       } else {
1357         path_list.push_back(vvec_t()); path.swap(path_list.back());
1358       }
1359     }
1360
1361     populateVectorFromList(path_list, paths);
1362     populateVectorFromList(cut_list, cuts);
1363
1364     // now only loops should remain in the graph.
1365     while (vertex_graph.size()) {
1366       detail::VVSMap::iterator p = vertex_graph.begin();
1367       carve::mesh::MeshSet<3>::vertex_t *v = (*p).first;
1368       CARVE_ASSERT((*p).second.size() == 2);
1369
1370       std::vector<carve::mesh::MeshSet<3>::vertex_t *> path;
1371       path.clear();
1372       path.push_back(v);
1373
1374       for (;;) {
1375         CARVE_ASSERT(p != vertex_graph.end());
1376         // pick a connected vertex to move to.
1377
1378         carve::mesh::MeshSet<3>::vertex_t *n = *((*p).second.begin());
1379         detail::VVSMap::iterator q = vertex_graph.find(n);
1380
1381         // remove the link.
1382         (*p).second.erase(n);
1383         (*q).second.erase(v);
1384
1385         // move on.
1386         v = n;
1387         path.push_back(v);
1388
1389         if ((*p).second.size() == 0) vertex_graph.erase(p);
1390         if ((*q).second.size() == 0) vertex_graph.erase(q);
1391
1392         p = q;
1393
1394         if (v == path[0]) break;
1395       }
1396
1397       loop_list.push_back(vvec_t()); path.swap(loop_list.back());
1398     }
1399  
1400     populateVectorFromList(loop_list, loops);
1401   }
1402
1403
1404
1405   template<typename T>
1406   std::string ptrstr(const T *ptr) {
1407     std::ostringstream s;
1408     s << ptr;
1409     return s.str().substr(1);
1410   }
1411
1412 #if 0
1413   void dumpAsGraph(carve::mesh::MeshSet<3>::face_t *face,
1414                    const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &base_loop,
1415                    const carve::csg::V2Set &face_edges,
1416                    const carve::csg::V2Set &split_edges) {
1417     std::map<carve::mesh::MeshSet<3>::vertex_t *, carve::geom2d::P2> proj;
1418
1419     for (size_t i = 0; i < base_loop.size(); ++i) {
1420       proj[base_loop[i]] = face->project(base_loop[i]->v);
1421     }
1422     for (carve::csg::V2Set::const_iterator i = split_edges.begin(); i != split_edges.end(); ++i) {
1423       proj[(*i).first] = face->project((*i).first->v);
1424       proj[(*i).second] = face->project((*i).second->v);
1425     }
1426
1427     {
1428       carve::geom2d::P2 lo, hi;
1429       std::map<carve::mesh::MeshSet<3>::vertex_t *, carve::geom2d::P2>::iterator i;
1430       i = proj.begin();
1431       lo = hi = (*i).second;
1432       for (; i != proj.end(); ++i) {
1433         lo.x = std::min(lo.x, (*i).second.x); lo.y = std::min(lo.y, (*i).second.y);
1434         hi.x = std::max(hi.x, (*i).second.x); hi.y = std::max(hi.y, (*i).second.y);
1435       }
1436       for (i = proj.begin(); i != proj.end(); ++i) {
1437         (*i).second.x = ((*i).second.x - lo.x) / (hi.x - lo.x) * 10;
1438         (*i).second.y = ((*i).second.y - lo.y) / (hi.y - lo.y) * 10;
1439       }
1440     }
1441
1442     std::cerr << "graph G {\nnode [shape=circle,style=filled,fixedsize=true,width=\".1\",height=\".1\"];\nedge [len=4]\n";
1443     for (std::map<carve::mesh::MeshSet<3>::vertex_t *, carve::geom2d::P2>::iterator i = proj.begin(); i != proj.end(); ++i) {
1444       std::cerr << "   " << ptrstr((*i).first) << " [pos=\"" << (*i).second.x << "," << (*i).second.y << "!\"];\n";
1445     }
1446     for (carve::csg::V2Set::const_iterator i = face_edges.begin(); i != face_edges.end(); ++i) {
1447       std::cerr << "   " << ptrstr((*i).first) << " -- " << ptrstr((*i).second) << ";\n";
1448     }
1449     for (carve::csg::V2Set::const_iterator i = split_edges.begin(); i != split_edges.end(); ++i) {
1450       std::cerr << "   " << ptrstr((*i).first) << " -- " << ptrstr((*i).second) << " [color=\"blue\"];\n";
1451     }
1452     std::cerr << "};\n";
1453   }
1454 #endif
1455
1456   void generateOneFaceLoop(carve::mesh::MeshSet<3>::face_t *face,
1457                            const carve::csg::detail::Data &data,
1458                            const carve::csg::VertexIntersections &vertex_intersections,
1459                            carve::csg::CSG::Hooks &hooks,
1460                            std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > &face_loops) {
1461     using namespace carve::csg;
1462
1463     std::vector<carve::mesh::MeshSet<3>::vertex_t *> base_loop;
1464     std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > hole_loops;
1465
1466     /*bool face_edge_intersected = */assembleBaseLoop(face, data, base_loop, hooks);
1467
1468     detail::FV2SMap::const_iterator fse_iter = data.face_split_edges.find(face);
1469
1470     face_loops.clear();
1471
1472     if (fse_iter == data.face_split_edges.end()) {
1473       // simple case: input face is output face (possibly with the
1474       // addition of vertices at intersections).
1475       face_loops.push_back(base_loop);
1476       return;
1477     }
1478
1479     // complex case: input face is split into multiple output faces.
1480     V2Set face_edges;
1481
1482     for (size_t j = 0, je = base_loop.size() - 1; j < je; ++j) {
1483       face_edges.insert(std::make_pair(base_loop[j], base_loop[j + 1]));
1484     }
1485     face_edges.insert(std::make_pair(base_loop.back(), base_loop[0]));
1486
1487     // collect the split edges (as long as they're not on the perimeter)
1488     const detail::FV2SMap::mapped_type &fse = ((*fse_iter).second);
1489
1490     // split_edges contains all of the edges created by intersections
1491     // that aren't part of the perimeter of the face.
1492     V2Set split_edges;
1493
1494     for (detail::FV2SMap::mapped_type::const_iterator
1495            j = fse.begin(), je =  fse.end();
1496          j != je;
1497          ++j) {
1498       carve::mesh::MeshSet<3>::vertex_t *v1 = ((*j).first), *v2 = ((*j).second);
1499
1500       if (face_edges.find(std::make_pair(v1, v2)) == face_edges.end() &&
1501           face_edges.find(std::make_pair(v2, v1)) == face_edges.end()) {
1502         // If the edge isn't part of the face perimeter, add it to
1503         // split_edges.
1504         split_edges.insert(ordered_edge(v1, v2));
1505       }
1506     }
1507
1508     // face is unsplit.
1509     if (!split_edges.size()) {
1510       face_loops.push_back(base_loop);
1511       return;
1512     }
1513
1514 #if defined(CARVE_DEBUG)
1515     dumpAsGraph(face, base_loop, face_edges, split_edges);
1516 #endif
1517
1518 #if 0
1519     // old face splitting method.
1520     for (V2Set::const_iterator i = split_edges.begin(); i != split_edges.end(); ++i) {
1521       face_edges.insert(std::make_pair((*i).first, (*i).second));
1522       face_edges.insert(std::make_pair((*i).second, (*i).first));
1523     }
1524     splitFace(face, face_edges, face_loops, hole_loops, vertex_intersections);
1525
1526     if (hole_loops.size()) {
1527       mergeFacesAndHoles(face, face_loops, hole_loops, hooks);
1528     }
1529     return;
1530 #endif
1531
1532 #if defined(CARVE_DEBUG)
1533     std::cerr << "### split_edges.size(): " << split_edges.size() << std::endl;
1534 #endif
1535     if (split_edges.size() == 1) {
1536       // handle the common case of a face that's split by a single edge.
1537       carve::mesh::MeshSet<3>::vertex_t *v1 = split_edges.begin()->first;
1538       carve::mesh::MeshSet<3>::vertex_t *v2 = split_edges.begin()->second;
1539
1540       std::vector<carve::mesh::MeshSet<3>::vertex_t *>::iterator vi1 = std::find(base_loop.begin(), base_loop.end(), v1);
1541       std::vector<carve::mesh::MeshSet<3>::vertex_t *>::iterator vi2 = std::find(base_loop.begin(), base_loop.end(), v2);
1542
1543       if (vi1 != base_loop.end() && vi2 != base_loop.end()) {
1544         // this is an inserted edge that connects two points on the base loop. nice and simple.
1545         if (vi2 < vi1) std::swap(vi1, vi2);
1546
1547         size_t loop1_size = vi2 - vi1 + 1;
1548         size_t loop2_size = base_loop.size() + 2 - loop1_size;
1549
1550         std::vector<carve::mesh::MeshSet<3>::vertex_t *> l1;
1551         std::vector<carve::mesh::MeshSet<3>::vertex_t *> l2;
1552
1553         l1.reserve(loop1_size);
1554         l2.reserve(loop2_size);
1555
1556         std::copy(vi1, vi2+1, std::back_inserter(l1));
1557         std::copy(vi2, base_loop.end(), std::back_inserter(l2));
1558         std::copy(base_loop.begin(), vi1+1, std::back_inserter(l2));
1559
1560         CARVE_ASSERT(l1.size() == loop1_size);
1561         CARVE_ASSERT(l2.size() == loop2_size);
1562
1563         face_loops.push_back(l1);
1564         face_loops.push_back(l2);
1565
1566         return;
1567       }
1568
1569       // Consider handling cases where one end of the edge touches the
1570       // perimeter, and where neither end does.
1571     }
1572
1573     std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > paths;
1574     std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > cuts;
1575     std::vector<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > loops;
1576
1577     // Take the split edges and compose them into a set of paths and
1578     // loops. Loops are edge paths that do not touch the boundary, or
1579     // any other path or loop - they are holes cut out of the centre
1580     // of the face. Paths are made up of all the other edge segments,
1581     // and start and end at the face perimeter, or where they meet
1582     // another path (sometimes both cases will be true).
1583     composeEdgesIntoPaths(split_edges, base_loop, paths, cuts, loops);
1584
1585 #if defined(CARVE_DEBUG)
1586     std::cerr << "###   paths.size(): " << paths.size() << std::endl;
1587     std::cerr << "###   cuts.size():  " << cuts.size() << std::endl;
1588     std::cerr << "###   loops.size(): " << loops.size() << std::endl;
1589 #endif
1590
1591     if (!paths.size()) {
1592       // No complex paths.
1593       face_loops.push_back(base_loop);
1594     } else {
1595       if (processCrossingEdges(face, vertex_intersections, hooks, base_loop, paths, face_loops)) {
1596         // Worked.
1597       } else {
1598         // complex case - fall back to old edge tracing code.
1599 #if defined(CARVE_DEBUG)
1600         std::cerr << "### processCrossingEdges failed. Falling back to edge tracing code" << std::endl;
1601 #endif
1602         for (size_t i = 0; i < paths.size(); ++i) {
1603           for (size_t j = 0; j < paths[i].size() - 1; ++j) {
1604             face_edges.insert(std::make_pair(paths[i][j], paths[i][j+1]));
1605             face_edges.insert(std::make_pair(paths[i][j+1], paths[i][j]));
1606           }
1607         }
1608         splitFace(face, face_edges, face_loops, hole_loops, vertex_intersections);
1609       }
1610     }
1611
1612     // Now merge cuts and loops into face loops.
1613
1614     // every cut creates a hole.
1615     for (size_t i = 0; i < cuts.size(); ++i) {
1616       hole_loops.push_back(std::vector<carve::mesh::MeshSet<3>::vertex_t *>());
1617       hole_loops.back().reserve(2 * cuts[i].size() - 2);
1618       std::copy(cuts[i].begin(), cuts[i].end(), std::back_inserter(hole_loops.back()));
1619       if (cuts[i].size() > 2) {
1620         std::copy(cuts[i].rbegin() + 1, cuts[i].rend() - 1, std::back_inserter(hole_loops.back()));
1621       }
1622     }
1623
1624     // every loop creates a hole and a corresponding face.
1625     for (size_t i = 0; i < loops.size(); ++i) {
1626       hole_loops.push_back(std::vector<carve::mesh::MeshSet<3>::vertex_t *>());
1627       hole_loops.back().reserve(loops[i].size()-1);
1628       std::copy(loops[i].begin(), loops[i].end()-1, std::back_inserter(hole_loops.back()));
1629
1630       face_loops.push_back(std::vector<carve::mesh::MeshSet<3>::vertex_t *>());
1631       face_loops.back().reserve(loops[i].size()-1);
1632       std::copy(loops[i].rbegin()+1, loops[i].rend(), std::back_inserter(face_loops.back()));
1633
1634       std::vector<carve::geom2d::P2> projected;
1635       projected.reserve(face_loops.back().size());
1636       for (size_t i = 0; i < face_loops.back().size(); ++i) {
1637         projected.push_back(face->project(face_loops.back()[i]->v));
1638       }
1639
1640       if (carve::geom2d::signedArea(projected) > 0.0) {
1641         std::swap(face_loops.back(), hole_loops.back());
1642       }
1643     }
1644
1645     // if there are holes, then they need to be merged with faces.
1646     if (hole_loops.size()) {
1647       mergeFacesAndHoles(face, face_loops, hole_loops, hooks);
1648     }
1649   }
1650 }
1651
1652
1653
1654 /** 
1655  * \brief Build a set of face loops for all (split) faces of a Polyhedron.
1656  * 
1657  * @param[in] poly The polyhedron to process
1658  * @param[in] data Internal intersection data
1659  * @param[out] face_loops_out The resulting face loops
1660  * 
1661  * @return The number of edges generated.
1662  */
1663 size_t carve::csg::CSG::generateFaceLoops(carve::mesh::MeshSet<3> *poly,
1664                                           const detail::Data &data,
1665                                           FaceLoopList &face_loops_out) {
1666   static carve::TimingName FUNC_NAME("CSG::generateFaceLoops()");
1667   carve::TimingBlock block(FUNC_NAME);
1668   size_t generated_edges = 0;
1669   std::vector<carve::mesh::MeshSet<3>::vertex_t *> base_loop;
1670   std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> > face_loops;
1671   
1672   for (carve::mesh::MeshSet<3>::face_iter i = poly->faceBegin(); i != poly->faceEnd(); ++i) {
1673     carve::mesh::MeshSet<3>::face_t *face = (*i);
1674
1675 #if defined(CARVE_DEBUG)
1676     double in_area = 0.0, out_area = 0.0;
1677
1678     {
1679       std::vector<carve::mesh::MeshSet<3>::vertex_t *> base_loop;
1680       assembleBaseLoop(face, data, base_loop);
1681
1682       {
1683         std::vector<carve::geom2d::P2> projected;
1684         projected.reserve(base_loop.size());
1685         for (size_t n = 0; n < base_loop.size(); ++n) {
1686           projected.push_back(face->project(base_loop[n]->v));
1687         }
1688
1689         in_area = carve::geom2d::signedArea(projected);
1690         std::cerr << "### in_area=" << in_area << std::endl;
1691       }
1692     }
1693 #endif
1694
1695     generateOneFaceLoop(face, data, vertex_intersections, hooks, face_loops);
1696
1697 #if defined(CARVE_DEBUG)
1698     {
1699       V2Set face_edges;
1700
1701       std::vector<carve::mesh::MeshSet<3>::vertex_t *> base_loop;
1702       assembleBaseLoop(face, data, base_loop);
1703
1704       for (size_t j = 0, je = base_loop.size() - 1; j < je; ++j) {
1705         face_edges.insert(std::make_pair(base_loop[j+1], base_loop[j]));
1706       }
1707       face_edges.insert(std::make_pair(base_loop[0], base_loop.back()));
1708       for (std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> >::const_iterator fli = face_loops.begin(); fli != face_loops.end(); ++ fli) {
1709
1710         {
1711           std::vector<carve::geom2d::P2> projected;
1712           projected.reserve((*fli).size());
1713           for (size_t n = 0; n < (*fli).size(); ++n) {
1714             projected.push_back(face->project((*fli)[n]->v));
1715           }
1716
1717           double area = carve::geom2d::signedArea(projected);
1718           std::cerr << "### loop_area[" << std::distance((std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> >::const_iterator)face_loops.begin(), fli) << "]=" << area << std::endl;
1719           out_area += area;
1720         }
1721
1722         const std::vector<carve::mesh::MeshSet<3>::vertex_t *> &fl = *fli;
1723         for (size_t j = 0, je = fl.size() - 1; j < je; ++j) {
1724           face_edges.insert(std::make_pair(fl[j], fl[j+1]));
1725         }
1726         face_edges.insert(std::make_pair(fl.back(), fl[0]));
1727       }
1728       for (V2Set::const_iterator j = face_edges.begin(); j != face_edges.end(); ++j) {
1729         if (face_edges.find(std::make_pair((*j).second, (*j).first)) == face_edges.end()) {
1730           std::cerr << "### error: unmatched edge [" << (*j).first << "-" << (*j).second << "]" << std::endl;
1731         }
1732       }
1733       std::cerr << "### out_area=" << out_area << std::endl;
1734       if (out_area != in_area) {
1735         std::cerr << "### error: area does not match. delta = " << (out_area - in_area) << std::endl;
1736         // CARVE_ASSERT(fabs(out_area - in_area) < 1e-5);
1737       }
1738     }
1739 #endif
1740
1741     // now record all the resulting face loops.
1742 #if defined(CARVE_DEBUG)
1743     std::cerr << "### ======" << std::endl;
1744 #endif
1745     for (std::list<std::vector<carve::mesh::MeshSet<3>::vertex_t *> >::const_iterator
1746            f = face_loops.begin(), fe = face_loops.end();
1747          f != fe;
1748          ++f) {
1749 #if defined(CARVE_DEBUG)
1750       std::cerr << "### loop:";
1751       for (size_t i = 0; i < (*f).size(); ++i) {
1752         std::cerr << " " << (*f)[i];
1753       }
1754       std::cerr << std::endl;
1755 #endif
1756
1757       face_loops_out.append(new FaceLoop(face, *f));
1758       generated_edges += (*f).size();
1759     }
1760 #if defined(CARVE_DEBUG)
1761     std::cerr << "### ======" << std::endl;
1762 #endif
1763   }
1764   return generated_edges;
1765 }