Silence some annoying warnings when doing full build with strict flags
[blender.git] / extern / carve / lib / mesh.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/mesh.hpp>
23 #include <carve/mesh_impl.hpp>
24 #include <carve/rtree.hpp>
25
26 #include <carve/poly.hpp>
27
28 namespace {
29   inline double CALC_X(const carve::geom::plane<3> &p, double y, double z) { return -(p.d + p.N.y * y + p.N.z * z) / p.N.x; }
30   inline double CALC_Y(const carve::geom::plane<3> &p, double x, double z) { return -(p.d + p.N.x * x + p.N.z * z) / p.N.y; }
31   inline double CALC_Z(const carve::geom::plane<3> &p, double x, double y) { return -(p.d + p.N.x * x + p.N.y * y) / p.N.z; }
32
33   carve::geom::vector<2> _project_1(const carve::geom::vector<3> &v) {
34     return carve::geom::VECTOR(v.z, v.y);
35   }
36
37   carve::geom::vector<2> _project_2(const carve::geom::vector<3> &v) {
38     return carve::geom::VECTOR(v.x, v.z);
39   }
40
41   carve::geom::vector<2> _project_3(const carve::geom::vector<3> &v) {
42     return carve::geom::VECTOR(v.y, v.x);
43   }
44
45   carve::geom::vector<2> _project_4(const carve::geom::vector<3> &v) {
46     return carve::geom::VECTOR(v.y, v.z);
47   }
48
49   carve::geom::vector<2> _project_5(const carve::geom::vector<3> &v) {
50     return carve::geom::VECTOR(v.z, v.x);
51   }
52
53   carve::geom::vector<2> _project_6(const carve::geom::vector<3> &v) {
54     return carve::geom::VECTOR(v.x, v.y);
55   }
56   
57   carve::geom::vector<3> _unproject_1(const carve::geom::vector<2> &p, const carve::geom3d::Plane &plane) {
58     return carve::geom::VECTOR(CALC_X(plane, p.y, p.x), p.y, p.x);
59   }
60
61   carve::geom::vector<3> _unproject_2(const carve::geom::vector<2> &p, const carve::geom3d::Plane &plane) {
62     return carve::geom::VECTOR(p.x, CALC_Y(plane, p.x, p.y), p.y);
63   }
64
65   carve::geom::vector<3> _unproject_3(const carve::geom::vector<2> &p, const carve::geom3d::Plane &plane) {
66     return carve::geom::VECTOR(p.y, p.x, CALC_Z(plane, p.y, p.x));
67   }
68
69   carve::geom::vector<3> _unproject_4(const carve::geom::vector<2> &p, const carve::geom3d::Plane &plane) {
70     return carve::geom::VECTOR(CALC_X(plane, p.x, p.y), p.x, p.y);
71   }
72
73   carve::geom::vector<3> _unproject_5(const carve::geom::vector<2> &p, const carve::geom3d::Plane &plane) {
74     return carve::geom::VECTOR(p.y, CALC_Y(plane, p.y, p.x), p.x);
75   }
76
77   carve::geom::vector<3> _unproject_6(const carve::geom::vector<2> &p, const carve::geom3d::Plane &plane) {
78     return carve::geom::VECTOR(p.x, p.y, CALC_Z(plane, p.x, p.y));
79   }
80
81   static carve::geom::vector<2> (*project_tab[2][3])(const carve::geom::vector<3> &) = {
82     { &_project_1, &_project_2, &_project_3 },
83     { &_project_4, &_project_5, &_project_6 }
84   };
85
86   static carve::geom::vector<3> (*unproject_tab[2][3])(const carve::geom::vector<2> &, const carve::geom3d::Plane &) = {
87     { &_unproject_1, &_unproject_2, &_unproject_3 },
88     { &_unproject_4, &_unproject_5, &_unproject_6 }
89   };
90
91 }
92
93 namespace carve {
94   namespace mesh {
95
96
97
98     template<unsigned ndim>
99     typename Face<ndim>::project_t Face<ndim>::getProjector(bool positive_facing, int axis) const {
100       return NULL;
101     }
102
103
104
105     template<>
106     Face<3>::project_t Face<3>::getProjector(bool positive_facing, int axis) const {
107       return project_tab[positive_facing ? 1 : 0][axis];
108     }
109
110
111
112     template<unsigned ndim>
113     typename Face<ndim>::unproject_t Face<ndim>::getUnprojector(bool positive_facing, int axis) const {
114       return NULL;
115     }
116
117
118
119     template<>
120     Face<3>::unproject_t Face<3>::getUnprojector(bool positive_facing, int axis) const {
121       return unproject_tab[positive_facing ? 1 : 0][axis];
122     }
123
124
125
126     template<unsigned ndim>
127     bool Face<ndim>::containsPoint(const vector_t &p) const {
128       if (!carve::math::ZERO(carve::geom::distance(plane, p))) return false;
129       // return pointInPolySimple(vertices, projector(), (this->*project)(p));
130       std::vector<carve::geom::vector<2> > verts;
131       getProjectedVertices(verts);
132       return carve::geom2d::pointInPoly(verts, project(p)).iclass != carve::POINT_OUT;
133     }
134
135
136
137     template<unsigned ndim>
138     bool Face<ndim>::containsPointInProjection(const vector_t &p) const {
139       std::vector<carve::geom::vector<2> > verts;
140       getProjectedVertices(verts);
141       return carve::geom2d::pointInPoly(verts, project(p)).iclass != carve::POINT_OUT;
142     }
143
144
145
146     template<unsigned ndim>
147     bool Face<ndim>::simpleLineSegmentIntersection(
148         const carve::geom::linesegment<ndim> &line,
149         vector_t &intersection) const {
150       if (!line.OK()) return false;
151
152       carve::mesh::MeshSet<3>::vertex_t::vector_t p;
153       carve::IntersectionClass intersects =
154         carve::geom3d::lineSegmentPlaneIntersection(plane, line, p);
155       if (intersects == carve::INTERSECT_NONE || intersects == carve::INTERSECT_BAD) {
156         return false;
157       }
158
159       std::vector<carve::geom::vector<2> > verts;
160       getProjectedVertices(verts);
161       if (carve::geom2d::pointInPolySimple(verts, project(p))) {
162         intersection = p;
163         return true;
164       }
165       return false;
166     }
167
168
169
170     template<unsigned ndim>
171     IntersectionClass Face<ndim>::lineSegmentIntersection(const carve::geom::linesegment<ndim> &line,
172                                                           vector_t &intersection) const {
173       if (!line.OK()) return INTERSECT_NONE;
174
175   
176       vector_t p;
177       IntersectionClass intersects = carve::geom3d::lineSegmentPlaneIntersection(plane, line, p);
178       if (intersects == INTERSECT_NONE || intersects == INTERSECT_BAD) {
179         return intersects;
180       }
181
182       std::vector<carve::geom::vector<2> > verts;
183       getProjectedVertices(verts);
184       carve::geom2d::PolyInclusionInfo pi = carve::geom2d::pointInPoly(verts, project(p));
185       switch (pi.iclass) {
186       case POINT_VERTEX:
187         intersection = p;
188         return INTERSECT_VERTEX;
189
190       case POINT_EDGE:
191         intersection = p;
192         return INTERSECT_EDGE;
193
194       case POINT_IN:
195         intersection = p;
196         return INTERSECT_FACE;
197       
198       case POINT_OUT:
199         return INTERSECT_NONE;
200
201       default:
202         break;
203       }
204       return INTERSECT_NONE;
205     }
206
207
208
209     template<unsigned ndim>
210     Face<ndim> *Face<ndim>::closeLoop(typename Face<ndim>::edge_t *start) {
211       edge_t *e = start;
212       std::vector<edge_t *> loop_edges;
213       do {
214         CARVE_ASSERT(e->rev == NULL);
215         loop_edges.push_back(e);
216         e = e->perimNext();
217       } while (e != start);
218
219       const size_t N = loop_edges.size();
220       for (size_t i = 0; i < N; ++i) {
221         loop_edges[i]->rev = new edge_t(loop_edges[i]->v2(), NULL);
222       }
223
224       for (size_t i = 0; i < N; ++i) {
225         edge_t *e1 = loop_edges[i]->rev;
226         edge_t *e2 = loop_edges[(i+1)%N]->rev;
227         e1->prev = e2;
228         e2->next = e1;
229       }
230
231       Face *f = new Face(start->rev);
232
233       CARVE_ASSERT(f->n_edges == N);
234
235       return f;
236     }
237
238
239
240     namespace detail {
241
242
243
244       bool FaceStitcher::EdgeOrderData::Cmp::operator()(const EdgeOrderData &a, const EdgeOrderData &b) const {
245         int v = carve::geom3d::compareAngles(edge_dir, base_dir, a.face_dir, b.face_dir);
246
247 #if defined(CARVE_DEBUG)
248         {
249           double da = carve::geom3d::antiClockwiseAngle(base_dir, a.face_dir, edge_dir);
250           double db = carve::geom3d::antiClockwiseAngle(base_dir, b.face_dir, edge_dir);
251           int v_cmp = 0;
252           if (da < db) v_cmp = -1;
253           if (db < da) v_cmp = +1;
254           if (v_cmp != v) {
255             std::cerr << "v= " << v << " v_cmp= " << v_cmp << " da= " << da << " db= " << db << "  edge_dir=" << edge_dir << " base_dir=" << base_dir << " a=" << a.face_dir << " b=" << b.face_dir << std::endl;
256           }
257         }
258 #endif
259
260         if (v < 0) return true;
261         if (v == 0) {
262           if (a.is_reversed && !b.is_reversed) return true;
263           if (a.is_reversed == b.is_reversed) {
264             return a.group_id < b.group_id;
265           }
266         }
267         return false;
268       }
269
270
271
272       void FaceStitcher::matchSimpleEdges() {
273         // join faces that share an edge, where no other faces are incident.
274         for (edge_map_t::iterator i = edges.begin(); i != edges.end(); ++i) {
275           const vpair_t &ev = (*i).first;
276           edge_map_t::iterator j = edges.find(vpair_t(ev.second, ev.first));
277           if (j == edges.end()) {
278             for (edgelist_t::iterator k = (*i).second.begin(); k != (*i).second.end(); ++k) {
279               is_open[ (*k)->face->id] = true;
280             }
281           } else if ((*i).second.size() != 1 || (*j).second.size() != 1) {
282             std::swap(complex_edges[(*i).first], (*i).second);
283           } else {
284             // simple edge.
285             edge_t *a = (*i).second.front();
286             edge_t *b = (*j).second.front();
287             if (a < b) {
288               // every simple edge pair is encountered twice. only merge once.
289               a->rev = b;
290               b->rev = a;
291               face_groups.merge_sets(a->face->id, b->face->id);
292             }
293           }
294         }
295       }
296
297
298
299       size_t FaceStitcher::faceGroupID(const Face<3> *face) {
300         return face_groups.find_set_head(face->id);
301       }
302
303
304
305       size_t FaceStitcher::faceGroupID(const Edge<3> *edge) {
306         return face_groups.find_set_head(edge->face->id);
307       }
308
309
310
311       void FaceStitcher::orderForwardAndReverseEdges(std::vector<std::vector<Edge<3> *> > &efwd,
312                                                      std::vector<std::vector<Edge<3> *> > &erev,
313                                                      std::vector<std::vector<EdgeOrderData> > &result) {
314         const size_t Nfwd = efwd.size();
315         const size_t Nrev = erev.size();
316         const size_t N = efwd[0].size();
317
318         result.resize(N);
319
320         for (size_t i = 0; i < N; ++i) {
321           Edge<3> *base = efwd[0][i];
322
323           result[i].reserve(Nfwd + Nrev);
324           for (size_t j = 0; j < Nfwd; ++j) {
325             result[i].push_back(EdgeOrderData(efwd[j][i], j, false));
326             CARVE_ASSERT(efwd[0][i]->v1() == efwd[j][i]->v1());
327             CARVE_ASSERT(efwd[0][i]->v2() == efwd[j][i]->v2());
328           }
329           for (size_t j = 0; j < Nrev; ++j) {
330             result[i].push_back(EdgeOrderData(erev[j][i], j, true));
331             CARVE_ASSERT(erev[0][i]->v1() == erev[j][i]->v1());
332             CARVE_ASSERT(erev[0][i]->v2() == erev[j][i]->v2());
333           }
334
335           geom::vector<3> sort_dir;
336           if (opts.opt_avoid_cavities) {
337             sort_dir = base->v1()->v - base->v2()->v;
338           } else {
339             sort_dir = base->v2()->v - base->v1()->v;
340           }
341
342           std::sort(result[i].begin(), result[i].end(), EdgeOrderData::Cmp(sort_dir, result[i][0].face_dir));
343         }
344       }
345
346
347
348       void FaceStitcher::edgeIncidentGroups(const vpair_t &e,
349                                             const edge_map_t &all_edges,
350                                             std::pair<std::set<size_t>, std::set<size_t> > &groups) {
351         groups.first.clear();
352         groups.second.clear();
353         edge_map_t::const_iterator i;
354
355         i = all_edges.find(e);
356         if (i != all_edges.end()) {
357           for (edgelist_t::const_iterator j = (*i).second.begin(); j != (*i).second.end(); ++j) {
358             groups.first.insert(faceGroupID(*j));
359           }
360         }
361
362         i = all_edges.find(vpair_t(e.second, e.first));
363         if (i != all_edges.end()) {
364           for (edgelist_t::const_iterator j = (*i).second.begin(); j != (*i).second.end(); ++j) {
365             groups.second.insert(faceGroupID(*j));
366           }
367         }
368       }
369
370
371
372       void FaceStitcher::buildEdgeGraph(const edge_map_t &all_edges) {
373         for (edge_map_t::const_iterator i = all_edges.begin();
374              i != all_edges.end();
375              ++i) {
376           edge_graph[(*i).first.first].insert((*i).first.second);
377         }
378       }
379
380
381
382       void FaceStitcher::extractPath(std::vector<const vertex_t *> &path) {
383         path.clear();
384
385         edge_graph_t::iterator iter = edge_graph.begin();
386
387         
388         const vertex_t *init = (*iter).first;
389         const vertex_t *next = *(*iter).second.begin();
390         const vertex_t *prev = NULL;
391         const vertex_t *vert = init;
392
393         while ((*iter).second.size() == 2) {
394           prev = *std::find_if((*iter).second.begin(),
395                                (*iter).second.end(),
396                                std::bind2nd(std::not_equal_to<const vertex_t *>(), next));
397           next = vert;
398           vert = prev;
399           iter = edge_graph.find(vert);
400           CARVE_ASSERT(iter != edge_graph.end());
401           if (vert == init) break;
402         }
403         init = vert;
404
405         std::vector<const edge_t *> efwd;
406         std::vector<const edge_t *> erev;
407
408         edge_map_t::iterator edgeiter;
409         edgeiter = complex_edges.find(vpair_t(vert, next));
410         std::copy((*edgeiter).second.begin(), (*edgeiter).second.end(), std::back_inserter(efwd));
411
412         edgeiter = complex_edges.find(vpair_t(next, vert));
413         std::copy((*edgeiter).second.begin(), (*edgeiter).second.end(), std::back_inserter(erev));
414
415         path.push_back(vert);
416
417         prev = vert;
418         vert = next;
419         path.push_back(vert);
420         iter = edge_graph.find(vert);
421         CARVE_ASSERT(iter != edge_graph.end());
422
423         while (vert != init && (*iter).second.size() == 2) {
424           next = *std::find_if((*iter).second.begin(),
425                                (*iter).second.end(),
426                                std::bind2nd(std::not_equal_to<const vertex_t *>(), prev));
427
428           edgeiter = complex_edges.find(vpair_t(vert, next));
429           if ((*edgeiter).second.size() != efwd.size()) goto done;
430
431           for (size_t i = 0; i < efwd.size(); ++i) {
432             Edge<3> *e_next = efwd[i]->perimNext();
433             if (e_next->v2() != next) goto done;
434             efwd[i] = e_next;
435           }
436
437           edgeiter = complex_edges.find(vpair_t(next, vert));
438           if ((*edgeiter).second.size() != erev.size()) goto done;
439
440           for (size_t i = 0; i < erev.size(); ++i) {
441             Edge<3> *e_prev = erev[i]->perimPrev();
442             if (e_prev->v1() != next) goto done;
443             erev[i] = e_prev;
444           }
445
446           prev = vert;
447           vert = next;
448           path.push_back(vert);
449           iter = edge_graph.find(vert);
450           CARVE_ASSERT(iter != edge_graph.end());
451         }
452       done:;
453       }
454
455
456
457       void FaceStitcher::removePath(const std::vector<const vertex_t *> &path) {
458         for (size_t i = 1; i < path.size() - 1; ++i) {
459           edge_graph.erase(path[i]);
460         }
461
462         edge_graph[path[0]].erase(path[1]);
463         if (edge_graph[path[0]].size() == 0) {
464           edge_graph.erase(path[0]);
465         }
466
467         edge_graph[path[path.size()-1]].erase(path[path.size()-2]);
468         if (edge_graph[path[path.size()-1]].size() == 0) {
469           edge_graph.erase(path[path.size()-1]);
470         }
471       }
472
473
474
475       void FaceStitcher::reorder(std::vector<EdgeOrderData> &ordering,
476                                  size_t grp) {
477         if (!ordering[0].is_reversed && ordering[0].group_id == grp) return;
478         for (size_t i = 1; i < ordering.size(); ++i) {
479           if (!ordering[i].is_reversed && ordering[i].group_id == grp) {
480             std::vector<EdgeOrderData> temp;
481             temp.reserve(ordering.size());
482             std::copy(ordering.begin() + i, ordering.end(), std::back_inserter(temp));
483             std::copy(ordering.begin(), ordering.begin() + i, std::back_inserter(temp));
484             std::copy(temp.begin(), temp.end(), ordering.begin());
485             return;
486           }
487         }
488       }
489
490
491
492       struct lt_second {
493         template<typename pair_t>
494         bool operator()(const pair_t &a, const pair_t &b) const {
495           return a.second < b.second;
496         }
497       };
498
499
500
501       void FaceStitcher::fuseEdges(std::vector<Edge<3> *> &fwd,
502                                    std::vector<Edge<3> *> &rev) {
503         for (size_t i = 0; i < fwd.size(); ++i) {
504           fwd[i]->rev = rev[i];
505           rev[i]->rev = fwd[i];
506           face_groups.merge_sets(fwd[i]->face->id, rev[i]->face->id);
507         }
508       }
509
510
511
512       void FaceStitcher::joinGroups(std::vector<std::vector<Edge<3> *> > &efwd,
513                                     std::vector<std::vector<Edge<3> *> > &erev,
514                                     size_t fwd_grp,
515                                     size_t rev_grp) {
516         fuseEdges(efwd[fwd_grp], erev[rev_grp]);
517       }
518
519
520
521       void FaceStitcher::matchOrderedEdges(const std::vector<std::vector<EdgeOrderData> >::iterator begin,
522                                            const std::vector<std::vector<EdgeOrderData> >::iterator end,
523                                            std::vector<std::vector<Edge<3> *> > &efwd,
524                                            std::vector<std::vector<Edge<3> *> > &erev) {
525         typedef std::unordered_map<std::pair<size_t, size_t>, size_t> pair_counts_t;
526         for (;;) {
527           pair_counts_t pair_counts;
528
529           for (std::vector<std::vector<EdgeOrderData> >::iterator i = begin; i != end; ++i) {
530             std::vector<EdgeOrderData> &e = *i;
531             for (size_t j = 0; j < e.size(); ++j) {
532               if (!e[j].is_reversed && e[(j+1)%e.size()].is_reversed) {
533                 pair_counts[std::make_pair(e[j].group_id,
534                                            e[(j+1)%e.size()].group_id)]++;
535               }
536             }
537           }
538
539           if (!pair_counts.size()) break;
540
541           std::vector<std::pair<size_t, std::pair<size_t, size_t> > > counts;
542           counts.reserve(pair_counts.size());
543           for (pair_counts_t::iterator iter = pair_counts.begin(); iter != pair_counts.end(); ++iter) {
544             counts.push_back(std::make_pair((*iter).second, (*iter).first));
545           }
546           std::make_heap(counts.begin(), counts.end());
547
548           std::set<size_t> rem_fwd, rem_rev;
549
550           while (counts.size()) {
551             std::pair<size_t, size_t> join = counts.front().second;
552             std::pop_heap(counts.begin(), counts.end());
553             counts.pop_back();
554             if (rem_fwd.find(join.first) != rem_fwd.end()) continue;
555             if (rem_rev.find(join.second) != rem_rev.end()) continue;
556
557             size_t g1 = join.first;
558             size_t g2 = join.second;
559
560             joinGroups(efwd, erev, g1, g2);
561
562             for (std::vector<std::vector<EdgeOrderData> >::iterator i = begin; i != end; ++i) {
563               (*i).erase(std::remove_if((*i).begin(), (*i).end(), EdgeOrderData::TestGroups(g1, g2)), (*i).end());
564             }
565
566             rem_fwd.insert(g1);
567             rem_rev.insert(g2);
568           }
569         }
570       }
571
572
573
574       void FaceStitcher::resolveOpenEdges() {
575         // Remove open regions of mesh. Doing this may make additional
576         // edges simple (for example, removing a fin from the edge of
577         // a cube), and may also expose more open mesh regions. In the
578         // latter case, the process must be repeated to deal with the
579         // newly uncovered regions.
580         std::unordered_set<size_t> open_groups;
581
582         for (size_t i = 0; i < is_open.size(); ++i) {
583           if (is_open[i]) open_groups.insert(face_groups.find_set_head(i));
584         }
585
586         while (!open_groups.empty()) {
587           std::list<vpair_t> edge_0, edge_1;
588
589           for (edge_map_t::iterator i = complex_edges.begin(); i != complex_edges.end(); ++i) {
590             bool was_modified = false;
591             for(edgelist_t::iterator j = (*i).second.begin(); j != (*i).second.end(); ) {
592               if (open_groups.find(faceGroupID(*j)) != open_groups.end()) {
593                 j = (*i).second.erase(j);
594                 was_modified = true;
595               } else {
596                 ++j;
597               }
598             }
599             if (was_modified) {
600               if ((*i).second.empty()) {
601                 edge_0.push_back((*i).first);
602               } else if ((*i).second.size() == 1) {
603                 edge_1.push_back((*i).first);
604               }
605             }
606           }
607
608           for (std::list<vpair_t>::iterator i = edge_1.begin(); i != edge_1.end(); ++i) {
609             vpair_t e1 = *i;
610             edge_map_t::iterator e1i = complex_edges.find(e1);
611             if (e1i == complex_edges.end()) continue;
612             vpair_t e2 = vpair_t(e1.second, e1.first);
613             edge_map_t::iterator e2i = complex_edges.find(e2);
614             CARVE_ASSERT(e2i != complex_edges.end()); // each complex edge should have a mate.
615
616             if ((*e2i).second.size() == 1) {
617               // merge newly simple edges, delete both from complex_edges.
618               edge_t *a = (*e1i).second.front();
619               edge_t *b = (*e2i).second.front();
620               a->rev = b;
621               b->rev = a;
622               face_groups.merge_sets(a->face->id, b->face->id);
623               complex_edges.erase(e1i);
624               complex_edges.erase(e2i);
625             }
626           }
627
628           open_groups.clear();
629
630           for (std::list<vpair_t>::iterator i = edge_0.begin(); i != edge_0.end(); ++i) {
631             vpair_t e1 = *i;
632             edge_map_t::iterator e1i = complex_edges.find(e1);
633             vpair_t e2 = vpair_t(e1.second, e1.first);
634             edge_map_t::iterator e2i = complex_edges.find(e2);
635             if (e2i == complex_edges.end()) {
636               // This could occur, for example, when two faces share
637               // an edge in the same direction, but are both not
638               // touching anything else. Both get removed by the open
639               // group removal code, leaving an edge map with zero
640               // edges. The edge in the opposite direction does not
641               // exist, because there's no face that adjoins either of
642               // the two open faces.
643               continue;
644             }
645
646             for (edgelist_t::iterator j = (*e2i).second.begin(); j != (*e2i).second.end(); ++j) {
647               open_groups.insert(faceGroupID(*j));
648             }
649             complex_edges.erase(e1i);
650             complex_edges.erase(e2i);
651           }
652         }
653       }
654
655
656
657       void FaceStitcher::extractConnectedEdges(std::vector<const vertex_t *>::iterator begin,
658                                                std::vector<const vertex_t *>::iterator end,
659                                                std::vector<std::vector<Edge<3> *> > &efwd,
660                                                std::vector<std::vector<Edge<3> *> > &erev) {
661         const size_t N = std::distance(begin, end) - 1;
662
663         std::vector<const vertex_t *>::iterator e1, e2;
664         e1 = e2 = begin; ++e2;
665         vpair_t start_f = vpair_t(*e1, *e2);
666         vpair_t start_r = vpair_t(*e2, *e1);
667
668         const size_t Nfwd = complex_edges[start_f].size();
669         const size_t Nrev = complex_edges[start_r].size();
670
671         size_t j;
672         edgelist_t::iterator ji;
673
674         efwd.clear(); efwd.resize(Nfwd);
675         erev.clear(); erev.resize(Nrev);
676
677         for (j = 0, ji = complex_edges[start_f].begin();
678              ji != complex_edges[start_f].end();
679              ++j, ++ji) {
680           efwd[j].reserve(N);
681           efwd[j].push_back(*ji);
682         }
683
684         for (j = 0, ji = complex_edges[start_r].begin();
685              ji != complex_edges[start_r].end();
686              ++j, ++ji) {
687           erev[j].reserve(N);
688           erev[j].push_back(*ji);
689         }
690
691         std::vector<Edge<3> *> temp_f, temp_r;
692         temp_f.resize(Nfwd);
693         temp_r.resize(Nrev);
694
695         for (j = 1; j < N; ++j) {
696           ++e1; ++e2;
697           vpair_t ef = vpair_t(*e1, *e2);
698           vpair_t er = vpair_t(*e2, *e1);
699
700           if (complex_edges[ef].size() != Nfwd || complex_edges[ef].size() != Nrev) break;
701
702           for (size_t k = 0; k < Nfwd; ++k) {
703             Edge<3> *e_next = efwd[k].back()->perimNext();
704             CARVE_ASSERT(e_next == NULL || e_next->rev == NULL);
705             if (e_next == NULL || e_next->v2() != *e2) goto done;
706             CARVE_ASSERT(e_next->v1() == *e1);
707             CARVE_ASSERT(std::find(complex_edges[ef].begin(), complex_edges[ef].end(), e_next) != complex_edges[ef].end());
708             temp_f[k] = e_next;
709           }
710
711           for (size_t k = 0; k < Nrev; ++k) {
712             Edge<3> *e_next = erev[k].back()->perimPrev();
713             if (e_next == NULL || e_next->v1() != *e2) goto done;
714             CARVE_ASSERT(e_next->v2() == *e1);
715             CARVE_ASSERT(std::find(complex_edges[er].begin(), complex_edges[er].end(), e_next) != complex_edges[er].end());
716             temp_r[k] = e_next;
717           }
718
719           for (size_t k = 0; k < Nfwd; ++k) {
720             efwd[k].push_back(temp_f[k]);
721           }
722
723           for (size_t k = 0; k < Nrev; ++k) {
724             erev[k].push_back(temp_r[k]);
725           }
726         }
727       done:;
728       }
729
730
731
732       void FaceStitcher::construct() {
733         matchSimpleEdges();
734         if (!complex_edges.size()) return;
735
736         resolveOpenEdges();
737         if (!complex_edges.size()) return;
738
739         buildEdgeGraph(complex_edges);
740
741         std::list<std::vector<const vertex_t *> > paths;
742
743         while (edge_graph.size()) {
744           paths.push_back(std::vector<const vertex_t *>());
745           extractPath(paths.back());
746           removePath(paths.back());
747         };
748
749
750         for (std::list<std::vector<const vertex_t *> >::iterator path = paths.begin(); path != paths.end(); ++path) {
751           for (size_t i = 0; i < (*path).size() - 1;) {
752             std::vector<std::vector<Edge<3> *> > efwd, erev;
753
754             extractConnectedEdges((*path).begin() + i, (*path).end(), efwd, erev);
755
756             std::vector<std::vector<EdgeOrderData> > orderings;
757             orderForwardAndReverseEdges(efwd, erev, orderings);
758
759             matchOrderedEdges(orderings.begin(), orderings.end(), efwd, erev);
760             i += efwd[0].size();
761           }
762         }
763       }
764
765       FaceStitcher::FaceStitcher(const MeshOptions &_opts) : opts(_opts) {
766       }
767     }
768   }
769
770
771
772
773   // construct a MeshSet from a Polyhedron, maintaining on the
774   // connectivity information in the Polyhedron.
775   mesh::MeshSet<3> *meshFromPolyhedron(const poly::Polyhedron *poly, int manifold_id) {
776     typedef mesh::Vertex<3> vertex_t;
777     typedef mesh::Edge<3> edge_t;
778     typedef mesh::Face<3> face_t;
779     typedef mesh::Mesh<3> mesh_t;
780     typedef mesh::MeshSet<3> meshset_t;
781
782     std::vector<vertex_t> vertex_storage;
783     vertex_storage.reserve(poly->vertices.size());
784     for (size_t i = 0; i < poly->vertices.size(); ++i) {
785       vertex_storage.push_back(vertex_t(poly->vertices[i].v));
786     }
787
788     std::vector<std::vector<face_t *> > faces;
789     faces.resize(poly->manifold_is_closed.size());
790
791     std::unordered_map<std::pair<size_t, size_t>, std::list<edge_t *> > vertex_to_edge;
792
793     std::vector<vertex_t *> vert_ptrs;
794     for (size_t i = 0; i < poly->faces.size(); ++i) {
795       const poly::Polyhedron::face_t &src = poly->faces[i];
796       if (manifold_id != -1 && src.manifold_id != manifold_id) continue;
797       vert_ptrs.clear();
798       vert_ptrs.reserve(src.nVertices());
799       for (size_t j = 0; j < src.nVertices(); ++j) {
800         size_t vi = poly->vertexToIndex_fast(src.vertex(j));
801         vert_ptrs.push_back(&vertex_storage[vi]);
802       }
803       face_t *face = new face_t(vert_ptrs.begin(), vert_ptrs.end());
804       face->id = src.manifold_id;
805       faces[src.manifold_id].push_back(face);
806
807       edge_t *edge = face->edge;
808       do {
809         vertex_to_edge[std::make_pair(size_t(edge->v1() - &vertex_storage[0]),
810                                       size_t(edge->v2() - &vertex_storage[0]))].push_back(edge);
811         edge = edge->next;
812       } while (edge != face->edge);
813     }
814
815     // copy connectivity from Polyhedron.
816     for (size_t i = 0; i < poly->edges.size(); ++i) {
817       const poly::Polyhedron::edge_t &src = poly->edges[i];
818       size_t v1i = poly->vertexToIndex_fast(src.v1);
819       size_t v2i = poly->vertexToIndex_fast(src.v2);
820
821       std::list<edge_t *> &efwd = vertex_to_edge[std::make_pair(v1i, v2i)];
822       std::list<edge_t *> &erev = vertex_to_edge[std::make_pair(v2i, v1i)];
823
824       const std::vector<const poly::Polyhedron::face_t *> &facepairs = poly->connectivity.edge_to_face[i];
825       for (size_t j = 0; j < facepairs.size(); j += 2) {
826         const poly::Polyhedron::face_t *fa, *fb;
827         fa = facepairs[j];
828         fb = facepairs[j+1];
829         if (!fa || !fb) continue;
830         CARVE_ASSERT(fa->manifold_id == fb->manifold_id);
831         if (manifold_id != -1 && fa->manifold_id != manifold_id) continue;
832
833         std::list<edge_t *>::iterator efwdi, erevi;
834         for (efwdi = efwd.begin(); efwdi != efwd.end() && (*efwdi)->face->id != (size_t)fa->manifold_id; ++efwdi);
835         for (erevi = erev.begin(); erevi != erev.end() && (*erevi)->face->id != (size_t)fa->manifold_id; ++erevi);
836         CARVE_ASSERT(efwdi != efwd.end() && erevi != erev.end());
837
838         (*efwdi)->rev = (*erevi);
839         (*erevi)->rev = (*efwdi);
840       }
841     }
842
843     std::vector<mesh_t *> meshes;
844     meshes.reserve(faces.size());
845     for (size_t i = 0; i < faces.size(); ++i) {
846       if (faces[i].size()) {
847         meshes.push_back(new mesh_t(faces[i]));
848       }
849     }
850
851     return new meshset_t(vertex_storage, meshes);
852   }
853
854
855
856   static void copyMeshFaces(const mesh::Mesh<3> *mesh,
857                             size_t manifold_id,
858                             const mesh::Vertex<3> *Vbase,
859                             poly::Polyhedron *poly,
860                             std::unordered_map<std::pair<size_t, size_t>, std::list<mesh::Edge<3> *> > &edges,
861                             std::unordered_map<const mesh::Face<3> *, size_t> &face_map) {
862     std::vector<const poly::Polyhedron::vertex_t *> vert_ptr;
863     for (size_t f = 0; f < mesh->faces.size(); ++f) {
864       mesh::Face<3> *src = mesh->faces[f];
865       vert_ptr.clear();
866       vert_ptr.reserve(src->nVertices());
867       mesh::Edge<3> *e = src->edge;
868       do {
869         vert_ptr.push_back(&poly->vertices[e->vert - Vbase]);
870         edges[std::make_pair(e->v1() - Vbase, e->v2() - Vbase)].push_back(e);
871         e = e->next;
872       } while (e != src->edge);
873       
874       face_map[src] = poly->faces.size();;
875       
876       poly->faces.push_back(poly::Polyhedron::face_t(vert_ptr));
877       poly->faces.back().manifold_id = manifold_id;
878       poly->faces.back().owner = poly;
879     }
880   }
881
882
883
884   // construct a Polyhedron from a MeshSet
885   poly::Polyhedron *polyhedronFromMesh(const mesh::MeshSet<3> *mesh, int manifold_id) {
886     typedef poly::Polyhedron::vertex_t vertex_t;
887     typedef poly::Polyhedron::edge_t edge_t;
888     typedef poly::Polyhedron::face_t face_t;
889
890     poly::Polyhedron *poly = new poly::Polyhedron();
891     const mesh::Vertex<3> *Vbase = &mesh->vertex_storage[0];
892
893     poly->vertices.reserve(mesh->vertex_storage.size());
894     for (size_t i = 0; i < mesh->vertex_storage.size(); ++i) {
895       poly->vertices.push_back(vertex_t(mesh->vertex_storage[i].v));
896       poly->vertices.back().owner = poly;
897     }
898
899     size_t n_faces = 0;
900     if (manifold_id == -1) {
901       poly->manifold_is_closed.resize(mesh->meshes.size());
902       poly->manifold_is_negative.resize(mesh->meshes.size());
903       for (size_t m = 0; m < mesh->meshes.size(); ++m) {
904         n_faces += mesh->meshes[m]->faces.size();
905         poly->manifold_is_closed[m] = mesh->meshes[m]->isClosed();
906         poly->manifold_is_negative[m] = mesh->meshes[m]->isNegative();
907       }
908     } else {
909       poly->manifold_is_closed.resize(1);
910       poly->manifold_is_negative.resize(1);
911       n_faces = mesh->meshes[manifold_id]->faces.size();
912       poly->manifold_is_closed[manifold_id] = mesh->meshes[manifold_id]->isClosed();
913       poly->manifold_is_negative[manifold_id] = mesh->meshes[manifold_id]->isNegative();
914     }
915
916     std::unordered_map<std::pair<size_t, size_t>, std::list<mesh::Edge<3> *> > edges;
917     std::unordered_map<const mesh::Face<3> *, size_t> face_map;
918     poly->faces.reserve(n_faces);
919
920     if (manifold_id == -1) {
921       for (size_t m = 0; m < mesh->meshes.size(); ++m) {
922         copyMeshFaces(mesh->meshes[m], m, Vbase, poly, edges, face_map);
923       }
924     } else {
925       copyMeshFaces(mesh->meshes[manifold_id], 0, Vbase, poly, edges, face_map);
926     }
927
928     size_t n_edges = 0;
929     for (std::unordered_map<std::pair<size_t, size_t>, std::list<mesh::Edge<3> *> >::iterator i = edges.begin(); i != edges.end(); ++i) {
930       if ((*i).first.first < (*i).first.second || edges.find(std::make_pair((*i).first.second, (*i).first.first)) == edges.end()) {
931         n_edges++;
932       }
933     }
934
935     poly->edges.reserve(n_edges);
936     for (std::unordered_map<std::pair<size_t, size_t>, std::list<mesh::Edge<3> *> >::iterator i = edges.begin(); i != edges.end(); ++i) {
937       if ((*i).first.first < (*i).first.second ||
938           edges.find(std::make_pair((*i).first.second, (*i).first.first)) == edges.end()) {
939         poly->edges.push_back(edge_t(&poly->vertices[(*i).first.first],
940                                      &poly->vertices[(*i).first.second],
941                                      poly));
942       }
943     }
944
945     poly->initVertexConnectivity();
946
947     // build edge entries for face.
948     for (size_t f = 0; f < poly->faces.size(); ++f) {
949       face_t &face = poly->faces[f];
950       size_t N = face.nVertices();
951       for (size_t v = 0; v < N; ++v) {
952         size_t v1i = poly->vertexToIndex_fast(face.vertex(v));
953         size_t v2i = poly->vertexToIndex_fast(face.vertex((v+1)%N));
954         std::vector<const edge_t *> found_edge;
955         std::set_intersection(poly->connectivity.vertex_to_edge[v1i].begin(), poly->connectivity.vertex_to_edge[v1i].end(),
956                               poly->connectivity.vertex_to_edge[v2i].begin(), poly->connectivity.vertex_to_edge[v2i].end(),
957                               std::back_inserter(found_edge));
958         CARVE_ASSERT(found_edge.size() == 1);
959         face.edge(v) = found_edge[0];
960       }
961     }
962
963     poly->connectivity.edge_to_face.resize(poly->edges.size());
964
965     for (size_t i = 0; i < poly->edges.size(); ++i) {
966       size_t v1i = poly->vertexToIndex_fast(poly->edges[i].v1);
967       size_t v2i = poly->vertexToIndex_fast(poly->edges[i].v2);
968       std::list<mesh::Edge<3> *> &efwd = edges[std::make_pair(v1i, v2i)];
969       std::list<mesh::Edge<3> *> &erev = edges[std::make_pair(v1i, v2i)];
970
971       for (std::list<mesh::Edge<3> *>::iterator j = efwd.begin(); j != efwd.end(); ++j) {
972         mesh::Edge<3> *edge = *j;
973         if (face_map.find(edge->face) != face_map.end()) {
974           poly->connectivity.edge_to_face[i].push_back(&poly->faces[face_map[edge->face]]);
975           if (edge->rev == NULL) {
976             poly->connectivity.edge_to_face[i].push_back(NULL);
977           } else {
978             poly->connectivity.edge_to_face[i].push_back(&poly->faces[face_map[edge->rev->face]]);
979           }
980         }
981       }
982       for (std::list<mesh::Edge<3> *>::iterator j = erev.begin(); j != erev.end(); ++j) {
983         mesh::Edge<3> *edge = *j;
984         if (face_map.find(edge->face) != face_map.end()) {
985           if (edge->rev == NULL) {
986             poly->connectivity.edge_to_face[i].push_back(NULL);
987             poly->connectivity.edge_to_face[i].push_back(&poly->faces[face_map[edge->face]]);
988           }
989         }
990       }
991
992     }
993
994     poly->initSpatialIndex();
995
996     // XXX: at this point, manifold_is_negative is not set up. This
997     // info should be computed/stored in Mesh instances.
998
999     return poly;
1000   }
1001
1002
1003
1004 }
1005
1006
1007
1008 // explicit instantiation for 2D case.
1009 // XXX: do not compile because of a missing definition for fitPlane in the 2d case.
1010
1011 // template class carve::mesh::Vertex<2>;
1012 // template class carve::mesh::Edge<2>;
1013 // template class carve::mesh::Face<2>;
1014 // template class carve::mesh::Mesh<2>;
1015 // template class carve::mesh::MeshSet<2>;
1016
1017 // explicit instantiation for 3D case.
1018 template class carve::mesh::Vertex<3>;
1019 template class carve::mesh::Edge<3>;
1020 template class carve::mesh::Face<3>;
1021 template class carve::mesh::Mesh<3>;
1022 template class carve::mesh::MeshSet<3>;
1023
1024
1025
1026 carve::PointClass carve::mesh::classifyPoint(
1027     const carve::mesh::MeshSet<3> *meshset,
1028     const carve::geom::RTreeNode<3, carve::mesh::Face<3> *> *face_rtree,
1029     const carve::geom::vector<3> &v,
1030     bool even_odd,
1031     const carve::mesh::Mesh<3> *mesh,
1032     const carve::mesh::Face<3> **hit_face) {
1033
1034   if (hit_face) *hit_face = NULL;
1035
1036 #if defined(DEBUG_CONTAINS_VERTEX)
1037   std::cerr << "{containsVertex " << v << "}" << std::endl;
1038 #endif
1039
1040   if (!face_rtree->bbox.containsPoint(v)) {
1041 #if defined(DEBUG_CONTAINS_VERTEX)
1042     std::cerr << "{final:OUT(aabb short circuit)}" << std::endl;
1043 #endif
1044     // XXX: if the top level manifolds are negative, this should be POINT_IN.
1045     // for the moment, this only works for a single manifold.
1046     if (meshset->meshes.size() == 1 && meshset->meshes[0]->isNegative()) {
1047       return POINT_IN;
1048     }
1049     return POINT_OUT;
1050   }
1051
1052   std::vector<carve::mesh::Face<3> *> near_faces;
1053   face_rtree->search(v, std::back_inserter(near_faces));
1054
1055   for (size_t i = 0; i < near_faces.size(); i++) {
1056     if (mesh != NULL && mesh != near_faces[i]->mesh) continue;
1057
1058     // XXX: Do allow the tested vertex to be ON an open
1059     // manifold. This was here originally because of the
1060     // possibility of an open manifold contained within a closed
1061     // manifold.
1062
1063     // if (!near_faces[i]->mesh->isClosed()) continue;
1064
1065     if (near_faces[i]->containsPoint(v)) {
1066 #if defined(DEBUG_CONTAINS_VERTEX)
1067       std::cerr << "{final:ON(hits face " << near_faces[i] << ")}" << std::endl;
1068 #endif
1069       if (hit_face) *hit_face = near_faces[i];
1070       return POINT_ON;
1071     }
1072   }
1073
1074   double ray_len = face_rtree->bbox.extent.length() * 2;
1075
1076
1077   std::vector<std::pair<const carve::mesh::Face<3> *, carve::geom::vector<3> > > manifold_intersections;
1078
1079   for (;;) {
1080     double a1 = random() / double(RAND_MAX) * M_TWOPI;
1081     double a2 = random() / double(RAND_MAX) * M_TWOPI;
1082
1083     carve::geom3d::Vector ray_dir = carve::geom::VECTOR(sin(a1) * sin(a2), cos(a1) * sin(a2), cos(a2));
1084
1085 #if defined(DEBUG_CONTAINS_VERTEX)
1086     std::cerr << "{testing ray: " << ray_dir << "}" << std::endl;
1087 #endif
1088
1089     carve::geom::vector<3> v2 = v + ray_dir * ray_len;
1090
1091     bool failed = false;
1092     carve::geom::linesegment<3> line(v, v2);
1093     carve::geom::vector<3> intersection;
1094
1095     near_faces.clear();
1096     manifold_intersections.clear();
1097     face_rtree->search(line, std::back_inserter(near_faces));
1098
1099     for (unsigned i = 0; !failed && i < near_faces.size(); i++) {
1100       if (mesh != NULL && mesh != near_faces[i]->mesh) continue;
1101
1102       if (!near_faces[i]->mesh->isClosed()) continue;
1103
1104       switch (near_faces[i]->lineSegmentIntersection(line, intersection)) {
1105       case INTERSECT_FACE: {
1106
1107 #if defined(DEBUG_CONTAINS_VERTEX)
1108         std::cerr << "{intersects face: " << near_faces[i]
1109                   << " dp: " << dot(ray_dir, near_faces[i]->plane.N) << "}" << std::endl;
1110 #endif
1111
1112         if (!even_odd && fabs(dot(ray_dir, near_faces[i]->plane.N)) < EPSILON) {
1113
1114 #if defined(DEBUG_CONTAINS_VERTEX)
1115           std::cerr << "{failing(small dot product)}" << std::endl;
1116 #endif
1117
1118           failed = true;
1119           break;
1120         }
1121         manifold_intersections.push_back(std::make_pair(near_faces[i], intersection));
1122         break;
1123       }
1124       case INTERSECT_NONE: {
1125         break;
1126       }
1127       default: {
1128
1129 #if defined(DEBUG_CONTAINS_VERTEX)
1130         std::cerr << "{failing(degenerate intersection)}" << std::endl;
1131 #endif
1132         failed = true;
1133         break;
1134       }
1135       }
1136     }
1137
1138     if (!failed) {
1139       if (even_odd) {
1140         return (manifold_intersections.size() & 1) ? POINT_IN : POINT_OUT;
1141       }
1142
1143 #if defined(DEBUG_CONTAINS_VERTEX)
1144       std::cerr << "{intersections ok [count:"
1145                 << manifold_intersections.size()
1146                 << "], sorting}"
1147                 << std::endl;
1148 #endif
1149
1150       carve::geom3d::sortInDirectionOfRay(ray_dir,
1151                                           manifold_intersections.begin(),
1152                                           manifold_intersections.end(),
1153                                           carve::geom3d::vec_adapt_pair_second());
1154
1155       std::map<const carve::mesh::Mesh<3> *, int> crossings;
1156
1157       for (size_t i = 0; i < manifold_intersections.size(); ++i) {
1158         const carve::mesh::Face<3> *f = manifold_intersections[i].first;
1159         if (dot(ray_dir, f->plane.N) < 0.0) {
1160           crossings[f->mesh]++;
1161         } else {
1162           crossings[f->mesh]--;
1163         }
1164       }
1165
1166 #if defined(DEBUG_CONTAINS_VERTEX)
1167       for (std::map<const carve::mesh::Mesh<3> *, int>::const_iterator i = crossings.begin(); i != crossings.end(); ++i) {
1168         std::cerr << "{mesh " << (*i).first << " crossing count: " << (*i).second << "}" << std::endl;
1169       }
1170 #endif
1171
1172       for (size_t i = 0; i < manifold_intersections.size(); ++i) {
1173         const carve::mesh::Face<3> *f = manifold_intersections[i].first;
1174
1175 #if defined(DEBUG_CONTAINS_VERTEX)
1176         std::cerr << "{intersection at "
1177                   << manifold_intersections[i].second
1178                   << " mesh: "
1179                   << f->mesh
1180                   << " count: "
1181                   << crossings[f->mesh]
1182                   << "}"
1183                   << std::endl;
1184 #endif
1185
1186         if (crossings[f->mesh] < 0) {
1187           // inside this manifold.
1188
1189 #if defined(DEBUG_CONTAINS_VERTEX)
1190           std::cerr << "{final:IN}" << std::endl;
1191 #endif
1192
1193           return POINT_IN;
1194         } else if (crossings[f->mesh] > 0) {
1195           // outside this manifold, but it's an infinite manifold. (for instance, an inverted cube)
1196
1197 #if defined(DEBUG_CONTAINS_VERTEX)
1198           std::cerr << "{final:OUT}" << std::endl;
1199 #endif
1200
1201           return POINT_OUT;
1202         }
1203       }
1204
1205 #if defined(DEBUG_CONTAINS_VERTEX)
1206       std::cerr << "{final:OUT(default)}" << std::endl;
1207 #endif
1208
1209       return POINT_OUT;
1210     }
1211   }
1212 }
1213
1214
1215