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