Subversion Repositories gelsvn

Rev

Rev 631 | Rev 636 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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