Subversion Repositories gelsvn

Rev

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

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