Subversion Repositories gelsvn

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
660 khor 1
/* ----------------------------------------------------------------------- *
2
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
3
 * Copyright (C) the authors and DTU Informatics
4
 * For license and list of authors, see ../../doc/intro.pdf
5
 * ----------------------------------------------------------------------- */
6
#include <iterator>
7
#include "Manifold.h"
8
 
9
#include <iostream>
10
#include <vector>
11
#include <map>
12
#include <iterator>
13
 
14
#include "../Geometry/TriMesh.h"
15
 
16
namespace HMesh
17
{
18
 
19
    using namespace std;
20
    using namespace Geometry;
21
    using namespace CGLA;
22
 
23
    namespace
24
    {
25
        /************************************************************
26
		 * Edgekeys and Edges for halfedge identification during build
27
		 ************************************************************/
28
        class EdgeKey
29
        {
30
        public:
31
            EdgeKey(VertexID va, VertexID vb)
32
            {
33
                if(va < vb){
34
                    v0 = va;
35
                    v1 = vb;
36
                }
37
                else{
38
                    v0 = vb;
39
                    v1 = va;
40
                }
41
            }
42
 
43
            bool operator<(const EdgeKey& k2) const
44
            {
45
                if(v0 < k2.v0){
46
                    return true;
47
                }
48
                else if( k2.v0 < v0){
49
                    return false;
50
                }
51
                else{
52
                    return v1 < k2.v1;
53
                }
54
            }
55
        private:
56
            VertexID v0;
57
            VertexID v1;
58
        };
59
 
60
        struct Edge
61
        {
62
            HalfEdgeID h0;
63
            HalfEdgeID h1;
64
            int count;
65
            Edge() : count(0){}
66
        };
67
    }
68
 
69
    /*********************************************
70
	 * Public functions
71
	 *********************************************/
72
    void Manifold::build(const TriMesh& mesh)
73
    {
74
        // A vector of 3's - used to tell build how many indices each face consists of
75
        vector<int> faces(mesh.geometry.no_faces(), 3);
76
 
77
        build_template( static_cast<size_t>(mesh.geometry.no_vertices()),
78
					   reinterpret_cast<const float*>(&mesh.geometry.vertex(0)),
79
					   static_cast<size_t>(faces.size()),
80
					   static_cast<const int*>(&faces[0]),
81
					   reinterpret_cast<const int*>(&mesh.geometry.face(0)));
82
    }
83
 
84
    void Manifold::build(   size_t no_vertices,
85
						 const float* vertvec,
86
						 size_t no_faces,
87
						 const int* facevec,
88
						 const int* indices)
89
    {
90
        build_template(no_vertices, vertvec, no_faces, facevec, indices);
91
    }
92
 
93
    void Manifold::build(   size_t no_vertices,
94
						 const double* vertvec,
95
						 size_t no_faces,
96
						 const int* facevec,
97
						 const int* indices)
98
    {
99
        build_template(no_vertices, vertvec, no_faces, facevec, indices);
100
    }
101
 
102
    FaceID Manifold::add_face(std::vector<Manifold::Vec> points)
103
    {
104
        int F = points.size();
105
        vector<int> indices;
106
        for(size_t i=0;i<points.size(); ++i)
107
            indices.push_back(i);
108
        FaceID fid = *faces_end();
109
        build(points.size(), reinterpret_cast<double*>(&points[0]), 1, &F, &indices[0]);
110
        return fid;
111
    }
112
 
113
    bool Manifold::remove_face(FaceID fid)
114
    {
115
        if(!in_use(fid))
116
            return false;
117
 
118
        HalfEdgeID he = kernel.last(fid);
119
        HalfEdgeID last = he;
120
 
121
        vector<HalfEdgeID> halfedges;
122
        do
123
        {
124
            halfedges.push_back(he);
125
            kernel.set_face(he, InvalidFaceID);
126
            he = kernel.next(he);
127
        }
128
        while(he != last);
129
 
130
        vector<HalfEdgeID> halfedges_garbage;
131
        vector<VertexID> vertices_garbage;
132
        for(size_t i=0;i<halfedges.size(); ++i)
133
        {
134
            HalfEdgeID h = halfedges[i];
135
            HalfEdgeID hopp = kernel.opp(h);
136
            if(kernel.face(hopp) == InvalidFaceID)
137
            {
138
                halfedges_garbage.push_back(h);
139
                VertexID v0 = kernel.vert(hopp);
140
                if(valency(*this, v0) <= 2 )
141
                    vertices_garbage.push_back(v0);
142
                else {
143
                    link(kernel.prev(h), kernel.next(hopp));
144
                    kernel.set_out(v0, kernel.opp(kernel.prev(h)));
145
                }
146
                VertexID v1 = kernel.vert(h);
147
                if(valency(*this, v1)>2)
148
                {
149
                    link(kernel.prev(hopp),kernel.next(h));
150
                    kernel.set_out(v1, kernel.opp(kernel.prev(hopp)));
151
                }
152
 
153
            }
154
        }
155
        for(size_t i=0;i<halfedges_garbage.size();++i){
156
            kernel.remove_halfedge(kernel.opp(halfedges_garbage[i]));
157
            kernel.remove_halfedge(halfedges_garbage[i]);
158
        }
159
        for(size_t i=0;i<vertices_garbage.size(); ++i)
160
            kernel.remove_vertex(vertices_garbage[i]);
161
 
162
        kernel.remove_face(fid);
163
        return true;
164
    }
165
 
166
 
167
    bool Manifold::remove_edge(HalfEdgeID hid)
168
    {
169
        if(!in_use(hid))
170
            return false;
171
 
172
        FaceID f0 = kernel.face(hid);
173
        FaceID f1 = kernel.face(kernel.opp(hid));
174
 
175
        remove_face(f0);
176
        remove_face(f1);
177
 
178
        return true;
179
    }
180
 
181
 
182
 
183
    bool Manifold::remove_vertex(VertexID vid)
184
    {
185
        if(!in_use(vid))
186
            return false;
187
 
188
        vector<FaceID> faces;
189
        int N = circulate_vertex_ccw(*this, vid, (std::function<void(FaceID)>)[&](FaceID f) {
190
            faces.push_back(f);
191
        });
192
        for(size_t i=0;i<N;++i)
193
            remove_face(faces[i]);
194
 
195
        return true;
196
    }
197
 
198
 
199
 
200
    void Manifold::collapse_edge(HalfEdgeID h, bool avg_vertices)
201
    {
202
        HalfEdgeID ho = kernel.opp(h);
203
        VertexID hv = kernel.vert(h);
204
        VertexID hov = kernel.vert(ho);
205
        HalfEdgeID hn = kernel.next(h);
206
        HalfEdgeID hp = kernel.prev(h);
207
        HalfEdgeID hon = kernel.next(ho);
208
        HalfEdgeID hop = kernel.prev(ho);
209
		FaceID f = kernel.face(h);
210
		FaceID fo = kernel.face(ho);
211
 
212
        // average the vertex positions
213
        pos(hv) = avg_vertices ? (0.5f * (pos(hov) + pos(hv))) : pos(hv);
214
 
215
        // update all halfedges pointing to hov to point to hv, effectively removing hov from all loops
216
        HalfEdgeID he = kernel.out(hov);
217
        HalfEdgeID last = he;
218
		do {
219
			assert(kernel.vert(kernel.opp(he)) == hov);
220
            kernel.set_vert(kernel.opp(he), hv);
221
            he = kernel.next(kernel.opp(he));
222
        }
223
		while(he != last);
224
        kernel.set_out(hv, hn);
225
 
226
        // link hp and hn, effectively removing h from opposite loop
227
        link(hp, hn);
228
        // make face owning h own hn instead
229
        if(kernel.face(h) != InvalidFaceID)
230
            kernel.set_last(f, hn);
231
 
232
        // link hop and hon, effectively removing h from opposite face loop
233
        link(hop, hon);
234
        // make opposite face owning h own hon instead
235
        if(kernel.face(ho) != InvalidFaceID)
236
            kernel.set_last(fo, hon);
237
 
238
        // remove the obsolete entities
239
        kernel.remove_vertex(hov);
240
        kernel.remove_halfedge(h);
241
        kernel.remove_halfedge(ho);
242
 
243
        // verify that remaining faces haven't become degenerate because of collapse
244
        remove_face_if_degenerate(hn);
245
        remove_face_if_degenerate(hon);
246
 
247
        // verify that v is sane after collapse
248
        ensure_boundary_consistency(hv);
249
    }
250
 
251
    FaceID Manifold::split_face_by_edge(FaceID f, VertexID v0, VertexID v1)
252
    {
253
		assert(!connected(*this, v0, v1));
254
        HalfEdgeID he = kernel.last(f);
255
        HalfEdgeID last = he;
256
        int steps = 0;
257
        while(he != last || steps == 0){
258
            ++steps;
259
            he = kernel.next(he);
260
        }
261
 
262
        // make sure this is not a triangle
263
        assert(steps > 3);
264
        // make sure we are not trying to connect a vertex to itself
265
        assert(v0 != v1);
266
 
267
		HalfEdgeID h0 = kernel.out(v0);
268
		for(Walker w = walker(v0); !w.full_circle(); w = w.circulate_vertex_cw()){
269
			if(w.face() == f){
270
				h0 = w.halfedge();
271
				break;
272
			}
273
		}
274
		assert(kernel.face(h0) != InvalidFaceID);
275
        assert(kernel.face(h0) == f);
276
 
277
        // the halfedge belonging to f, going out from v0, is denoted h. Move along h until we hit v1.
278
        // now we have the halfedge which belongs to f and points to v1.
279
        HalfEdgeID h = h0;
280
        while(kernel.vert(h) != v1){
281
            h = kernel.next(h);
282
        }
283
        assert(h != h0);
284
 
285
        // create a new halfedge ha which connects v1 and v0 closing the first loop.
286
        HalfEdgeID h1 = kernel.next(h);
287
        HalfEdgeID ha = kernel.add_halfedge();
288
        link(h, ha);
289
        link(ha, h0);
290
        kernel.set_face(ha, f);
291
        kernel.set_vert(ha, v0);
292
        kernel.set_last(f, ha);
293
 
294
        // create a new face, f2, and set all halfedges in the remaining part of the polygon to point to this face.
295
        h = h1;
296
        FaceID f2 = kernel.add_face();
297
        while(kernel.vert(h) != v0){
298
            kernel.set_face(h, f2);
299
            h = kernel.next(h);
300
        }
301
        kernel.set_face(h, f2);
302
        assert(h != h1);
303
 
304
        // create a new halfedge hb to connect v0 and v1.
305
        HalfEdgeID hb = kernel.add_halfedge();
306
        link(h, hb);
307
        link(hb, h1);
308
        kernel.set_face(hb, f2);
309
        kernel.set_vert(hb, v1);
310
        kernel.set_last(f2, hb);
311
 
312
        // complete the operation by gluing the two new halfedges
313
        glue(ha, hb);
314
 
315
        // assert sanity of operation
316
        assert(kernel.next(kernel.opp(kernel.prev(h1))) == h0);
317
        assert(kernel.next(kernel.opp(kernel.prev(h0))) == h1);
318
        assert(kernel.face(kernel.next(hb)) == f2);
319
        assert(kernel.face(kernel.next(kernel.next(hb))) == f2);
320
        assert(kernel.face(hb) == f2);
321
 
322
        // return handle to newly created face
323
        return f2;
324
    }
325
 
326
    VertexID Manifold::split_face_by_vertex(FaceID f)
327
    {
328
        //create the new vertex, with the barycenter of the face as position
329
        Manifold::Vec p(0.0f);
330
        HalfEdgeID last_he = kernel.last(f);
331
        HalfEdgeID he = last_he;
332
        int steps = 0;
333
 
334
        while(he != last_he || steps == 0){
335
            p += positions[kernel.vert(he)];
336
            ++steps;
337
            he = kernel.next(he);
338
        }
339
 
340
        p /= steps;
341
 
342
        VertexID v = kernel.add_vertex();
343
        positions[v]  = p;
344
 
345
        //circulate the face, create halfedges and connect to vertex
346
        vector<HalfEdgeID> eout;
347
        vector<HalfEdgeID> ein;
348
        last_he = kernel.last(f);
349
        he = last_he;
350
 
351
        do{
352
            HalfEdgeID hn = kernel.next(he);
353
 
354
            HalfEdgeID ho = kernel.add_halfedge();
355
            HalfEdgeID hi = kernel.add_halfedge();
356
 
357
            FaceID fn = kernel.add_face();
358
            kernel.set_face(hi, fn);
359
            kernel.set_vert(hi, v);
360
            kernel.set_face(ho, fn);
361
            kernel.set_vert(ho, kernel.vert(kernel.opp(he)));
362
            kernel.set_face(he, fn);
363
            kernel.set_last(fn, ho);
364
 
365
            link(hi, ho);
366
            link(ho, he);
367
            link(he, hi);
368
 
369
            eout.push_back(ho);
370
            ein.push_back(hi);
371
 
372
            he = hn;
373
        }
374
        while(he != last_he);
375
 
376
        // glue new halfedges together
377
        size_t N = eout.size();
378
        for(size_t i = 0; i < N; ++i){
379
            glue(ein[i], eout[(i+1)%N]);
380
        }
381
        kernel.set_out(v, eout[0]);
382
 
383
        //remove the now replaced face
384
        kernel.remove_face(f);
385
 
386
        return v;
387
    }
388
    VertexID Manifold::split_edge(HalfEdgeID h)
389
    {
390
        HalfEdgeID ho = kernel.opp(h);
391
        VertexID v = kernel.vert(h);
392
        VertexID vo = kernel.vert(ho);
393
 
394
        //create the new vertex with middle of edge as position and update connectivity
395
        VertexID vn = kernel.add_vertex();
396
        positions[vn] = .5f * (positions[v] + positions[vo]);
397
        kernel.set_out(vn, h);
398
 
399
        //create two necessary halfedges, and update connectivity
400
        HalfEdgeID hn = kernel.add_halfedge();
401
        HalfEdgeID hno = kernel.add_halfedge();
402
 
403
        kernel.set_out(vo, hn);
404
        kernel.set_out(v, hno);
405
 
406
        glue(h, hno);
407
        glue(hn, ho);
408
 
409
        link(kernel.prev(h), hn);
410
        link(hn, h);
411
        kernel.set_vert(hn, vn);
412
 
413
        link(kernel.prev(ho), hno);
414
        link(hno, ho);
415
        kernel.set_vert(hno, vn);
416
 
417
        // update faces in case of non boundary edge
418
        if(kernel.face(h) != InvalidFaceID)
419
            kernel.set_last(kernel.face(h), hn);
420
 
421
        if(kernel.face(ho) != InvalidFaceID)
422
            kernel.set_last(kernel.face(ho), ho);
423
 
424
        // update new edge with faces from existing edge
425
        kernel.set_face(hn, kernel.face(h));
426
        kernel.set_face(hno, kernel.face(ho));
427
 
428
        //if split occurs on a boundary, consistency must be ensured.
429
        ensure_boundary_consistency(vn);
430
        ensure_boundary_consistency(v);
431
        ensure_boundary_consistency(vo);
432
 
433
        return vn;
434
    }
435
 
436
    size_t link_intersection(const Manifold& m, VertexID v0, VertexID v1, vector<VertexID>& lisect)
437
    {
438
        // get the one-ring of v0
439
        vector<VertexID> link0;
440
        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID vn) {
441
            link0.push_back(vn);
442
        });
443
 
444
        // get the one-ring of v1
445
        vector<VertexID> link1;
446
        circulate_vertex_ccw(m, v1, (std::function<void(VertexID)>)[&](VertexID vn) {
447
            link1.push_back(vn);
448
        });
449
 
450
        // sort the vertices of the two rings
451
        sort(link0.begin(), link0.end());
452
        sort(link1.begin(), link1.end());
453
 
454
        // get the intersection of the shared vertices from both rings
455
        std::back_insert_iterator<vector<VertexID> > lii(lisect);
456
        set_intersection(link0.begin(), link0.end(),
457
						 link1.begin(), link1.end(),
458
						 lii);
459
 
460
        return lisect.size();
461
    }
462
 
463
    bool Manifold::stitch_boundary_edges(HalfEdgeID h0, HalfEdgeID h1)
464
    {
465
        // Cannot stitch an edge with itself
466
        if(h0 == h1)
467
            return false;
468
 
469
        // Only stitch a pair of boundary edges.
470
        if(kernel.face(h0) == InvalidFaceID && kernel.face(h1) == InvalidFaceID)
471
        {
472
            HalfEdgeID h0o = kernel.opp(h0);
473
            HalfEdgeID h1o = kernel.opp(h1);
474
            VertexID v0a = kernel.vert(h0);
475
            VertexID v0b = kernel.vert(kernel.opp(h1));
476
            VertexID v1b = kernel.vert(h1);
477
            VertexID v1a = kernel.vert(kernel.opp(h0));
478
 
479
            // Case below implies that h0 and h1 are the same edge with different ID
480
            // That should not happen.
481
            assert(!(v0a == v0b && v1a == v1b));
482
 
483
            if(connected(*this, v0a, v0b))
484
                return false;
485
            if(connected(*this, v1a, v1b))
486
                return false;
487
 
488
 
489
            if(v0b != v0a)
490
            {
491
 
492
                // Check the link intersection v0a and v0b are welded together
493
                // if they share a neighbouring vertex, it will appear twice in the combined
494
                // one ring unless it v1a and v1a==v1b
495
                vector<VertexID> lisect;
496
                if(link_intersection(*this, v0a, v0b, lisect))
497
                {
498
                    vector<VertexID>::iterator iter;
499
 
500
                    if(v1a == v1b)
501
                    {
502
                        iter = find(lisect.begin(), lisect.end(), v1a);
503
                        if(iter != lisect.end())
504
                            lisect.erase(iter);
505
                    }
506
                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h0)));
507
                    if(iter != lisect.end())
508
                        lisect.erase(iter);
509
                    if(lisect.size() > 0)
510
                        return false;
511
                }
512
            }
513
 
514
            if(v1a != v1b)
515
            {
516
                // Check the same for the other endpoints.
517
                vector<VertexID> lisect;
518
                if(link_intersection(*this, v1a, v1b, lisect))
519
                {
520
                    vector<VertexID>::iterator iter;
521
 
522
                    if(v0a == v0b)
523
                    {
524
                        iter = find(lisect.begin(), lisect.end(), v1a);
525
                        if(iter != lisect.end())
526
                            lisect.erase(iter);
527
                    }
528
                    iter = find(lisect.begin(), lisect.end(), kernel.vert(kernel.next(h1)));
529
                    if(iter != lisect.end())
530
                        lisect.erase(iter);
531
                    if(lisect.size() > 0)
532
                        return false;
533
                }
534
            }
535
 
536
 
537
            if(v0b != v0a)
538
                circulate_vertex_ccw(*this, v0b, (std::function<void(Walker&)>)[&](Walker hew) {
539
                    kernel.set_vert(hew.opp().halfedge(), v0a);
540
                });
541
 
542
            if(v1b != v1a)
543
                circulate_vertex_ccw(*this, v1b, (std::function<void(Walker&)>)[&](Walker hew) {
544
                    kernel.set_vert(hew.opp().halfedge(), v1a);
545
                });
546
 
547
            if(v0a != v0b)
548
            {
549
                HalfEdgeID h1p = kernel.prev(h1);
550
                HalfEdgeID h0n = kernel.next(h0);
551
 
552
                if(kernel.next(h0n) == h1p)
553
                {
554
                    glue(kernel.opp(h0n), kernel.opp(h1p));
555
                    kernel.set_out(kernel.vert(h0n),kernel.opp(h0n));
556
                    kernel.remove_halfedge(h0n);
557
                    kernel.remove_halfedge(h1p);
558
                }
559
                else
560
                    link(h1p, h0n);
561
                kernel.remove_vertex(v0b);
562
            }
563
 
564
            if(v1a != v1b)
565
            {
566
                HalfEdgeID h0p = kernel.prev(h0);
567
                HalfEdgeID h1n = kernel.next(h1);
568
                if(kernel.next(h1n) == h0p)
569
                {
570
                    glue(kernel.opp(h0p), kernel.opp(h1n));
571
                    kernel.set_out(kernel.vert(h1n), kernel.opp(h1n));
572
                    kernel.remove_halfedge(h0p);
573
                    kernel.remove_halfedge(h1n);
574
                }
575
                else
576
                    link(h0p, h1n);
577
                kernel.remove_vertex(v1b);
578
            }
579
            glue(h0o, h1o);
580
 
581
            kernel.remove_halfedge(h0);
582
            kernel.remove_halfedge(h1);
583
 
584
            kernel.set_out(v1a, h1o);
585
            kernel.set_out(v0a, h0o);
586
 
587
            ensure_boundary_consistency(v1a);
588
            ensure_boundary_consistency(v0a);
589
 
590
            return true;
591
        }
592
        return false;
593
    }
594
 
595
 
596
 
597
    FaceID Manifold::merge_one_ring(VertexID v, float max_loop_length)
598
    {
599
        // If the vertex is either not in use or has just
600
        // one incident edge (or less), bail out.
601
        int val = valency(*this,v);
602
        if(!in_use(v) || val<2)
603
            return InvalidFaceID;
604
 
605
        // If the vertex is  a boundary vertex, we preparte the walker so that the
606
        // first face visited is not the invalid face outside the boundary. If the boundary
607
        // vertex is adjacent to only one vertex, there is little to do and we bail out.
608
        bool vertex_is_boundary = false;
609
        Walker hew = walker(v);
610
        if(boundary(*this, v))
611
        {
612
            if(val==2) return InvalidFaceID;
613
            vertex_is_boundary = true;
614
            hew = hew.circulate_vertex_ccw();
615
        }
616
 
617
        // Prepare some vectors for taking out the trash: We remove all old faces and all orphaned edges
618
        // and vertices
619
        vector<HalfEdgeID> garbage_halfedges;
620
        vector<FaceID> garbage_faces;
621
        vector<VertexID> garbage_vertices;
622
 
623
        vector<HalfEdgeID> loop; // The halfedges which form the outer loop of the merged one ring.
624
 
625
        // Below we loop over all faces adjacent to the vertex and add their halfedges to a big loop
626
        // which will form the loop of the new merged face. Below we remove from the loop edges
627
        // that appear twice (as each other's opposite).
628
        for(;!hew.full_circle(); hew = hew.circulate_vertex_ccw())
629
        {
630
            garbage_faces.push_back(hew.face());
631
            for(Walker hew2 = walker(hew.halfedge());
632
                !hew2.full_circle(); hew2 = hew2.circulate_face_ccw())
633
                loop.push_back(hew2.halfedge());
634
        }
635
 
636
 
637
        // Now merge the loops. We iteratively remove pairs of adjacent halfedges from the loop
638
        // if we find that the second is the opposite of the first since this is a degenerate
639
        // situation. However, we stop at four remaining halfedges since otherwise the loop degenerates
640
        // to zero after the next step if these four are also pairwise each other's opposites.
641
        int did_work;
642
        do
643
        {
644
            did_work = 0;
645
            vector<HalfEdgeID> loop_tmp(0);
646
            for(size_t i=0;i<loop.size();++i)
647
                if(walker(loop[i]).opp().halfedge() == loop[(i+1)%loop.size()])
648
                {
649
                    VertexID vid = walker(loop[i]).vertex();
650
                    if(vid != v)
651
                        garbage_vertices.push_back(walker(loop[i]).vertex());
652
                    garbage_halfedges.push_back(loop[i]);
653
                    garbage_halfedges.push_back(loop[(i+1)%loop.size()]);
654
                    loop[i] = InvalidHalfEdgeID;
655
                    loop[(i+1)%loop.size()] = InvalidHalfEdgeID;
656
                    ++did_work;
657
                    ++i;
658
                }
659
            for(size_t i=0;i<loop.size();++i)
660
                if(loop[i] != InvalidHalfEdgeID)
661
                    loop_tmp.push_back(loop[i]);
662
            loop = loop_tmp;
663
        } while(did_work > 0 && loop.size() > 4);
664
 
665
        // Check whether the loop is too long
666
        float loop_len=0.0;
667
        for(size_t i=0;i<loop.size();++i)
668
            loop_len += length(*this, loop[i]);
669
        if(loop_len > max_loop_length)
670
            return InvalidFaceID;
671
 
672
        // The following block checks wheteher the same halfedge appears twice. In this
673
        // case we fail since it means that the one ring contains the same face twice.
674
        vector<HalfEdgeID> loop_tmp = loop;
675
        sort(loop_tmp.begin(), loop_tmp.end());
676
        vector<HalfEdgeID>::iterator end_iter = unique(loop_tmp.begin(), loop_tmp.end());
677
        if(end_iter != loop_tmp.end())
678
            return InvalidFaceID;
679
 
680
        // Remove all faces and connected halfedges and the original vertex v.
681
        for(size_t i=0;i<garbage_vertices.size(); ++i)
682
            kernel.remove_vertex(garbage_vertices[i]);
683
        for(size_t i=0;i<garbage_faces.size(); ++i)
684
            kernel.remove_face(garbage_faces[i]);
685
        for(size_t i=0;i<garbage_halfedges.size(); ++i)
686
            kernel.remove_halfedge(garbage_halfedges[i]);
687
        if(!vertex_is_boundary)
688
            kernel.remove_vertex(v);
689
 
690
        // Create a new face for the merged one ring and link up all the halfedges
691
        // in the loop
692
        FaceID f = kernel.add_face();
693
        kernel.set_last(f,loop[0]);
694
        for(size_t i=0;i<loop.size(); ++i)
695
        {
696
            kernel.set_face(loop[i], f);
697
            Walker hw = walker(loop[i]);
698
            kernel.set_out(hw.vertex(),hw.opp().halfedge());
699
            link(loop[i],loop[(i+1)%loop.size()]);
700
            assert(hw.vertex() == walker(loop[(i+1)%loop.size()]).opp().vertex());
701
        }
702
 
703
        // Finally, we ensure boundary consitency for all vertices in the loop.
704
        for(size_t i=0;i<loop.size(); ++i)
705
        {
706
            Walker hw = walker(loop[i]);
707
            ensure_boundary_consistency(hw.vertex());
708
        }
709
 
710
        // Return the brand new merged face.
711
        return f;
712
    }
713
 
714
 
715
 
716
    bool Manifold::merge_faces(FaceID f, HalfEdgeID h)
717
    {
718
        //assert that we're merging a valid face with the corresponding halfedge
719
        assert(kernel.face(h) == f);
720
 
721
        HalfEdgeID ho = kernel.opp(h);
722
        FaceID fo = kernel.face(ho);
723
        HalfEdgeID hn = kernel.next(h);
724
        HalfEdgeID hon = kernel.next(ho);
725
 
726
        if(fo == f)
727
            return false;
728
 
729
        //boundary case
730
        if(kernel.vert(hn) == kernel.vert(hon))
731
            return false;
732
 
733
        HalfEdgeID hp = kernel.prev(h);
734
        HalfEdgeID hop = kernel.prev(ho);
735
        VertexID v = kernel.vert(h);
736
        VertexID vo = kernel.vert(ho);
737
 
738
        if(valency(*this, v) < 3 || valency(*this, vo) < 3)
739
            return false;
740
 
741
        link(hop, hn);
742
        link(hp, hon);
743
        kernel.set_out(vo, hon);
744
        kernel.set_out(v, hn);
745
        kernel.set_last(f, hn);
746
 
747
        HalfEdgeID hx = hon;
748
 
749
        assert(kernel.face(hx) == fo);
750
        while(kernel.face(hx) != f){
751
            kernel.set_face(hx, f);
752
            hx = kernel.next(hx);
753
        }
754
 
755
        ensure_boundary_consistency(v);
756
        ensure_boundary_consistency(vo);
757
 
758
        kernel.remove_halfedge(h);
759
        kernel.remove_halfedge(ho);
760
        kernel.remove_face(fo);
761
 
762
        return true;
763
    }
764
 
765
    FaceID Manifold::close_hole(HalfEdgeID h)
766
    {
767
        // invalid face is a hole
768
        if(kernel.face(h) == InvalidFaceID){
769
            FaceID f = kernel.add_face();
770
            kernel.set_last(f, h);
771
            do{
772
                kernel.set_face(h, f);
773
                h = kernel.next(h);
774
            }
775
            while(kernel.face(h) != f);
776
            return f;
777
        }
778
        return kernel.face(h);
779
    }
780
 
781
    VertexID Manifold::slit_vertex(VertexID v, HalfEdgeID h_in, HalfEdgeID h_out)
782
    {
783
        assert(kernel.face(h_in) != InvalidFaceID);
784
        assert(kernel.face(h_out) != InvalidFaceID);
785
        assert(kernel.opp(h_out) != h_in);
786
 
787
        VertexID v_new = kernel.add_vertex();
788
        pos(v_new) = pos(v);
789
        HalfEdgeID h = kernel.prev(h_out);
790
        kernel.set_vert(h, v_new);
791
        while ( h != h_in) {
792
            h = kernel.prev(kernel.opp(h));
793
            kernel.set_vert(h, v_new);
794
        }
795
 
796
        HalfEdgeID h_in_opp = kernel.opp(h_in);
797
        HalfEdgeID hn_in, hn_in_opp;
798
        if(kernel.face(h_in_opp) != InvalidFaceID)
799
        {
800
            hn_in = kernel.add_halfedge();
801
            kernel.set_face(hn_in, InvalidFaceID);
802
            glue(h_in_opp, hn_in);
803
 
804
            hn_in_opp = kernel.add_halfedge();
805
            kernel.set_face(hn_in_opp, InvalidFaceID);
806
            glue(h_in, hn_in_opp);
807
 
808
            link(hn_in_opp, hn_in);
809
 
810
            VertexID v_i = kernel.vert(h_in_opp);
811
            kernel.set_vert(hn_in_opp, v_i);
812
            kernel.set_out(v_i, hn_in);
813
        }
814
        else
815
        {
816
            hn_in_opp = h_in_opp;
817
            hn_in = kernel.prev(hn_in_opp);
818
            h_in_opp = kernel.opp(hn_in);
819
        }
820
 
821
        HalfEdgeID h_out_opp = kernel.opp(h_out);
822
        HalfEdgeID hn_out,hn_out_opp;
823
        if(kernel.face(h_out_opp) != InvalidFaceID)
824
        {
825
            hn_out = kernel.add_halfedge();
826
            kernel.set_face(hn_out, InvalidFaceID);
827
            glue(h_out_opp, hn_out);
828
 
829
            hn_out_opp = kernel.add_halfedge();
830
            kernel.set_face(hn_out_opp, InvalidFaceID);
831
            glue(h_out, hn_out_opp);
832
 
833
            link(hn_out, hn_out_opp);
834
 
835
            VertexID v_o = kernel.vert(h_out);
836
            kernel.set_vert(hn_out, v_o);
837
            kernel.set_out(v_o, hn_out_opp);
838
        }
839
        else
840
        {
841
            hn_out_opp = h_out_opp;
842
            hn_out = kernel.next(hn_out_opp);
843
            h_out_opp = kernel.opp(hn_out);
844
        }
845
 
846
        link(hn_out_opp, hn_in_opp);
847
        link(hn_in, hn_out);
848
 
849
        kernel.set_vert(hn_in, v);
850
        kernel.set_vert(hn_out_opp, v_new);
851
 
852
        kernel.set_out(v, hn_out);
853
        kernel.set_out(v_new, hn_in_opp);
854
 
855
        return v_new;
856
    }
857
 
858
 
859
    HalfEdgeID Manifold::slit_edges(VertexAttributeVector<int>& insel)
860
    {
861
        HalfEdgeID h;
862
        for(auto vid : vertices())
863
        {
864
            if(insel[vid])
865
            {
866
                HalfEdgeID h_in = InvalidHalfEdgeID, h_out = InvalidHalfEdgeID;
867
                Walker w = walker(vid);
868
                while(!w.full_circle())
869
                {
870
                    if(insel[w.vertex()]) {
871
                        if(h_in == InvalidHalfEdgeID) {
872
                            if(w.opp().face() == InvalidFaceID)
873
                                h_in = w.opp().next().opp().halfedge();
874
                            else
875
                                h_in = w.opp().halfedge();
876
                        }
877
                        else {
878
                            if(w.face() == InvalidFaceID)
879
                                h_out = w.prev().opp().halfedge();
880
                            else
881
                                h_out = w.halfedge();
882
                            break;
883
                        }
884
                    }
885
                    w = w.circulate_vertex_ccw();
886
                }
887
                if(h_in != InvalidHalfEdgeID &&
888
                   h_out != InvalidHalfEdgeID) {
889
                    VertexID v_new = slit_vertex(vid, h_in, h_out);
890
                    h = walker(v_new).halfedge();
891
                }
892
            }
893
        }
894
        return h;
895
    }
896
 
897
 
898
    void Manifold::flip_edge(HalfEdgeID h)
899
    {
900
        HalfEdgeID hn = kernel.next(h);
901
        HalfEdgeID hp = kernel.prev(h);
902
        HalfEdgeID ho = kernel.opp(h);
903
        HalfEdgeID hon = kernel.next(ho);
904
        HalfEdgeID hop = kernel.prev(ho);
905
 
906
        FaceID hf = kernel.face(h);
907
        FaceID hof = kernel.face(ho);
908
 
909
        VertexID hv = kernel.vert(h);
910
        VertexID hnv = kernel.vert(hn);
911
        VertexID hov = kernel.vert(ho);
912
        VertexID honv = kernel.vert(hon);
913
 
914
        // update face connectivity of halfedges changing face
915
        kernel.set_face(hop, hf);
916
        kernel.set_face(hp, hof);
917
 
918
        // connectivity of faces with halfedges of flipped edge
919
        kernel.set_last(hf, h);
920
        kernel.set_last(hof, ho);
921
 
922
        // connectivity of halfedges of first face after flip
923
        link(hn, h);
924
        link(h, hop);
925
        link(hop, hn);
926
 
927
        // connectivity of halfedges of second face after flip
928
        link(hon, ho);
929
        link(ho, hp);
930
        link(hp, hon);
931
 
932
        // connectivity of flipped edge and corresponding vertices
933
        kernel.set_vert(h, honv);
934
        kernel.set_vert(ho, hnv);
935
 
936
        if(kernel.out(hv) == ho)
937
            kernel.set_out(hv, hn);
938
        if(kernel.out(hov) == h)
939
            kernel.set_out(hov, hon);
940
 
941
        //
942
        //        // if the flip occurs next to a boundary, ensure the boundary consistency
943
        //        ensure_boundary_consistency(hv);
944
        //        ensure_boundary_consistency(hov);
945
    }
946
 
947
 
948
    /**********************************************
949
     * Private functions
950
     **********************************************/
951
 
952
    template<typename size_type, typename float_type, typename int_type>
953
    void Manifold::build_template(  size_type no_vertices,
954
                                  const float_type* vertvec,
955
                                  size_type no_faces,
956
                                  const int_type* facevec,
957
                                  const int_type* indices)
958
    {
959
        vector<VertexID> vids(no_vertices);
960
 
961
        // create vertices corresponding to positions stored in vertvec
962
        for(size_t i=0;i<no_vertices;++i)
963
        {
964
            const float_type* v0 = &vertvec[i*3];
965
            pos(vids[i] = kernel.add_vertex()) = Manifold::Vec(v0[0], v0[1], v0[2]);
966
        }
967
 
968
        //map over the created edges - needed to preserve edge uniqueness
969
        typedef map<EdgeKey, Edge> EdgeMap;
970
        EdgeMap edge_map;
971
 
972
        // counter that jumps between faces in indices vector
973
        int_type n  = 0;
974
 
975
        // create faces correspponding to faces stored in facevec
976
        for(size_type i = 0; i < no_faces; ++i){
977
            //amount of vertices in current face
978
            size_type N = facevec[i];
979
            vector<HalfEdgeID> fh;
980
 
981
            //each face indice results corresponds to 1 edge, 2 halfedges
982
            for(size_type j = 0; j < N; ++j){
983
                // ensure indice integrity
984
 
985
                assert(static_cast<size_type>(indices[j + n]) < no_vertices);
986
                assert(static_cast<size_type>(indices[(j + 1) % N + n]) < no_vertices);
987
 
988
 
989
                // each iteration uses two indices from the face
990
                VertexID v0 = vids[static_cast<size_type>(indices[j + n])];
991
                VertexID v1 = vids[static_cast<size_type>(indices[(j + 1) % N + n])];
992
 
993
                // create key and search map for edge
994
                EdgeKey ek(v0, v1);
995
                EdgeMap::iterator em_iter = edge_map.find(ek);
996
 
997
                // current edge has not been created
998
                if(em_iter == edge_map.end()){
999
                    // create edge for map
1000
                    Edge e;
1001
                    e.h0 = kernel.add_halfedge();
1002
                    e.h1 = kernel.add_halfedge();
1003
                    e.count = 1;
1004
 
1005
                    // glue operation: 1 edge = 2 glued halfedges
1006
                    glue(e.h0, e.h1);
1007
 
1008
                    // update vertices with their outgoing halfedges
1009
                    kernel.set_out(v0, e.h0);
1010
                    kernel.set_out(v1, e.h1);
1011
 
1012
                    // update halfedges with the vertices they point to
1013
                    kernel.set_vert(e.h0, v1);
1014
                    kernel.set_vert(e.h1, v0);
1015
 
1016
                    // update map
1017
                    edge_map[ek] = e;
1018
 
1019
                    // update vector of halfedges belonging to current face
1020
                    fh.push_back(e.h0);
1021
                }
1022
                else{
1023
                    // update current face with halfedge from edge
1024
                    fh.push_back(em_iter->second.h1);
1025
 
1026
                    // asserting that a halfedge is visited exactly twice;
1027
                    // once for each face on either side of the edge.
1028
                    em_iter->second.count++;
1029
                    assert(em_iter->second.count == 2);
1030
                }
1031
            }
1032
 
1033
            FaceID fid = kernel.add_face();
1034
            for(size_type j = 0; j < N; ++j){
1035
                // update halfedge with face
1036
                kernel.set_face(fh[j], fid);
1037
 
1038
                // link operation: link two halfedges in the same face
1039
                link(fh[j], fh[(j + 1) % N]);
1040
            }
1041
            //update face with the first halfedge created
1042
            kernel.set_last(fid, fh[0]);
1043
 
1044
            // step to first index of next face
1045
            n += N;
1046
        }
1047
 
1048
        // test for unused vertices
1049
        for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
1050
            assert( (*v) != InvalidVertexID);
1051
            if(kernel.out(*v) == InvalidHalfEdgeID)
1052
                kernel.remove_vertex(*v);
1053
        }
1054
 
1055
        // boundary check while avoiding unused vertices
1056
        for(VertexIDIterator v = vertices_begin(); v != vertices_end(); ++v){
1057
            if((*v) != InvalidVertexID)
1058
                ensure_boundary_consistency(*v);
1059
        }
1060
    }
1061
    void Manifold::link(HalfEdgeID h0, HalfEdgeID h1)
1062
    {
1063
        kernel.set_next(h0, h1);
1064
        kernel.set_prev(h1, h0);
1065
    }
1066
    void Manifold::glue(HalfEdgeID h0, HalfEdgeID h1)
1067
    {
1068
        kernel.set_opp(h0, h1);
1069
        kernel.set_opp(h1, h0);
1070
    }
1071
    void Manifold::ensure_boundary_consistency(VertexID v)
1072
    {
1073
        // boundary consistency check by performing two vertex circulations
1074
        HalfEdgeID h = kernel.out(v);
1075
        HalfEdgeID last = h;
1076
 
1077
        // step 1: circle through edges pointing away from vertex until reaching a null face
1078
        while(kernel.face(h) != InvalidFaceID){
1079
            h = kernel.opp(kernel.prev(h));
1080
            if(h == last) // We came full circle - vertex not boundary - return.
1081
                return;
1082
        }
1083
        // null face encountered, we update our vertex with half edge index and prepare for step 2
1084
        kernel.set_out(v, h);
1085
        HalfEdgeID ho = kernel.opp(h);
1086
 
1087
        // step 2: circle through edges pointing towards vertex until reaching a null face
1088
        while(kernel.face(ho) != InvalidFaceID){
1089
            ho = kernel.opp(kernel.next(ho));
1090
        }
1091
        // null face encountered again, we update our edge with vertex index
1092
        kernel.set_vert(ho, v);
1093
 
1094
        // remaining step is to make the in and out going edges link to each other
1095
        link(ho, h);
1096
    }
1097
    void Manifold::remove_face_if_degenerate(HalfEdgeID h)
1098
    {
1099
        // face is degenerate if there is only two halfedges in face loop
1100
        if(kernel.next(kernel.next(h)) == h)
1101
        {
1102
            HalfEdgeID hn = kernel.next(h);
1103
            HalfEdgeID ho = kernel.opp(h);
1104
            HalfEdgeID hno = kernel.opp(hn);
1105
            VertexID hv = kernel.vert(h);
1106
            VertexID hnv = kernel.vert(hn);
1107
            FaceID f = kernel.face(h);
1108
 
1109
            assert(ho != hn);
1110
            assert(h != hno);
1111
            assert(hv != hnv);
1112
 
1113
            // glue opposite halfedge to halfedge opposite next halfedge, obsoleting h and hn from loop
1114
            glue(ho, hno);
1115
 
1116
            // make v and vn point to opposite and next opposite halfedge, obsoleting h and hn from loop
1117
            kernel.set_out(hnv, hno);
1118
            kernel.set_out(hv, ho);
1119
 
1120
            // if face owning h is valid, remove face
1121
            if(f != InvalidFaceID)
1122
                kernel.remove_face(f);
1123
            // remove the two invalid halfedges and the invalid face loop
1124
            kernel.remove_halfedge(h);
1125
            kernel.remove_halfedge(hn);
1126
 
1127
            // verify that v and vn is sane after removing the degenerate face
1128
            ensure_boundary_consistency(hv);
1129
            ensure_boundary_consistency(hnv);
1130
        }
1131
    }
1132
 
1133
    vector<HalfEdgeID> Manifold::bridge_faces(FaceID f0, FaceID f1, const vector<pair<VertexID, VertexID> >& pairs)
1134
    {
1135
        // Let N be the number of vertex pairs to connect by edges
1136
        size_t N = pairs.size();
1137
 
1138
        // We now create N pairs of edges, glue each pair and let them point to the appropriate
1139
        // vertices.
1140
        vector<HalfEdgeID> new_halfedges(N);
1141
        vector<HalfEdgeID> new_halfedges_opp(N);
1142
        for(size_t i=0;i<N; ++i)
1143
        {
1144
            new_halfedges[i] = kernel.add_halfedge();
1145
            new_halfedges_opp[i] = kernel.add_halfedge();
1146
            glue(new_halfedges[i], new_halfedges_opp[i]);
1147
            kernel.set_vert(new_halfedges[i], pairs[i].second);
1148
            kernel.set_vert(new_halfedges_opp[i], pairs[i].first);
1149
        }
1150
 
1151
        // We now maintain some halfedge indices to right before
1152
        // and right after the point we are trying to connect on
1153
        // each of the two loops.
1154
        HalfEdgeID h0p = kernel.last(f0);
1155
        HalfEdgeID h1n = kernel.last(f1);
1156
 
1157
        // loop over all connections and link the new halfedges with the old
1158
        // ones.
1159
        for(size_t i=0;i<N; ++i)
1160
        {
1161
            while(kernel.vert(h0p) != pairs[i].first)
1162
                h0p = kernel.next(h0p);
1163
            while(kernel.vert(kernel.prev(h1n)) != pairs[i].second)
1164
                h1n = kernel.prev(h1n);
1165
 
1166
            HalfEdgeID h0n = kernel.next(h0p);
1167
            HalfEdgeID h1p = kernel.prev(h1n);
1168
 
1169
            link(h0p, new_halfedges[i]);
1170
            link(new_halfedges[i],h1n);
1171
            link(h1p, new_halfedges_opp[i]);
1172
            link(new_halfedges_opp[i],h0n);
1173
 
1174
            h0p = new_halfedges_opp[i];
1175
            h1n = new_halfedges_opp[i];
1176
        }
1177
 
1178
        // Create the faces and their connections
1179
        for(size_t i=0;i<N; ++i)
1180
        {
1181
            HalfEdgeID last = new_halfedges[i];
1182
            FaceID f = kernel.add_face();
1183
            kernel.set_last(f, last);
1184
            HalfEdgeID h = last;
1185
            do
1186
            {
1187
                kernel.set_face(h, f);
1188
                h = kernel.next(h);
1189
            } while(h != last);
1190
        }
1191
 
1192
        // Delete the old faces that were bridged.
1193
        kernel.remove_face(f0);
1194
        kernel.remove_face(f1);
1195
 
1196
        new_halfedges.insert(new_halfedges.end(), new_halfedges_opp.begin(), new_halfedges_opp.end());
1197
        return new_halfedges;
1198
    }
1199
 
1200
    /***************************************************
1201
     * Namespace functions
1202
     ***************************************************/
1203
    bool valid(const Manifold& m)
1204
    {
1205
        // Verify components of halfedges
1206
        for(HalfEdgeIDIterator h = m.halfedges_begin(); h != m.halfedges_end(); ++h){
1207
            Walker j = m.walker(*h);
1208
 
1209
            if(j.vertex() == InvalidVertexID){
1210
                cout << "Halfedge lacks vert" << endl;
1211
                return false;
1212
            }
1213
            if(j.next().halfedge() == InvalidHalfEdgeID){
1214
                cout << "Halfedge lacks next" << endl;
1215
                return false;
1216
            }
1217
            if(j.prev().halfedge() == InvalidHalfEdgeID){
1218
                cout << "Halfedge lacks prev" << endl;
1219
                return false;
1220
            }
1221
            if(j.opp().halfedge() == InvalidHalfEdgeID){
1222
                cout << "Halfedge lacks opp" << endl;
1223
                return false;
1224
            }
1225
 
1226
        }
1227
        // Verify components of vertices
1228
        for(VertexIDIterator v = m.vertices_begin(); v != m.vertices_end(); ++v){
1229
            vector<VertexID> link;
1230
 
1231
            // circulate the halfedges of vertex
1232
            for(Walker j = m.walker(*v); !j.full_circle(); j = j.circulate_vertex_cw()){
1233
                // test halfedges around v
1234
                if(j.halfedge() == InvalidHalfEdgeID){
1235
                    cout << "Vertex circulation produced invalid halfedge" << endl;
1236
                    return false;
1237
                }
1238
                VertexID ring_v = j.vertex();
1239
 
1240
                // test one-ring for multiple occurences of vertex
1241
                if(find(link.begin(), link.end(), ring_v) != link.end()){
1242
                    cout << "Vertex appears two times in one-ring of vertex" << endl;
1243
                    return false;
1244
                }
1245
                link.push_back(ring_v);
1246
 
1247
                // test for infinite loop around vertex
1248
                if(static_cast<size_t>(j.no_steps()) > m.no_vertices()){
1249
                    cout << "Vertex loop contains more vertices than manifold" << endl;
1250
                    return false;
1251
                }
1252
            }
1253
 
1254
            // test one_ring size for boundary consistency
1255
            if(link.size() <= 2){
1256
                Walker j = m.walker(*v);
1257
 
1258
                if(j.face() != InvalidFaceID && j.opp().face() != InvalidFaceID)
1259
                {
1260
                    if(link.size()==1)
1261
                        cout << "Vertex contains only a single incident edge" << endl;
1262
                    //
1263
                    //                    cout << "Vertex contains only " << link.size() <<" incident edges" << endl;
1264
                }
1265
                else
1266
                    cout << "Boundary vertex contains only " << link.size() <<" incident edges" << endl;
1267
            }
1268
        }
1269
        // verify components of faces
1270
        for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
1271
            // count edges on face
1272
            Walker j = m.walker(*f);
1273
 
1274
            for(; !j.full_circle(); j = j.circulate_face_cw()){
1275
                // test that all halfedges in faces bind properly to their face
1276
                if(j.face() != *f){
1277
                    cout << "Face is inconsistent, halfedge is not bound to face" << endl;
1278
                    return false;
1279
                }
1280
            }
1281
            // test faces for valid geometrical properties
1282
            if(j.no_steps() < 3){
1283
                cout << "Face contains less than 3 edges" << endl;
1284
                return false;
1285
            }
1286
            // test for infinite loop around face
1287
            if(j.no_steps() > m.no_halfedges() * 0.5f){
1288
                cout << "Face loop contains more halfedges than manifold" << endl;
1289
                return false;
1290
            }
1291
        }
1292
        return true;
1293
    }
1294
 
1295
    void bbox(const Manifold& m, Manifold::Vec& pmin, Manifold::Vec& pmax)
1296
    {
1297
        if(m.no_vertices()==0)
1298
            return;
1299
        VertexIDIterator v = m.vertices_begin();
1300
        pmin = pmax = m.pos(*v);
1301
        ++v;
1302
        for(; v != m.vertices_end(); ++v){
1303
            pmin = v_min(m.pos(*v), pmin);
1304
            pmax = v_max(m.pos(*v), pmax);
1305
        }
1306
    }
1307
 
1308
    void bsphere(const Manifold& m, Manifold::Vec& c, float& r)
1309
    {
1310
        Manifold::Vec p0, p7;
1311
        bbox(m, p0, p7);
1312
        Manifold::Vec rad = (p7 - p0) * 0.5f;
1313
        c = p0 + rad;
1314
        r = rad.length();
1315
    }
1316
 
1317
 
1318
 
1319
    bool precond_collapse_edge(const Manifold& m, HalfEdgeID h)
1320
    {
1321
        Walker hew = m.walker(h);
1322
        HalfEdgeID ho = hew.opp().halfedge();
1323
        VertexID v0 = hew.opp().vertex();
1324
        VertexID v1 = hew.vertex();
1325
 
1326
        // get the one-ring of v0
1327
        vector<VertexID> link0;
1328
        vector<FaceID> faces0;
1329
        for(Walker vj = m.walker(h);
1330
            !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
1331
            link0.push_back(vj.vertex());
1332
            if(vj.halfedge() != h)
1333
                faces0.push_back(vj.face());
1334
        }
1335
        assert(link0.size() > 1);
1336
 
1337
        // get the one-ring of v1
1338
        vector<VertexID> link1;
1339
        vector<FaceID> faces1;
1340
        for(Walker vj = m.walker(ho);
1341
            !vj.full_circle(); vj = vj.circulate_vertex_ccw()){
1342
            link1.push_back(vj.vertex());
1343
            if(vj.halfedge() != ho)
1344
                faces1.push_back(vj.face());
1345
        }
1346
        assert(link1.size() > 1);
1347
 
1348
        // sort the vertices of the two rings
1349
        sort(link0.begin(), link0.end());
1350
        sort(link1.begin(), link1.end());
1351
 
1352
        // get the intersection of the shared vertices from both rings
1353
        vector<VertexID> lisect;
1354
        std::back_insert_iterator<vector<VertexID> > lii(lisect);
1355
        set_intersection(link0.begin(), link0.end(),
1356
                         link1.begin(), link1.end(),
1357
                         lii);
1358
 
1359
        // sort the faces of the two rings
1360
        sort(faces0.begin(), faces0.end());
1361
        sort(faces1.begin(), faces1.end());
1362
 
1363
        // get the intersection of the shared faces from both rings
1364
        vector<FaceID> fisect;
1365
        std::back_insert_iterator<vector<FaceID> > fii(fisect);
1366
        set_intersection(faces0.begin(), faces0.end(),
1367
                         faces1.begin(), faces1.end(),
1368
                         fii);
1369
        if(fisect.size() > 0)
1370
            return false;
1371
 
1372
        int k = 0;
1373
        // if the adjacent face is a triangle (see 2)
1374
        if(hew.next().next().next().halfedge() == h){
1375
            VertexID v = hew.next().vertex();
1376
 
1377
            // valency test (see 5)
1378
            if(valency(m, v) < 3)
1379
                return false;
1380
 
1381
            // remove the vertex shared by the two rings from the intersection
1382
            vector<VertexID>::iterator iter;
1383
            iter = find(lisect.begin(), lisect.end(), v);
1384
            assert(iter != lisect.end());
1385
            lisect.erase(iter);
1386
            ++k;
1387
        }
1388
        // if the adjacent face is a triangle (see 2)
1389
        if(hew.opp().next().next().next().halfedge() == hew.opp().halfedge()){
1390
            VertexID v = hew.opp().next().vertex();
1391
 
1392
            // valency test (see 5)
1393
            if(valency(m, v) < 3)
1394
                return false;
1395
 
1396
            // remove the vertex shared by the two rings from the intersection
1397
            vector<VertexID>::iterator iter;
1398
            iter = find(lisect.begin(), lisect.end(), v);
1399
            assert(iter != lisect.end());
1400
            lisect.erase(iter);
1401
            ++k;
1402
        }
1403
        // double edge test (see 3)
1404
        if(lisect.size() != 0)
1405
            return false;
1406
 
1407
        // tetrahedon test (see 4)
1408
        if(k == 2 && (link0.size() + link1.size() == 6))
1409
            return false;
1410
 
1411
        // test that we do not merge holes (see 6)
1412
        if(boundary(m, v0) && boundary(m, v1) && hew.face() != InvalidFaceID && hew.opp().face() != InvalidFaceID)
1413
            return false;
1414
 
1415
        return true;
1416
    }
1417
 
1418
    bool precond_flip_edge(const Manifold& m, HalfEdgeID h)
1419
    {
1420
        Walker j = m.walker(h);
1421
 
1422
        FaceID hf = j.face();
1423
        FaceID hof = j.opp().face();
1424
 
1425
        // boundary case
1426
        if(hf == InvalidFaceID || hof == InvalidFaceID)
1427
           return false;
1428
 
1429
 
1430
        // We can only flip an edge if both incident polygons are triangles.
1431
        if(no_edges(m, hf) != 3 || no_edges(m, hof) !=3)
1432
            return false;
1433
 
1434
 
1435
        // non boundary vertices with a valency of less than 4(less than 3 after operation) degenerates mesh.
1436
        VertexID hv = j.vertex();
1437
        VertexID hov = j.opp().vertex();
1438
        if((valency(m, hv) < 4 && !boundary(m, hv)) || (valency(m, hov) < 4 && !boundary(m, hov))){
1439
            return false;
1440
        }
1441
 
1442
        // Disallow flip if vertices being connected already are.
1443
        VertexID hnv = j.next().vertex();
1444
        VertexID honv = j.opp().next().vertex();
1445
        if(connected(m, hnv, honv)){
1446
           return false;
1447
        }
1448
 
1449
        return true;
1450
    }
1451
 
1452
    bool boundary(const Manifold& m, VertexID v)
1453
    {
1454
        Walker j  = m.walker(v);
1455
        return boundary(m, j.halfedge());
1456
    }
1457
 
1458
    int valency(const Manifold& m, VertexID v)
1459
    {
1460
        return circulate_vertex_ccw(m,v, (std::function<void(Walker&)>)[](Walker){});
1461
    }
1462
 
1463
    Manifold::Vec normal(const Manifold& m, VertexID v)
1464
    {
1465
        Manifold::Vec p0 = m.pos(v);
1466
        vector<Manifold::Vec> one_ring;
1467
 
1468
        // run through outgoing edges, and store them normalized
1469
        circulate_vertex_ccw(m, v, (std::function<void(VertexID)>)[&](VertexID vn) {
1470
            Manifold::Vec edge = m.pos(vn) - p0;
1471
            double l = length(edge);
1472
            if(l > 0.0)
1473
                one_ring.push_back(edge/l);
1474
        });
1475
        int N = one_ring.size();
1476
        if(N<2)
1477
            return Manifold::Vec(0);
1478
 
1479
        size_t N_count = N;
1480
        size_t N_start = 0;
1481
        if(boundary(m, v))
1482
            N_start = 1;
1483
 
1484
        // sum up the normals of each face surrounding the vertex
1485
        Manifold::Vec n(0);
1486
        for(size_t i = N_start; i < N_count; ++i){
1487
            Manifold::Vec e0 = one_ring[i];
1488
            Manifold::Vec e1 = one_ring[(i+1) % N];
1489
 
1490
            Manifold::Vec n_part = normalize(cross(e0, e1));
1491
            n += n_part * acos(max(-1.0, min(1.0, dot(e0, e1))));
1492
        }
1493
 
1494
        // normalize and return the normal
1495
        float sqr_l = sqr_length(n);
1496
        if(sqr_l > 0.0f) return n / sqrt(sqr_l);
1497
 
1498
        return n;
1499
    }
1500
 
1501
 
1502
    bool connected(const Manifold& m, VertexID v0, VertexID v1)
1503
    {
1504
        bool c=false;
1505
        circulate_vertex_ccw(m, v0, (std::function<void(VertexID)>)[&](VertexID v){ c |= (v==v1);});
1506
        return c;
1507
    }
1508
 
1509
 
1510
    int no_edges(const Manifold& m, FaceID f)
1511
    {
1512
        return circulate_face_ccw(m, f, (std::function<void(Walker&)>)[](Walker w){});
1513
    }
1514
 
1515
    Manifold::Vec normal(const Manifold& m, FaceID f)
1516
    {
1517
        vector<Manifold::Vec> v;
1518
 
1519
        int k= circulate_face_ccw(m, f, (std::function<void(VertexID)>)[&](VertexID vid) {
1520
            v.push_back(m.pos(vid));
1521
        });
1522
 
1523
        Manifold::Vec norm(0);
1524
        for(int i=0;i<k;++i)
1525
        {
1526
            norm[0] += (v[i][1]-v[(i+1)%k][1])*(v[i][2]+v[(i+1)%k][2]);
1527
            norm[1] += (v[i][2]-v[(i+1)%k][2])*(v[i][0]+v[(i+1)%k][0]);
1528
            norm[2] += (v[i][0]-v[(i+1)%k][0])*(v[i][1]+v[(i+1)%k][1]);
1529
        }
1530
        float l = sqr_length(norm);
1531
        if(l>0.0f)
1532
            norm /= sqrt(l);
1533
        return norm;
1534
    }
1535
 
1536
 
1537
    double area(const Manifold& m, FaceID fid)
1538
    {
1539
        // Get all projected vertices
1540
        vector<Manifold::Vec> vertices;
1541
        int N = circulate_face_ccw(m, fid, (std::function<void(VertexID)>)[&](VertexID vid) {
1542
            vertices.push_back(m.pos(vid));
1543
        });
1544
 
1545
 
1546
        double area = 0;
1547
        Manifold::Vec norm = normal(m,fid);
1548
        for(int i = 1; i < N-1; ++i)
1549
            area += 0.5 * dot(norm,cross(vertices[i]-vertices[0], vertices[(i+1 )]-vertices[0]));
1550
        return area;
1551
    }
1552
 
1553
    Manifold::Vec centre(const Manifold& m, FaceID f)
1554
    {
1555
        Manifold::Vec c(0);
1556
        int n = circulate_face_ccw(m, f, (std::function<void(VertexID)>)[&](VertexID v) {c+=m.pos(v);});
1557
        return c / n;
1558
    }
1559
 
1560
    double perimeter(const Manifold& m, FaceID f)
1561
    {
1562
        double l=0.0;
1563
        circulate_face_ccw(m, f, (std::function<void(HalfEdgeID)>)[&](HalfEdgeID h) { l+= length(m, h);});
1564
        return l;
1565
    }
1566
 
1567
    bool boundary(const Manifold& m, HalfEdgeID h)
1568
    {
1569
        Walker w = m.walker(h);
1570
        return w.face() == InvalidFaceID || w.opp().face() == InvalidFaceID;
1571
    }
1572
 
1573
    double length(const Manifold& m, HalfEdgeID h)
1574
    {
1575
        Walker w = m.walker(h);
1576
        return (m.pos(w.vertex()) - m.pos(w.opp().vertex())).length();
1577
    }
1578
}