Subversion Repositories gelsvn

Rev

Details | Last modification | View Log | RSS feed

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