Subversion Repositories gelsvn

Rev

Rev 561 | Rev 601 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

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