Subversion Repositories gelsvn

Rev

Rev 601 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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