Subversion Repositories gelsvn

Rev

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

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