Subversion Repositories gelsvn

Rev

Rev 589 | Details | Compare with Previous | Last modification | View Log | RSS feed

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