WebSVN – gelsvn – Diff – /branches/cpp11-devel/src/HMesh/triangulate.cpp


/* ----------------------------------------------------------------------- *
 * This file is part of GEL, http://www.imm.dtu.dk/GEL
 * Copyright (C) the authors and DTU Informatics
 * For license and list of authors, see ../../doc/intro.pdf
 * ----------------------------------------------------------------------- */
 
#include "triangulate.h"
 
#include <queue>
#include <vector>
#include <iterator>
#include <cassert>
 
#include "../CGLA/Vec3d.h"
 
#include "Manifold.h"
#include "AttributeVector.h"
 
namespace HMesh
{
    using namespace std;
    using namespace CGLA;
 
    void get_candidates(const Manifold& m, VertexID v, vector<HalfEdgeID>& candidates)
    {
        vector<VertexID> vertices;
        vector<HalfEdgeID> hedges;
 
        Walker wv = m.walker(v);
        for(;!wv.full_circle(); wv = wv.circulate_vertex_cw()){
            vertices.push_back(wv.vertex());
            hedges.push_back(wv.halfedge());
        }
        int N = wv.no_steps();
        vector<VertexID> vertices_check(vertices);
        assert(N == vertices.size());
 
        for(int i = N - 1; i >= 0; --i){
            for(Walker w = m.walker(hedges[i]).next(); w.vertex() != vertices[(i+N-1)%N]; w = w.next()){
                if(find(vertices_check.begin(), vertices_check.end(), w.vertex()) == vertices_check.end())
                    candidates.push_back(w.halfedge());       
            }
        }
    }
 
    float curv(const Vec3d& p, vector<Vec3d>& vec)
    {
        size_t N = vec.size();
        std::vector<Vec3d> normals;
        for(size_t i = 0; i < N; ++i){
            Vec3d norm = normalize(cross(vec[i]-p, vec[(i+1)%N]-p));
            normals.push_back(norm);
        }
        float alpha = 0;
        for(size_t i = 0; i < N; ++i)
            alpha += (vec[i]-p).length()*acos(dot(normals[i],normals[(i+1)%N]));
 
        return alpha;
    }
 
    float get_badness(const Manifold& m, VertexID v, VertexID n)
    {
        vector<HalfEdgeID> hedges;
 
        Walker wv = m.walker(v);
        for(; !wv.full_circle(); wv = wv.circulate_vertex_cw())
            hedges.push_back(wv.halfedge());
 
        vector<Vec3d> one_ring;
        vector<Vec3d> one_ring_n;
        int N = wv.no_steps();
        for(int i = N - 1; i >= 0; --i){
            for(Walker w = m.walker(hedges[i]).next(); w.vertex() != v; w = w.next()){
                one_ring.push_back(m.pos(w.vertex()));
                if(w.vertex() != n)
                    one_ring_n.push_back(m.pos(w.vertex()));
            }
        }
        return curv(m.pos(v), one_ring) - curv(m.pos(v), one_ring_n);
    }
 
 
    const CGLA::Vec3d get_normal(const Manifold& m, VertexID v)
    {
 
        vector<HalfEdgeID> hedges;
 
        Walker wv = m.walker(v);
        for(; !wv.full_circle(); wv = wv.circulate_vertex_cw())
            hedges.push_back(wv.halfedge());
 
        vector<Vec3d> one_ring;
        size_t N = wv.no_steps();
        for(int i = N - 1; i >= 0; --i){
            for(Walker w = m.walker(hedges[i]).next(); w.vertex() != v; w = w.next())
                one_ring.push_back(m.pos(w.vertex()));   
        }
 
        Vec3d norm(0);
        N = one_ring.size();
        Vec3d p0 = m.pos(v);
        for(size_t i = 0; i < N; ++i){
            Vec3d p1 = one_ring[i];
            Vec3d p2 = one_ring[(i+1) % N];
            Vec3d e0 = normalize(p1 - p0);
            Vec3d e1 = normalize(p2 - p0);
            norm += cross(e0, e1) * acos(dot(e0, e1));
        }
        return normalize(norm);
    }
 
    void triangulate_by_vertex_face_split(Manifold& m)
    {
        vector<FaceID> fv;
        fv.reserve(m.no_faces());
        copy(m.faces_begin(), m.faces_end(), back_inserter(fv));
 
        for(size_t i = 0; i < fv.size(); ++i)
            m.split_face_by_vertex(fv[i]);
    }
 
    void triangulate_by_edge_face_split(Manifold& m)
    {
        vector<FaceID> fv;
        fv.reserve(m.no_faces());
        copy(m.faces_begin(), m.faces_end(), back_inserter(fv));
 
        for(size_t i = 0; i < fv.size(); ++i)
            triangulate_face_by_edge_split(m, fv[i]);
    }
 
    struct PotentialEdge
    {
        int time_tag;
        float badness;
        FaceID f;
        VertexID v0, v1;
    };
 
    bool operator>(const PotentialEdge& e0, const PotentialEdge& e1)
    { return e0.badness>e1.badness; }
 
    typedef std::priority_queue<PotentialEdge,
        std::vector<PotentialEdge>,
        std::greater<PotentialEdge> > 
        PotentialEdgeQueue;
 
    void insert_potential_edges(const Manifold& m, VertexAttributeVector<int>& vtouched, VertexID v, PotentialEdgeQueue& pot_edges)
    {
        vector<Vec3d> one_ring;
        vector<HalfEdgeID> candidates;
        get_candidates(m, v, candidates);
 
        Vec3d p0 = m.pos(v);
        Vec3d norm = normal(m, v);
        int n = candidates.size();
        for(int i = 0; i < n; ++i){
            Walker w = m.walker(candidates[i]);
            VertexID v_n = w.vertex();
            Vec3d edir = normalize(m.pos(v_n) - p0);
            Vec3d norm_n = normal(m, v_n);
            float bad = sqr(dot(edir, norm));
            float bad_n = sqr(dot(edir, norm_n));
 
            PotentialEdge pot;
 
            /* So if the edge between two vertices is not orthogonal to 
            their normals, the badness is increased. Badness is also
            increased if the normals are very different. */
 
            pot.badness = bad+bad_n - dot(norm_n,norm);
            pot.time_tag = vtouched[v];
            pot.v0 = v;
            pot.v1 = w.vertex();
            pot.f = w.face();
 
            pot_edges.push(pot);
        }
    }
 
    void curvature_triangulate(Manifold& m)
    {
        PotentialEdgeQueue pot_edges;
        VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);
 
        // Create candidates
        for(VertexIDIterator v = m.vertices_begin(); v!= m.vertices_end(); ++v)
            insert_potential_edges(m, vtouched, *v, pot_edges);
 
        while(!pot_edges.empty()){
            const PotentialEdge& pot_edge = pot_edges.top();
            // Record all the vertices of the face. We need to 
            // recompute the candidates.
            std::vector<VertexID> reeval_vec;
 
            for(Walker w = m.walker(pot_edge.f); !w.full_circle(); w = w.circulate_face_cw())
                reeval_vec.push_back(w.vertex());
 
            // Make sure that the vertex has not been touched since 
            // we created the potential edge. If the vertex has been
            // touched the candidate edge has either (a) been processed,
            // (b) received a lower priority or (c) become invalid.
            if(pot_edge.time_tag == vtouched[pot_edge.v0]){
                vector<Vec3d> one_ring;
                vector<HalfEdgeID> candidates;
 
                m.split_face_by_edge(pot_edge.f, pot_edge.v0, pot_edge.v1);
 
                // Recompute priorities.
                for(size_t i = 0; i < reeval_vec.size(); ++i){
                    VertexID& v = reeval_vec[i];
                    ++vtouched[v];
                    insert_potential_edges(m, vtouched, v, pot_edges);
                }
 
            }
            pot_edges.pop();
        }
 
    }
 
    void shortest_edge_triangulate(Manifold& m)
    {
        int work;
        do{
            // For every face.
            work = 0;
            for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){
                // Create a vector of vertices.
                vector<VertexID> verts;
                for(Walker w = m.walker(*f); !w.full_circle(); w = w.circulate_face_ccw())
				{
					FaceID fa = w.face();
					FaceID fb = *f;
					assert(fa==fb);
                    verts.push_back(w.vertex());
				}
                // If there are just three we are done.
                if(verts.size() == 3) continue;
 
                // Find vertex pairs that may be connected.
                vector<pair<int,int> > vpairs;
                const int N = verts.size();
                for(int i = 0; i < N - 2; ++i){
                    for(int j = i + 2; j < N; ++j){
                        if(verts[i] != verts[j] && !connected(m, verts[i], verts[j]))
                            vpairs.push_back(pair<int,int>(i, j));
                    }
                }
                if(vpairs.empty()){
                    cout << "Warning: could not triangulate a face." 
                        << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
                    continue;
                }
 
                /* For all vertex pairs, find the edge lengths. Combine the
                vertices forming the shortest edge. */
 
                float min_len=FLT_MAX;
                int min_k = -1;
                for(size_t k = 0; k < vpairs.size(); ++k){
                    int i = vpairs[k].first;
                    int j = vpairs[k].second;
                    float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));
 
                    if(len<min_len){
                        min_len = len;
                        min_k = k;
                    }
                }
                assert(min_k != -1);
 
                {
                    // Split faces along edge whose midpoint is closest to isovalue
                    int i = vpairs[min_k].first;
                    int j = vpairs[min_k].second;
                    m.split_face_by_edge(*f, verts[i], verts[j]);
                }
                ++work;
 
            }
        }
        while(work);
    }
 
    void triangulate_face_by_edge_split(Manifold& m, FaceID f)
    {
        if(no_edges(m, f)<=3)
            return;
        
        Walker w = m.walker(f);
 
        // as long as f is not a triangle
        while(w.next().next().next().halfedge() != w.halfedge()){
            // assert that face has at least 3 edges
            // f is split into triangle from first three vertices, and becomes remaining polygon in next iteration
            assert(w.next().next().halfedge() != w.halfedge());
            VertexID v0 = w.vertex();
            VertexID v1 = w.next().next().vertex();
            FaceID f_old = f;
            if(v0 != v1 && !connected(m, v0, v1))
                f = m.split_face_by_edge(f, v0, v1);
            if(f == f_old)
                return;
        }
    }
}
 

Rev 607	Rev 630
1	`/* ----------------------------------------------------------------------- *`	1	`/* ----------------------------------------------------------------------- *`
2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`	2	`* This file is part of GEL, http://www.imm.dtu.dk/GEL`
3	`* Copyright (C) the authors and DTU Informatics`	3	`* Copyright (C) the authors and DTU Informatics`
4	`* For license and list of authors, see ../../doc/intro.pdf`	4	`* For license and list of authors, see ../../doc/intro.pdf`
5	`* ----------------------------------------------------------------------- */`	5	`* ----------------------------------------------------------------------- */`
6		6
7	`#include "triangulate.h"`	7	`#include "triangulate.h"`
8		8
9	`#include <queue>`	9	`#include <queue>`
10	`#include <vector>`	10	`#include <vector>`
11	`#include <iterator>`	11	`#include <iterator>`
12	`#include <cassert>`	12	`#include <cassert>`
13		13
14	`#include "../CGLA/Vec3d.h"`	14	`#include "../CGLA/Vec3d.h"`
15		15
16	`#include "Manifold.h"`	16	`#include "Manifold.h"`
17	`#include "AttributeVector.h"`	17	`#include "AttributeVector.h"`
18		18
19	`namespace HMesh`	19	`namespace HMesh`
20	`{`	20	`{`
21	`using namespace std;`	21	`using namespace std;`
22	`using namespace CGLA;`	22	`using namespace CGLA;`
23		23
24	`void get_candidates(const Manifold& m, VertexID v, vector<HalfEdgeID>& candidates)`	24	`void get_candidates(const Manifold& m, VertexID v, vector<HalfEdgeID>& candidates)`
25	`{`	25	`{`
26	`vector<VertexID> vertices;`	26	`vector<VertexID> vertices;`
27	`vector<HalfEdgeID> hedges;`	27	`vector<HalfEdgeID> hedges;`
28		28
29	`Walker wv = m.walker(v);`	29	`Walker wv = m.walker(v);`
30	`for(;!wv.full_circle(); wv = wv.circulate_vertex_cw()){`	30	`for(;!wv.full_circle(); wv = wv.circulate_vertex_cw()){`
31	`vertices.push_back(wv.vertex());`	31	`vertices.push_back(wv.vertex());`
32	`hedges.push_back(wv.halfedge());`	32	`hedges.push_back(wv.halfedge());`
33	`}`	33	`}`
34	`int N = wv.no_steps();`	34	`int N = wv.no_steps();`
35	`vector<VertexID> vertices_check(vertices);`	35	`vector<VertexID> vertices_check(vertices);`
36	`assert(N == vertices.size());`	36	`assert(N == vertices.size());`
37		37
38	`for(int i = N - 1; i >= 0; --i){`	38	`for(int i = N - 1; i >= 0; --i){`
39	`for(Walker w = m.walker(hedges[i]).next(); w.vertex() != vertices[(i+N-1)%N]; w = w.next()){`	39	`for(Walker w = m.walker(hedges[i]).next(); w.vertex() != vertices[(i+N-1)%N]; w = w.next()){`
40	`if(find(vertices_check.begin(), vertices_check.end(), w.vertex()) == vertices_check.end())`	40	`if(find(vertices_check.begin(), vertices_check.end(), w.vertex()) == vertices_check.end())`
41	`candidates.push_back(w.halfedge());`	41	`candidates.push_back(w.halfedge());`
42	`}`	42	`}`
43	`}`	43	`}`
44	`}`	44	`}`
45		45
46	`float curv(const Vec3d& p, vector<Vec3d>& vec)`	46	`float curv(const Vec3d& p, vector<Vec3d>& vec)`
47	`{`	47	`{`
48	`size_t N = vec.size();`	48	`size_t N = vec.size();`
49	`std::vector<Vec3d> normals;`	49	`std::vector<Vec3d> normals;`
50	`for(size_t i = 0; i < N; ++i){`	50	`for(size_t i = 0; i < N; ++i){`
51	`Vec3d norm = normalize(cross(vec[i]-p, vec[(i+1)%N]-p));`	51	`Vec3d norm = normalize(cross(vec[i]-p, vec[(i+1)%N]-p));`
52	`normals.push_back(norm);`	52	`normals.push_back(norm);`
53	`}`	53	`}`
54	`float alpha = 0;`	54	`float alpha = 0;`
55	`for(size_t i = 0; i < N; ++i)`	55	`for(size_t i = 0; i < N; ++i)`
56	`alpha += (vec[i]-p).length()*acos(dot(normals[i],normals[(i+1)%N]));`	56	`alpha += (vec[i]-p).length()*acos(dot(normals[i],normals[(i+1)%N]));`
57		57
58	`return alpha;`	58	`return alpha;`
59	`}`	59	`}`
60		60
61	`float get_badness(const Manifold& m, VertexID v, VertexID n)`	61	`float get_badness(const Manifold& m, VertexID v, VertexID n)`
62	`{`	62	`{`
63	`vector<HalfEdgeID> hedges;`	63	`vector<HalfEdgeID> hedges;`
64		64
65	`Walker wv = m.walker(v);`	65	`Walker wv = m.walker(v);`
66	`for(; !wv.full_circle(); wv = wv.circulate_vertex_cw())`	66	`for(; !wv.full_circle(); wv = wv.circulate_vertex_cw())`
67	`hedges.push_back(wv.halfedge());`	67	`hedges.push_back(wv.halfedge());`
68		68
69	`vector<Vec3d> one_ring;`	69	`vector<Vec3d> one_ring;`
70	`vector<Vec3d> one_ring_n;`	70	`vector<Vec3d> one_ring_n;`
71	`int N = wv.no_steps();`	71	`int N = wv.no_steps();`
72	`for(int i = N - 1; i >= 0; --i){`	72	`for(int i = N - 1; i >= 0; --i){`
73	`for(Walker w = m.walker(hedges[i]).next(); w.vertex() != v; w = w.next()){`	73	`for(Walker w = m.walker(hedges[i]).next(); w.vertex() != v; w = w.next()){`
74	`one_ring.push_back(m.pos(w.vertex()));`	74	`one_ring.push_back(m.pos(w.vertex()));`
75	`if(w.vertex() != n)`	75	`if(w.vertex() != n)`
76	`one_ring_n.push_back(m.pos(w.vertex()));`	76	`one_ring_n.push_back(m.pos(w.vertex()));`
77	`}`	77	`}`
78	`}`	78	`}`
79	`return curv(m.pos(v), one_ring) - curv(m.pos(v), one_ring_n);`	79	`return curv(m.pos(v), one_ring) - curv(m.pos(v), one_ring_n);`
80	`}`	80	`}`
81		81
82		82
83	`const CGLA::Vec3d get_normal(const Manifold& m, VertexID v)`	83	`const CGLA::Vec3d get_normal(const Manifold& m, VertexID v)`
84	`{`	84	`{`
85		85
86	`vector<HalfEdgeID> hedges;`	86	`vector<HalfEdgeID> hedges;`
87		87
88	`Walker wv = m.walker(v);`	88	`Walker wv = m.walker(v);`
89	`for(; !wv.full_circle(); wv = wv.circulate_vertex_cw())`	89	`for(; !wv.full_circle(); wv = wv.circulate_vertex_cw())`
90	`hedges.push_back(wv.halfedge());`	90	`hedges.push_back(wv.halfedge());`
91		91
92	`vector<Vec3d> one_ring;`	92	`vector<Vec3d> one_ring;`
93	`size_t N = wv.no_steps();`	93	`size_t N = wv.no_steps();`
94	`for(int i = N - 1; i >= 0; --i){`	94	`for(int i = N - 1; i >= 0; --i){`
95	`for(Walker w = m.walker(hedges[i]).next(); w.vertex() != v; w = w.next())`	95	`for(Walker w = m.walker(hedges[i]).next(); w.vertex() != v; w = w.next())`
96	`one_ring.push_back(m.pos(w.vertex()));`	96	`one_ring.push_back(m.pos(w.vertex()));`
97	`}`	97	`}`
98		98
99	`Vec3d norm(0);`	99	`Vec3d norm(0);`
100	`N = one_ring.size();`	100	`N = one_ring.size();`
101	`Vec3d p0 = m.pos(v);`	101	`Vec3d p0 = m.pos(v);`
102	`for(size_t i = 0; i < N; ++i){`	102	`for(size_t i = 0; i < N; ++i){`
103	`Vec3d p1 = one_ring[i];`	103	`Vec3d p1 = one_ring[i];`
104	`Vec3d p2 = one_ring[(i+1) % N];`	104	`Vec3d p2 = one_ring[(i+1) % N];`
105	`Vec3d e0 = normalize(p1 - p0);`	105	`Vec3d e0 = normalize(p1 - p0);`
106	`Vec3d e1 = normalize(p2 - p0);`	106	`Vec3d e1 = normalize(p2 - p0);`
107	`norm += cross(e0, e1) * acos(dot(e0, e1));`	107	`norm += cross(e0, e1) * acos(dot(e0, e1));`
108	`}`	108	`}`
109	`return normalize(norm);`	109	`return normalize(norm);`
110	`}`	110	`}`
111		111
112	`void triangulate_by_vertex_face_split(Manifold& m)`	112	`void triangulate_by_vertex_face_split(Manifold& m)`
113	`{`	113	`{`
114	`vector<FaceID> fv;`	114	`vector<FaceID> fv;`
115	`fv.reserve(m.no_faces());`	115	`fv.reserve(m.no_faces());`
116	`copy(m.faces_begin(), m.faces_end(), back_inserter(fv));`	116	`copy(m.faces_begin(), m.faces_end(), back_inserter(fv));`
117		117
118	`for(size_t i = 0; i < fv.size(); ++i)`	118	`for(size_t i = 0; i < fv.size(); ++i)`
119	`m.split_face_by_vertex(fv[i]);`	119	`m.split_face_by_vertex(fv[i]);`
120	`}`	120	`}`
121		121
122	`void triangulate_by_edge_face_split(Manifold& m)`	122	`void triangulate_by_edge_face_split(Manifold& m)`
123	`{`	123	`{`
124	`vector<FaceID> fv;`	124	`vector<FaceID> fv;`
125	`fv.reserve(m.no_faces());`	125	`fv.reserve(m.no_faces());`
126	`copy(m.faces_begin(), m.faces_end(), back_inserter(fv));`	126	`copy(m.faces_begin(), m.faces_end(), back_inserter(fv));`
127		127
128	`for(size_t i = 0; i < fv.size(); ++i)`	128	`for(size_t i = 0; i < fv.size(); ++i)`
129	`triangulate_face_by_edge_split(m, fv[i]);`	129	`triangulate_face_by_edge_split(m, fv[i]);`
130	`}`	130	`}`
131		131
132	`struct PotentialEdge`	132	`struct PotentialEdge`
133	`{`	133	`{`
134	`int time_tag;`	134	`int time_tag;`
135	`float badness;`	135	`float badness;`
136	`FaceID f;`	136	`FaceID f;`
137	`VertexID v0, v1;`	137	`VertexID v0, v1;`
138	`};`	138	`};`
139		139
140	`bool operator>(const PotentialEdge& e0, const PotentialEdge& e1)`	140	`bool operator>(const PotentialEdge& e0, const PotentialEdge& e1)`
141	`{ return e0.badness>e1.badness; }`	141	`{ return e0.badness>e1.badness; }`
142		142
143	`typedef std::priority_queue<PotentialEdge,`	143	`typedef std::priority_queue<PotentialEdge,`
144	`std::vector<PotentialEdge>,`	144	`std::vector<PotentialEdge>,`
145	`std::greater<PotentialEdge> >`	145	`std::greater<PotentialEdge> >`
146	`PotentialEdgeQueue;`	146	`PotentialEdgeQueue;`
147		147
148	`void insert_potential_edges(const Manifold& m, VertexAttributeVector<int>& vtouched, VertexID v, PotentialEdgeQueue& pot_edges)`	148	`void insert_potential_edges(const Manifold& m, VertexAttributeVector<int>& vtouched, VertexID v, PotentialEdgeQueue& pot_edges)`
149	`{`	149	`{`
150	`vector<Vec3d> one_ring;`	150	`vector<Vec3d> one_ring;`
151	`vector<HalfEdgeID> candidates;`	151	`vector<HalfEdgeID> candidates;`
152	`get_candidates(m, v, candidates);`	152	`get_candidates(m, v, candidates);`
153		153
154	`Vec3d p0 = m.pos(v);`	154	`Vec3d p0 = m.pos(v);`
155	`Vec3d norm = normal(m, v);`	155	`Vec3d norm = normal(m, v);`
156	`int n = candidates.size();`	156	`int n = candidates.size();`
157	`for(int i = 0; i < n; ++i){`	157	`for(int i = 0; i < n; ++i){`
158	`Walker w = m.walker(candidates[i]);`	158	`Walker w = m.walker(candidates[i]);`
159	`VertexID v_n = w.vertex();`	159	`VertexID v_n = w.vertex();`
160	`Vec3d edir = normalize(m.pos(v_n) - p0);`	160	`Vec3d edir = normalize(m.pos(v_n) - p0);`
161	`Vec3d norm_n = normal(m, v_n);`	161	`Vec3d norm_n = normal(m, v_n);`
162	`float bad = sqr(dot(edir, norm));`	162	`float bad = sqr(dot(edir, norm));`
163	`float bad_n = sqr(dot(edir, norm_n));`	163	`float bad_n = sqr(dot(edir, norm_n));`
164		164
165	`PotentialEdge pot;`	165	`PotentialEdge pot;`
166		166
167	`/* So if the edge between two vertices is not orthogonal to`	167	`/* So if the edge between two vertices is not orthogonal to`
168	`their normals, the badness is increased. Badness is also`	168	`their normals, the badness is increased. Badness is also`
169	`increased if the normals are very different. */`	169	`increased if the normals are very different. */`
170		170
171	`pot.badness = bad+bad_n - dot(norm_n,norm);`	171	`pot.badness = bad+bad_n - dot(norm_n,norm);`
172	`pot.time_tag = vtouched[v];`	172	`pot.time_tag = vtouched[v];`
173	`pot.v0 = v;`	173	`pot.v0 = v;`
174	`pot.v1 = w.vertex();`	174	`pot.v1 = w.vertex();`
175	`pot.f = w.face();`	175	`pot.f = w.face();`
176		176
177	`pot_edges.push(pot);`	177	`pot_edges.push(pot);`
178	`}`	178	`}`
179	`}`	179	`}`
180		180
181	`void curvature_triangulate(Manifold& m)`	181	`void curvature_triangulate(Manifold& m)`
182	`{`	182	`{`
183	`PotentialEdgeQueue pot_edges;`	183	`PotentialEdgeQueue pot_edges;`
184	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`	184	`VertexAttributeVector<int> vtouched(m.allocated_vertices(), 0);`
185		185
186	`// Create candidates`	186	`// Create candidates`
187	`for(VertexIDIterator v = m.vertices_begin(); v!= m.vertices_end(); ++v)`	187	`for(VertexIDIterator v = m.vertices_begin(); v!= m.vertices_end(); ++v)`
188	`insert_potential_edges(m, vtouched, *v, pot_edges);`	188	`insert_potential_edges(m, vtouched, *v, pot_edges);`
189		189
190	`while(!pot_edges.empty()){`	190	`while(!pot_edges.empty()){`
191	`const PotentialEdge& pot_edge = pot_edges.top();`	191	`const PotentialEdge& pot_edge = pot_edges.top();`
192	`// Record all the vertices of the face. We need to`	192	`// Record all the vertices of the face. We need to`
193	`// recompute the candidates.`	193	`// recompute the candidates.`
194	`std::vector<VertexID> reeval_vec;`	194	`std::vector<VertexID> reeval_vec;`
195		195
196	`for(Walker w = m.walker(pot_edge.f); !w.full_circle(); w = w.circulate_face_cw())`	196	`for(Walker w = m.walker(pot_edge.f); !w.full_circle(); w = w.circulate_face_cw())`
197	`reeval_vec.push_back(w.vertex());`	197	`reeval_vec.push_back(w.vertex());`
198		198
199	`// Make sure that the vertex has not been touched since`	199	`// Make sure that the vertex has not been touched since`
200	`// we created the potential edge. If the vertex has been`	200	`// we created the potential edge. If the vertex has been`
201	`// touched the candidate edge has either (a) been processed,`	201	`// touched the candidate edge has either (a) been processed,`
202	`// (b) received a lower priority or (c) become invalid.`	202	`// (b) received a lower priority or (c) become invalid.`
203	`if(pot_edge.time_tag == vtouched[pot_edge.v0]){`	203	`if(pot_edge.time_tag == vtouched[pot_edge.v0]){`
204	`vector<Vec3d> one_ring;`	204	`vector<Vec3d> one_ring;`
205	`vector<HalfEdgeID> candidates;`	205	`vector<HalfEdgeID> candidates;`
206		206
207	`m.split_face_by_edge(pot_edge.f, pot_edge.v0, pot_edge.v1);`	207	`m.split_face_by_edge(pot_edge.f, pot_edge.v0, pot_edge.v1);`
208		208
209	`// Recompute priorities.`	209	`// Recompute priorities.`
210	`for(size_t i = 0; i < reeval_vec.size(); ++i){`	210	`for(size_t i = 0; i < reeval_vec.size(); ++i){`
211	`VertexID& v = reeval_vec[i];`	211	`VertexID& v = reeval_vec[i];`
212	`++vtouched[v];`	212	`++vtouched[v];`
213	`insert_potential_edges(m, vtouched, v, pot_edges);`	213	`insert_potential_edges(m, vtouched, v, pot_edges);`
214	`}`	214	`}`
215		215
216	`}`	216	`}`
217	`pot_edges.pop();`	217	`pot_edges.pop();`
218	`}`	218	`}`
219		219
220	`}`	220	`}`
221		221
222	`void shortest_edge_triangulate(Manifold& m)`	222	`void shortest_edge_triangulate(Manifold& m)`
223	`{`	223	`{`
224	`int work;`	224	`int work;`
225	`do{`	225	`do{`
226	`// For every face.`	226	`// For every face.`
227	`work = 0;`	227	`work = 0;`
228	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){`	228	`for(FaceIDIterator f = m.faces_begin(); f != m.faces_end(); ++f){`
229	`// Create a vector of vertices.`	229	`// Create a vector of vertices.`
230	`vector<VertexID> verts;`	230	`vector<VertexID> verts;`
231	`for(Walker w = m.walker(*f); !w.full_circle(); w = w.circulate_face_ccw())`	231	`for(Walker w = m.walker(*f); !w.full_circle(); w = w.circulate_face_ccw())`
232	`{`	232	`{`
233	`FaceID fa = w.face();`	233	`FaceID fa = w.face();`
234	`FaceID fb = *f;`	234	`FaceID fb = *f;`
235	`assert(fa==fb);`	235	`assert(fa==fb);`
236	`verts.push_back(w.vertex());`	236	`verts.push_back(w.vertex());`
237	`}`	237	`}`
238	`// If there are just three we are done.`	238	`// If there are just three we are done.`
239	`if(verts.size() == 3) continue;`	239	`if(verts.size() == 3) continue;`
240		240
241	`// Find vertex pairs that may be connected.`	241	`// Find vertex pairs that may be connected.`
242	`vector<pair<int,int> > vpairs;`	242	`vector<pair<int,int> > vpairs;`
243	`const int N = verts.size();`	243	`const int N = verts.size();`
244	`for(int i = 0; i < N - 2; ++i){`	244	`for(int i = 0; i < N - 2; ++i){`
245	`for(int j = i + 2; j < N; ++j){`	245	`for(int j = i + 2; j < N; ++j){`
246	`if(verts[i] != verts[j] && !connected(m, verts[i], verts[j]))`	246	`if(verts[i] != verts[j] && !connected(m, verts[i], verts[j]))`
247	`vpairs.push_back(pair<int,int>(i, j));`	247	`vpairs.push_back(pair<int,int>(i, j));`
248	`}`	248	`}`
249	`}`	249	`}`
250	`if(vpairs.empty()){`	250	`if(vpairs.empty()){`
251	`cout << "Warning: could not triangulate a face."`	251	`cout << "Warning: could not triangulate a face."`
252	`<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;`	252	`<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;`
253	`continue;`	253	`continue;`
254	`}`	254	`}`
255		255
256	`/* For all vertex pairs, find the edge lengths. Combine the`	256	`/* For all vertex pairs, find the edge lengths. Combine the`
257	`vertices forming the shortest edge. */`	257	`vertices forming the shortest edge. */`
258		258
259	`float min_len=FLT_MAX;`	259	`float min_len=FLT_MAX;`
260	`int min_k = -1;`	260	`int min_k = -1;`
261	`for(size_t k = 0; k < vpairs.size(); ++k){`	261	`for(size_t k = 0; k < vpairs.size(); ++k){`
262	`int i = vpairs[k].first;`	262	`int i = vpairs[k].first;`
263	`int j = vpairs[k].second;`	263	`int j = vpairs[k].second;`
264	`float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));`	264	`float len = sqr_length(m.pos(verts[i]) - m.pos(verts[j]));`
265		265
266	`if(len<min_len){`	266	`if(len<min_len){`
267	`min_len = len;`	267	`min_len = len;`
268	`min_k = k;`	268	`min_k = k;`
269	`}`	269	`}`
270	`}`	270	`}`
271	`assert(min_k != -1);`	271	`assert(min_k != -1);`
272		272
273	`{`	273	`{`
274	`// Split faces along edge whose midpoint is closest to isovalue`	274	`// Split faces along edge whose midpoint is closest to isovalue`
275	`int i = vpairs[min_k].first;`	275	`int i = vpairs[min_k].first;`
276	`int j = vpairs[min_k].second;`	276	`int j = vpairs[min_k].second;`
277	`m.split_face_by_edge(*f, verts[i], verts[j]);`	277	`m.split_face_by_edge(*f, verts[i], verts[j]);`
278	`}`	278	`}`
279	`++work;`	279	`++work;`
280		280
281	`}`	281	`}`
282	`}`	282	`}`
283	`while(work);`	283	`while(work);`
284	`}`	284	`}`
285		285
286	`void triangulate_face_by_edge_split(Manifold& m, FaceID f)`	286	`void triangulate_face_by_edge_split(Manifold& m, FaceID f)`
287	`{`	287	`{`
288	`if(no_edges(m, f)<=3)`	288	`if(no_edges(m, f)<=3)`
289	`return;`	289	`return;`
290		290
291	`Walker w = m.walker(f);`	291	`Walker w = m.walker(f);`
292		292
293	`// as long as f is not a triangle`	293	`// as long as f is not a triangle`
294	`while(w.next().next().next().halfedge() != w.halfedge()){`	294	`while(w.next().next().next().halfedge() != w.halfedge()){`
295	`// assert that face has at least 3 edges`	295	`// assert that face has at least 3 edges`
296	`// f is split into triangle from first three vertices, and becomes remaining polygon in next iteration`	296	`// f is split into triangle from first three vertices, and becomes remaining polygon in next iteration`
297	`assert(w.next().next().halfedge() != w.halfedge());`	297	`assert(w.next().next().halfedge() != w.halfedge());`
298	`VertexID v0 = w.vertex();`	298	`VertexID v0 = w.vertex();`
299	`VertexID v1 = w.next().next().vertex();`	299	`VertexID v1 = w.next().next().vertex();`
300	`FaceID f_old = f;`	300	`FaceID f_old = f;`
301	`if(v0 != v1 && !connected(m, v0, v1))`	301	`if(v0 != v1 && !connected(m, v0, v1))`
302	`f = m.split_face_by_edge(f, v0, v1);`	302	`f = m.split_face_by_edge(f, v0, v1);`
303	`if(f == f_old)`	303	`if(f == f_old)`
304	`return;`	304	`return;`
305	`}`	305	`}`
306	`}`	306	`}`
307	`}`	307	`}`
308		308

Subversion Repositories gelsvn

(root)/branches/cpp11-devel/src/HMesh/triangulate.cpp – Rev 607 → 630