Bundle latest version of Carve library which shall resolve compilation issues with...
[blender.git] / extern / carve / include / carve / mesh_impl.hpp
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 #pragma once
19
20 #include <carve/geom2d.hpp>
21 #include <carve/geom3d.hpp>
22 #include <carve/djset.hpp>
23
24 #include <iostream>
25 #include <deque>
26
27 #include <stddef.h>
28
29 namespace carve {
30   namespace mesh {
31
32
33
34     namespace detail {
35       template<typename list_t>
36       struct list_iter_t {
37         typedef std::bidirectional_iterator_tag iterator_category;
38         typedef list_t                          value_type;
39         typedef ptrdiff_t                       difference_type;
40         typedef value_type &                    reference;
41         typedef value_type *                    pointer;
42
43         list_t *curr;
44         int pos;
45
46         list_iter_t() { }
47         list_iter_t(list_t *_curr, int _pos) : curr(_curr), pos(_pos) { }
48
49         list_iter_t operator++(int) { list_iter_t result(*this); ++pos; curr = curr->next; return result; }
50         list_iter_t operator--(int) { list_iter_t result(*this); --pos; curr = curr->prev; return result; }
51
52         list_iter_t operator++() { ++pos; curr = curr->next; return *this; }
53         list_iter_t operator--() { --pos; curr = curr->prev; return *this; }
54
55         bool operator==(const list_iter_t &other) const { return curr == other.curr && pos == other.pos; }
56         bool operator!=(const list_iter_t &other) const { return curr != other.curr || pos != other.pos; }
57
58         reference operator*() { return *curr; }
59         pointer operator->() { return curr; }
60
61         int idx() const { return pos; }
62       };
63     }
64
65
66
67     template<unsigned ndim>
68     Edge<ndim> *Edge<ndim>::mergeFaces() {
69       if (rev == NULL) return NULL;
70
71       face_t *fwdface = face;
72       face_t *revface = rev->face;
73
74       size_t n_removed = 0;
75
76       Edge *splice_beg = this;
77       do {
78         splice_beg = splice_beg->prev;
79         ++n_removed;
80       } while (splice_beg != this &&
81                splice_beg->rev &&
82                splice_beg->next->rev->prev == splice_beg->rev);
83
84       if (splice_beg == this) {
85         // edge loops are completely matched.
86         return NULL;
87       }
88
89       Edge *splice_end = this;
90       do {
91         splice_end = splice_end->next;
92         ++n_removed;
93       } while (splice_end->rev &&
94                splice_end->prev->rev->next == splice_end->rev);
95
96       --n_removed;
97
98       Edge *link1_p = splice_beg;
99       Edge *link1_n = splice_beg->next->rev->next;
100
101       Edge *link2_p = splice_end->prev->rev->prev;
102       Edge *link2_n = splice_end;
103
104       CARVE_ASSERT(link1_p->face == fwdface);
105       CARVE_ASSERT(link1_n->face == revface);
106
107       CARVE_ASSERT(link2_p->face == revface);
108       CARVE_ASSERT(link2_n->face == fwdface);
109
110       Edge *left_loop = link1_p->next;
111
112       CARVE_ASSERT(left_loop->rev == link1_n->prev);
113
114       _link(link2_n->prev, link1_p->next);
115       _link(link1_n->prev, link2_p->next);
116
117       _link(link1_p, link1_n);
118       _link(link2_p, link2_n);
119
120       fwdface->edge = link1_p;
121
122       for (Edge *e = link1_n; e != link2_n; e = e->next) {
123         CARVE_ASSERT(e->face == revface);
124         e->face = fwdface;
125         fwdface->n_edges++;
126       }
127       for (Edge *e = link2_n; e != link1_n; e = e->next) {
128         CARVE_ASSERT(e->face == fwdface);
129       }
130
131       fwdface->n_edges -= n_removed;
132
133       revface->n_edges = 0;
134       revface->edge = NULL;
135
136       _setloopface(left_loop, NULL);
137       _setloopface(left_loop->rev, NULL);
138
139       return left_loop;
140     }
141
142
143
144     template<unsigned ndim>
145     Edge<ndim> *Edge<ndim>::removeHalfEdge() {
146       Edge *n = NULL;
147       if (face) {
148         face->n_edges--;
149       }
150
151       if (next == this) {
152         if (face) face->edge = NULL;
153       } else {
154         if (face && face->edge == this) face->edge = next;
155         next->prev = prev;
156         prev->next = next;
157         n = next;
158       }
159       delete this;
160       return n;
161     }
162
163
164
165     template<unsigned ndim>
166     Edge<ndim> *Edge<ndim>::removeEdge() {
167       if (rev) {
168         rev->removeHalfEdge();
169       }
170       return removeHalfEdge();
171     }
172
173
174
175     template<unsigned ndim>
176     void Edge<ndim>::unlink() {
177       if (rev) { rev->rev = NULL; rev = NULL; }
178       if (prev->rev) { prev->rev->rev = NULL; prev->rev = NULL; }
179
180       if (face) {
181         face->n_edges--;
182         if (face->edge == this) face->edge = next;
183         face = NULL;
184       }
185
186       next->prev = prev;
187       prev->next = next;
188
189       prev = next = this;
190     }
191
192
193
194     template<unsigned ndim>
195     void Edge<ndim>::insertBefore(Edge<ndim> *other) {
196       if (prev != this) unlink();
197       prev = other->prev;
198       next = other;
199       next->prev = this;
200       prev->next = this;
201
202       if (prev->rev) { prev->rev->rev = NULL;  prev->rev = NULL; }
203     }
204
205
206
207     template<unsigned ndim>
208     void Edge<ndim>::insertAfter(Edge<ndim> *other) {
209       if (prev != this) unlink();
210       next = other->next;
211       prev = other;
212       next->prev = this;
213       prev->next = this;
214
215       if (prev->rev) { prev->rev->rev = NULL;  prev->rev = NULL; }
216     }
217
218
219
220     template<unsigned ndim>
221     size_t Edge<ndim>::loopSize() const {
222       const Edge *e = this;
223       size_t n = 0;
224       do { e = e->next; ++n; } while (e != this);
225       return n;
226     }
227
228
229
230     template<unsigned ndim>
231     Edge<ndim> *Edge<ndim>::perimNext() const {
232       if (rev) return NULL;
233       Edge *e = next;
234       while(e->rev) {
235         e = e->rev->next;
236       }
237       return e;
238     }
239
240
241
242     template<unsigned ndim>
243     Edge<ndim> *Edge<ndim>::perimPrev() const {
244       if (rev) return NULL;
245       Edge *e = prev;
246       while(e->rev) {
247         e = e->rev->prev;
248       }
249       return e;
250     }
251
252
253
254     template<unsigned ndim>
255     Edge<ndim>::Edge(vertex_t *_vert, face_t *_face) :
256         vert(_vert), face(_face), prev(NULL), next(NULL), rev(NULL) {
257       prev = next = this;
258     }
259
260
261
262     template<unsigned ndim>
263     Edge<ndim>::~Edge() {
264     }
265
266
267
268     template<unsigned ndim>
269     typename Face<ndim>::aabb_t Face<ndim>::getAABB() const {
270       aabb_t aabb;
271       aabb.fit(begin(), end(), vector_mapping());
272       return aabb;
273     }
274
275
276
277     template<unsigned ndim>
278     bool Face<ndim>::recalc() {
279       if (!carve::geom3d::fitPlane(begin(), end(), vector_mapping(), plane)) {
280         return false;
281       }
282
283       int da = carve::geom::largestAxis(plane.N);
284       double A = carve::geom2d::signedArea(begin(), end(), projection_mapping(getProjector(false, da)));
285
286       if ((A < 0.0) ^ (plane.N.v[da] < 0.0)) {
287         plane.negate();
288       }
289
290       project = getProjector(plane.N.v[da] > 0, da);
291       unproject = getUnprojector(plane.N.v[da] > 0, da);
292
293       return true;
294     }
295
296
297
298     template<unsigned ndim>
299     void Face<ndim>::clearEdges() {
300       if (!edge) return;
301
302       edge_t *curr = edge;
303       do {
304         edge_t *next = curr->next;
305         delete curr;
306         curr = next;
307       } while (curr != edge);
308
309       edge = NULL;
310
311       n_edges = 0;
312     }
313
314
315
316     template<unsigned ndim>
317     template<typename iter_t>
318     void Face<ndim>::loopFwd(iter_t begin, iter_t end) {
319       clearEdges();
320       if (begin == end) return;
321       edge = new edge_t(*begin, this); ++n_edges; ++begin;
322       while (begin != end) {
323         edge_t *e = new edge_t(*begin, this);
324         e->insertAfter(edge->prev);
325         ++n_edges;
326         ++begin;
327       }
328     }
329
330
331
332     template<unsigned ndim>
333     template<typename iter_t>
334     void Face<ndim>::loopRev(iter_t begin, iter_t end) {
335       clearEdges();
336       if (begin == end) return;
337       edge = new edge_t(*begin, this); ++n_edges; ++begin;
338       while (begin != end) {
339         edge_t *e = new edge_t(*begin, this);
340         e->insertBefore(edge->next);
341         ++n_edges;
342         ++begin;
343       }
344     }
345
346
347
348     template<unsigned ndim>
349     template<typename iter_t>
350     void Face<ndim>::init(iter_t begin, iter_t end) {
351       loopFwd(begin, end);
352     }
353
354
355
356     template<unsigned ndim>
357     void Face<ndim>::init(vertex_t *a, vertex_t *b, vertex_t *c) {
358       clearEdges();
359       edge_t *ea = new edge_t(a, this);
360       edge_t *eb = new edge_t(b, this);
361       edge_t *ec = new edge_t(c, this);
362       eb->insertAfter(ea);
363       ec->insertAfter(eb);
364       edge = ea;
365       n_edges = 3;
366     }
367
368
369
370     template<unsigned ndim>
371     void Face<ndim>::init(vertex_t *a, vertex_t *b, vertex_t *c, vertex_t *d) {
372       clearEdges();
373       edge_t *ea = new edge_t(a, this);
374       edge_t *eb = new edge_t(b, this);
375       edge_t *ec = new edge_t(c, this);
376       edge_t *ed = new edge_t(d, this);
377       eb->insertAfter(ea);
378       ec->insertAfter(eb);
379       ed->insertAfter(ec);
380       edge = ea;
381       n_edges = 4;
382     }
383
384
385
386     template<unsigned ndim>
387     void Face<ndim>::getVertices(std::vector<vertex_t *> &verts) const {
388       verts.clear();
389       verts.reserve(n_edges);
390       const edge_t *e = edge;
391       do { verts.push_back(e->vert); e = e->next; } while (e != edge);
392     }
393
394
395
396     template<unsigned ndim>
397     void Face<ndim>::getProjectedVertices(std::vector<carve::geom::vector<2> > &verts) const {
398       verts.clear();
399       verts.reserve(n_edges);
400       const edge_t *e = edge;
401       do { verts.push_back(project(e->vert->v)); e = e->next; } while (e != edge);
402     }
403
404
405
406     template<unsigned ndim>
407     typename Face<ndim>::vector_t Face<ndim>::centroid() const {
408       vector_t v;
409       edge_t *e = edge;
410       do {
411         v += e->vert->v;
412         e = e->next;
413       } while(e != edge);
414       v /= n_edges;
415       return v;
416     }
417
418
419
420     template<unsigned ndim>
421     void Face<ndim>::canonicalize() {
422       edge_t *min = edge;
423       edge_t *e = edge;
424
425       do {
426         if (e->vert < min->vert) min = e;
427         e = e->next;
428       } while (e != edge);
429
430       edge = min;
431     }
432
433
434
435     template<unsigned ndim>
436     template<typename iter_t>
437     Face<ndim> *Face<ndim>::create(iter_t beg, iter_t end, bool reversed) const {
438       Face *r = new Face();
439
440       if (reversed) {
441         r->loopRev(beg, end);
442         r->plane = -plane;
443       } else {
444         r->loopFwd(beg, end);
445         r->plane = plane;
446       }
447
448       int da = carve::geom::largestAxis(r->plane.N);
449
450       r->project = r->getProjector(r->plane.N.v[da] > 0, da);
451       r->unproject = r->getUnprojector(r->plane.N.v[da] > 0, da);
452
453       return r;
454     }
455
456
457
458     template<unsigned ndim>
459     Face<ndim> *Face<ndim>::clone(const vertex_t *old_base,
460                                   vertex_t *new_base,
461                                   std::unordered_map<const edge_t *, edge_t *> &edge_map) const {
462       Face *r = new Face(*this);
463
464       edge_t *e = edge;
465       edge_t *r_p = NULL;
466       edge_t *r_e;
467       do {
468         r_e = new edge_t(e->vert - old_base + new_base, r);
469         edge_map[e] = r_e;
470         if (r_p) {
471           r_p->next = r_e;
472           r_e->prev = r_p;
473         } else {
474           r->edge = r_e;
475         }
476         r_p = r_e;
477
478         if (e->rev) {
479           typename std::unordered_map<const edge_t *, edge_t *>::iterator rev_i = edge_map.find(e->rev);
480           if (rev_i != edge_map.end()) {
481             r_e->rev = (*rev_i).second;
482             (*rev_i).second->rev = r_e;
483           }
484         }
485
486         e = e->next;
487       } while (e != edge);
488       r_e->next = r->edge;
489       r->edge->prev = r_e;
490       return r;
491     }
492
493
494
495     template<unsigned ndim>
496     Mesh<ndim>::Mesh(std::vector<face_t *> &_faces,
497                      std::vector<edge_t *> &_open_edges,
498                      std::vector<edge_t *> &_closed_edges,
499                      bool _is_negative) {
500       std::swap(faces, _faces);
501       std::swap(open_edges, _open_edges);
502       std::swap(closed_edges, _closed_edges);
503       is_negative = _is_negative;
504       meshset = NULL;
505
506       for (size_t i = 0; i < faces.size(); ++i) {
507         faces[i]->mesh = this;
508       }
509     }
510
511
512
513     namespace detail {
514       template<typename iter_t>
515       void FaceStitcher::initEdges(iter_t begin,
516                                    iter_t end) {
517         size_t c = 0;
518         for (iter_t i = begin; i != end; ++i) {
519           face_t *face = *i;
520           CARVE_ASSERT(face->mesh == NULL); // for the moment, can only insert a face into a mesh once.
521
522           face->id = c++;
523           edge_t *e = face->edge;
524           do {
525             edges[vpair_t(e->v1(), e->v2())].push_back(e);
526             e = e->next;
527             if (e->rev) { e->rev->rev = NULL; e->rev = NULL; }
528           } while (e != face->edge);
529         }
530         face_groups.init(c);
531         is_open.clear();
532         is_open.resize(c, false);
533       }
534
535       template<typename iter_t>
536       void FaceStitcher::build(iter_t begin,
537                                iter_t end,
538                                std::vector<Mesh<3> *> &meshes) {
539         // work out what set each face belongs to, and then construct
540         // mesh instances for each set of faces.
541         std::vector<size_t> index_set;
542         std::vector<size_t> set_size;
543         face_groups.get_index_to_set(index_set, set_size);
544
545         std::vector<std::vector<face_t *> > mesh_faces;
546         mesh_faces.resize(set_size.size());
547         for (size_t i = 0; i < set_size.size(); ++i) {
548           mesh_faces[i].reserve(set_size[i]);
549         }
550       
551         for (iter_t i = begin; i != end; ++i) {
552           face_t *face = *i;
553           mesh_faces[index_set[face->id]].push_back(face);
554         }
555
556         meshes.clear();
557         meshes.reserve(mesh_faces.size());
558         for (size_t i = 0; i < mesh_faces.size(); ++i) {
559           meshes.push_back(new Mesh<3>(mesh_faces[i]));
560         }
561       }
562
563       template<typename iter_t>
564       void FaceStitcher::create(iter_t begin,
565                                 iter_t end,
566                                 std::vector<Mesh<3> *> &meshes) {
567         initEdges(begin, end);
568         construct();
569         build(begin, end, meshes);
570       }
571     }
572
573
574
575     template<unsigned ndim>
576     void Mesh<ndim>::cacheEdges() {
577       closed_edges.clear();
578       open_edges.clear();
579
580       for (size_t i = 0; i < faces.size(); ++i) {
581         face_t *face = faces[i];
582         edge_t *e = face->edge;
583         do {
584           if (e->rev == NULL) {
585             open_edges.push_back(e);
586           } else if (e < e->rev) {
587             closed_edges.push_back(e);
588           }
589           e = e->next;
590         } while (e != face->edge);
591       }
592     }
593
594
595
596     template<unsigned ndim>
597     Mesh<ndim>::Mesh(std::vector<face_t *> &_faces) : faces(), open_edges(), closed_edges(), meshset(NULL) {
598       faces.swap(_faces);
599       for (size_t i = 0; i < faces.size(); ++i) {
600         faces[i]->mesh = this;
601       }
602       cacheEdges();
603       calcOrientation();
604     }
605
606
607
608     template<unsigned ndim>
609     void Mesh<ndim>::calcOrientation() {
610       if (open_edges.size() || !closed_edges.size()) {
611         is_negative = false;
612       } else {
613         edge_t *emin = closed_edges[0];
614         if (emin->rev->v1()->v < emin->v1()->v) emin = emin->rev;
615         for (size_t i = 1; i < closed_edges.size(); ++i) {
616           if (closed_edges[i]->v1()->v      < emin->v1()->v) emin = closed_edges[i];
617           if (closed_edges[i]->rev->v1()->v < emin->v1()->v) emin = closed_edges[i]->rev;
618         }
619
620         std::vector<face_t *> min_faces;
621         edge_t *e = emin;
622         do {
623           min_faces.push_back(e->face);
624           CARVE_ASSERT(e->rev != NULL);
625           e = e->rev->next;
626           CARVE_ASSERT(e->v1() == emin->v1());
627           CARVE_ASSERT(e->v1()->v <= e->v2()->v);
628         } while (e != emin);
629
630         double max_abs_x = 0.0;
631         for (size_t f = 0; f < min_faces.size(); ++f) {
632           if (fabs(min_faces[f]->plane.N.x) > fabs(max_abs_x)) max_abs_x = min_faces[f]->plane.N.x;
633         }
634         is_negative =  max_abs_x > 0.0;
635       }
636     }
637
638
639
640     template<unsigned ndim>
641     Mesh<ndim> *Mesh<ndim>::clone(const vertex_t *old_base,
642                                   vertex_t *new_base) const {
643       std::vector<face_t *> r_faces;
644       std::vector<edge_t *> r_open_edges;
645       std::vector<edge_t *> r_closed_edges;
646       std::unordered_map<const edge_t *, edge_t *> edge_map;
647
648       r_faces.reserve(faces.size());
649       r_open_edges.reserve(r_open_edges.size());
650       r_closed_edges.reserve(r_closed_edges.size());
651
652       for (size_t i = 0; i < faces.size(); ++i) {
653         r_faces.push_back(faces[i]->clone(old_base, new_base, edge_map));
654       }
655       for (size_t i = 0; i < closed_edges.size(); ++i) {
656         r_closed_edges.push_back(edge_map[closed_edges[i]]);
657         r_closed_edges.back()->rev = edge_map[closed_edges[i]->rev];
658       }
659       for (size_t i = 0; i < open_edges.size(); ++i) {
660         r_open_edges.push_back(edge_map[open_edges[i]]);
661       }
662
663       return new Mesh(r_faces, r_open_edges, r_closed_edges, is_negative);
664     }
665
666
667
668     template<unsigned ndim>
669     Mesh<ndim>::~Mesh() {
670       for (size_t i = 0; i < faces.size(); ++i) {
671         delete faces[i];
672       }
673     }
674
675
676
677     template<unsigned ndim>
678     template<typename iter_t>
679     void Mesh<ndim>::create(iter_t begin, iter_t end, std::vector<Mesh<ndim> *> &meshes) {
680       meshes.clear();
681     }
682
683
684
685     template<>
686     template<typename iter_t>
687     void Mesh<3>::create(iter_t begin, iter_t end, std::vector<Mesh<3> *> &meshes) {
688       detail::FaceStitcher().create(begin, end, meshes);
689     }
690
691
692
693     template<unsigned ndim>
694     template<typename iter_t>
695     void MeshSet<ndim>::_init_from_faces(iter_t begin, iter_t end) {
696       typedef std::unordered_map<const vertex_t *, size_t> map_t;
697       map_t vmap;
698
699       for (iter_t i = begin; i != end; ++i) {
700         face_t *f = *i;
701         edge_t *e = f->edge;
702         do {
703           typename map_t::const_iterator j = vmap.find(e->vert);
704           if (j == vmap.end()) {
705             size_t idx = vmap.size();
706             vmap[e->vert] = idx;
707           }
708           e = e->next;
709         } while (e != f->edge);
710       }
711
712       vertex_storage.resize(vmap.size());
713       for (typename map_t::const_iterator i = vmap.begin(); i != vmap.end(); ++i) {
714         vertex_storage[(*i).second].v = (*i).first->v;
715       }
716
717       for (iter_t i = begin; i != end; ++i) {
718         face_t *f = *i;
719         edge_t *e = f->edge;
720         do {
721           e->vert = &vertex_storage[vmap[e->vert]];
722           e = e->next;
723         } while (e != f->edge);
724       }
725
726       mesh_t::create(begin, end, meshes);
727
728       for (size_t i = 0; i < meshes.size(); ++i) {
729         meshes[i]->meshset = this;
730       }
731     }
732
733
734
735     template<unsigned ndim>
736     MeshSet<ndim>::MeshSet(const std::vector<typename MeshSet<ndim>::vertex_t::vector_t> &points,
737                            size_t n_faces,
738                            const std::vector<int> &face_indices) {
739       vertex_storage.reserve(points.size());
740       std::vector<face_t *> faces;
741       faces.reserve(n_faces);
742       for (size_t i = 0; i < points.size(); ++i) {
743         vertex_storage.push_back(vertex_t(points[i]));
744       }
745
746       std::vector<vertex_t *> v;
747       size_t p = 0;
748       for (size_t i = 0; i < n_faces; ++i) {
749         const size_t N = face_indices[p++];
750         v.clear();
751         v.reserve(N);
752         for (size_t j = 0; j < N; ++j) {
753           v.push_back(&vertex_storage[face_indices[p++]]);
754         }
755         faces.push_back(new face_t(v.begin(), v.end()));
756       }
757       CARVE_ASSERT(p == face_indices.size());
758       mesh_t::create(faces.begin(), faces.end(), meshes);
759
760       for (size_t i = 0; i < meshes.size(); ++i) {
761         meshes[i]->meshset = this;
762       }
763     }
764
765
766
767     template<unsigned ndim>
768     MeshSet<ndim>::MeshSet(std::vector<face_t *> &faces) {
769       _init_from_faces(faces.begin(), faces.end());
770     }
771
772
773
774     template<unsigned ndim>
775     MeshSet<ndim>::MeshSet(std::list<face_t *> &faces) {
776       _init_from_faces(faces.begin(), faces.end());
777     }
778
779
780
781     template<unsigned ndim>
782     MeshSet<ndim>::MeshSet(std::vector<vertex_t> &_vertex_storage,
783                            std::vector<mesh_t *> &_meshes) {
784       vertex_storage.swap(_vertex_storage);
785       meshes.swap(_meshes);
786
787       for (size_t i = 0; i < meshes.size(); ++i) {
788         meshes[i]->meshset = this;
789       }
790     }
791
792
793
794     template<unsigned ndim>
795     MeshSet<ndim>::MeshSet(std::vector<typename MeshSet<ndim>::mesh_t *> &_meshes) {
796       meshes.swap(_meshes);
797       std::unordered_map<vertex_t *, size_t> vert_idx;
798
799       for (size_t m = 0; m < meshes.size(); ++m) {
800         mesh_t *mesh = meshes[m];
801         CARVE_ASSERT(mesh->meshset == NULL);
802         mesh->meshset = this;
803         for (size_t f = 0; f < mesh->faces.size(); ++f) {
804           face_t *face = mesh->faces[f];
805           edge_t *edge = face->edge;
806           do {
807             vert_idx[edge->vert] = 0;
808             edge = edge->next;
809           } while (edge != face->edge);
810         }
811       }
812
813       vertex_storage.reserve(vert_idx.size());
814       for (typename std::unordered_map<vertex_t *, size_t>::iterator i = vert_idx.begin(); i != vert_idx.end(); ++i) {
815         (*i).second = vertex_storage.size();
816         vertex_storage.push_back(*(*i).first);
817       }
818
819       for (size_t m = 0; m < meshes.size(); ++m) {
820         mesh_t *mesh = meshes[m];
821         for (size_t f = 0; f < mesh->faces.size(); ++f) {
822           face_t *face = mesh->faces[f];
823           edge_t *edge = face->edge;
824           do {
825             size_t i = vert_idx[edge->vert];
826             edge->vert = &vertex_storage[i];
827             edge = edge->next;
828           } while (edge != face->edge);
829         }
830       }
831     }
832
833
834
835     template<unsigned ndim>
836     MeshSet<ndim> *MeshSet<ndim>::clone() const {
837       std::vector<vertex_t> r_vertex_storage = vertex_storage;
838       std::vector<mesh_t *> r_meshes;
839       for (size_t i = 0; i < meshes.size(); ++i) {
840         r_meshes.push_back(meshes[i]->clone(&vertex_storage[0], &r_vertex_storage[0]));
841       }
842
843       return new MeshSet(r_vertex_storage, r_meshes);
844     }
845
846
847
848     template<unsigned ndim>
849     MeshSet<ndim>::~MeshSet() {
850       for (size_t i = 0; i < meshes.size(); ++i) {
851         delete meshes[i];
852       }
853     }
854
855
856
857     template<unsigned ndim>
858     template<typename face_type>
859     MeshSet<ndim>::FaceIter<face_type>::FaceIter(const MeshSet<ndim> *_obj, size_t _mesh, size_t _face) : obj(_obj), mesh(_mesh), face(_face) {
860     }
861
862
863
864     template<unsigned ndim>
865     template<typename face_type>
866     void MeshSet<ndim>::FaceIter<face_type>::fwd(size_t n) {
867       if (mesh < obj->meshes.size()) {
868         face += n;
869         while (face >= obj->meshes[mesh]->faces.size()) {
870           face -= obj->meshes[mesh++]->faces.size();
871           if (mesh == obj->meshes.size()) { face = 0; break; }
872         }
873       }
874     }
875
876
877
878     template<unsigned ndim>
879     template<typename face_type>
880     void MeshSet<ndim>::FaceIter<face_type>::rev(size_t n) {
881       while (n > face) {
882         n -= face;
883         if (mesh == 0) { face = 0; return; }
884         face = obj->meshes[--mesh]->faces.size() - 1;
885       }
886       face -= n;
887     }
888
889
890
891     template<unsigned ndim>
892     template<typename face_type>
893     void MeshSet<ndim>::FaceIter<face_type>::adv(int n) {
894       if (n > 0) {
895         fwd((size_t)n);
896       } else if (n < 0) {
897         rev((size_t)-n);
898       }
899     }
900
901
902
903     template<unsigned ndim>
904     template<typename face_type>
905     typename MeshSet<ndim>::template FaceIter<face_type>::difference_type
906     MeshSet<ndim>::FaceIter<face_type>::operator-(const FaceIter &other) const {
907       CARVE_ASSERT(obj == other.obj);
908       if (mesh == other.mesh) return face - other.face;
909
910       size_t m = 0;
911       for (size_t i = std::min(mesh, other.mesh) + 1; i < std::max(mesh, other.mesh); ++i) {
912         m += obj->meshes[i]->faces.size();
913       }
914
915       if (mesh < other.mesh) {
916         return -(difference_type)((obj->meshes[mesh]->faces.size() - face) + m + other.face);
917       } else {
918         return +(difference_type)((obj->meshes[other.mesh]->faces.size() - other.face) + m + face);
919       }
920     }
921
922
923
924     template<typename order_t>
925     struct VPtrSort {
926       order_t order;
927
928       VPtrSort(const order_t &_order = order_t()) : order(_order) {}
929
930       template<unsigned ndim>
931       bool operator()(carve::mesh::Vertex<ndim> *a,
932                       carve::mesh::Vertex<ndim> *b) const {
933         return order(a->v, b->v);
934       }
935     };
936
937
938
939     template<unsigned ndim>
940     void MeshSet<ndim>::collectVertices() {
941       std::unordered_map<vertex_t *, size_t> vert_idx;
942
943       for (size_t m = 0; m < meshes.size(); ++m) {
944         mesh_t *mesh = meshes[m];
945
946         for (size_t f = 0; f < mesh->faces.size(); ++f) {
947           face_t *face = mesh->faces[f];
948           edge_t *edge = face->edge;
949           do {
950             vert_idx[edge->vert] = 0;
951             edge = edge->next;
952           } while (edge != face->edge);
953         }
954       }
955
956       std::vector<vertex_t> new_vertex_storage;
957       new_vertex_storage.reserve(vert_idx.size());
958       for (typename std::unordered_map<vertex_t *, size_t>::iterator
959              i = vert_idx.begin(); i != vert_idx.end(); ++i) {
960         (*i).second = new_vertex_storage.size();
961         new_vertex_storage.push_back(*(*i).first);
962       }
963
964       for (size_t m = 0; m < meshes.size(); ++m) {
965         mesh_t *mesh = meshes[m];
966         for (size_t f = 0; f < mesh->faces.size(); ++f) {
967           face_t *face = mesh->faces[f];
968           edge_t *edge = face->edge;
969           do {
970             size_t i = vert_idx[edge->vert];
971             edge->vert = &new_vertex_storage[i];
972             edge = edge->next;
973           } while (edge != face->edge);
974         }
975       }
976
977       std::swap(vertex_storage, new_vertex_storage);
978     }
979
980
981
982     template<unsigned ndim>
983     void MeshSet<ndim>::canonicalize() {
984       std::vector<vertex_t *> vptr;
985       std::vector<vertex_t *> vmap;
986       std::vector<vertex_t> vout;
987       const size_t N = vertex_storage.size();
988
989       vptr.reserve(N);
990       vout.reserve(N);
991       vmap.resize(N);
992
993       for (size_t i = 0; i != N; ++i) {
994         vptr.push_back(&vertex_storage[i]);
995       }
996       std::sort(vptr.begin(), vptr.end(), VPtrSort<std::less<typename vertex_t::vector_t> >());
997
998       for (size_t i = 0; i != N; ++i) {
999         vout.push_back(*vptr[i]);
1000         vmap[vptr[i] - &vertex_storage[0]] = &vout[i];
1001       }
1002
1003       for (face_iter i = faceBegin(); i != faceEnd(); ++i) {
1004         for (typename face_t::edge_iter_t j = (*i)->begin(); j != (*i)->end(); ++j) {
1005           (*j).vert = vmap[(*j).vert - &vertex_storage[0]];
1006         }
1007         (*i)->canonicalize();
1008       }
1009
1010       vertex_storage.swap(vout);
1011     }
1012
1013   }
1014 }