Subversion Repositories gelsvn

Rev

Rev 61 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
61 jab 1 
#include <cfloat>
2 
#include <algorithm>
3 
#include <fstream>
4 
#include <iostream>
5 
#include <list>
6 
 
7 
#include "CGLA/Mat2x3f.h"
8 
#include "CGLA/Mat2x2f.h"
9 
#include "HMesh/Manifold.h"
10 
#include "HMesh/VertexCirculator.h"
11 
#include "HMesh/FaceCirculator.h"
12 
#include "HMeshUtil/triangulate.h"
13 
#include "HMeshUtil/build_manifold.h"
14 
#include "HMeshUtil/x3d_save.h"
15 
#include "HMeshUtil/obj_save.h"
16 
#include "HMeshUtil/mesh_optimization.h"
17 
#include "HMeshUtil/smooth.h"
18 
 
19 
#include "Geometry/RGrid.h"
20 
#include "Geometry/TrilinFilter.h"
21 
#include "Geometry/GradientFilter.h"
22 
#include "Geometry/Neighbours.h"
23 
#include "Util/HashTable.h"
24 
#include "Util/HashKey.h"
25 
 
26 
#include "fair_polygonize.h"
27 
 
64 jab 28 
namespace Geometry
61 jab 29 
{
30 
	using namespace HMesh;
31 
	using namespace HMeshUtil;
32 
	using namespace CMP;
33 
	using namespace CGLA;
34 
	using namespace std;
35 
 
36 
	namespace
37 
	{
38 
 
39 
		const Vec3i edges[6][4] = {
40 
			{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},
41 
			{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},
42 
			{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},
43 
			{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},
44 
			{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},
45 
			{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};
46 
 
47 
		int dir_to_face(const Vec3i& dir)
48 
		{
49 
			assert(dot(dir,dir)==1);
50 
 
51 
			if(dir[0] != 0)
52 
				return (dir[0]<0) ? 0 : 1;
53 
			if(dir[1] != 0)
54 
				return (dir[1]<0) ? 2 : 3;
55 
			if(dir[2] != 0)
56 
				return (dir[2]<0) ? 4 : 5;
57 
		}
58 
 
59 
 
60 
		int dir_to_edge(int face, const Vec3i& edge)
61 
		{
62 
			for(int i=0;i<4;++i)
63 
				if(edges[face][i] == edge) return i;
64 
			assert(0);
65 
			return -1;
66 
		}
67 
 
68 
		struct CubeFace
69 
		{
70 
			bool used;
71 
			VertexIter vert;
72 
			HalfEdgeIter he[4];
73 
 
74 
			CubeFace();
75 
 
76 
			HalfEdgeIter& get_he(Manifold* m, int i)
77 
			{
78 
				if(he[i] == NULL_HALFEDGE_ITER)
79 
					he[i] = m->create_halfedge();
80 
				return he[i];
81 
			}
82 
		};
83 
 
84 
		CubeFace::CubeFace(): used(false), vert(NULL_VERTEX_ITER)
85 
		{
86 
			for(int i=0;i<4;++i)
87 
				he[i] = NULL_HALFEDGE_ITER;
88 
		}
89 
 
90 
 
91 
		struct Cube
92 
		{
93 
			Vec3i pos;
94 
			CubeFace faces[6];
95 
 
96 
			Cube(const Vec3i& _pos): pos(_pos) {}
97 
		};
98 
 
99 
 
100 
 
101 
		template<class T>
102 
		class Mesher
103 
		{
104 
			const RGrid<T>* const voxel_grid;
105 
			Manifold* m;
106 
			HashTable<HashKey3usi,Cube*> cube_hash;
107 
			list<Cube> cube_table;
108 
			std::vector<VertexIter> mul_vertices;
109 
 
110 
		public:
111 
 
112 
			Mesher(const RGrid<T>* _voxel_grid,
113 
						 Manifold* _m,
114 
						 float iso, bool inv_cont);
115 
 
116 
			void process_cubes();
117 
 
118 
			void create_faces();
119 
 
120 
			void triangulate(float iso);
121 
 
122 
			void remove_duplicates();
123 
 
124 
			void push_vertices(float iso);
125 
 
126 
		};
127 
 
128 
		template<class T>
129 
		inline bool comp(bool d, T a, T b)
130 
		{
131 
			return d? (a<b) : (b<a);
132 
		}
133 
 
134 
		template<class T>
135 
		Mesher<T>::Mesher(const RGrid<T>* _voxel_grid,
136 
											Manifold* _m, float iso, bool inv_cont):
137 
			voxel_grid(_voxel_grid),
138 
			m(_m)
139 
		{
140 
			for(int k=0;k<voxel_grid->get_dims()[ 2 ];++k)
141 
				for(int j=0;j<voxel_grid->get_dims()[ 1 ];++j)
142 
					for(int i=0;i<voxel_grid->get_dims()[ 0 ];++i)
143 
						{
144 
							Vec3i pi0(i,j,k);
145 
							float val0 = (*voxel_grid)[pi0];
146 
							if(comp(inv_cont,val0,iso))
147 
								{
148 
									bool interface_cell = false;
149 
									vector<VertexIter> vertices;
150 
									Cube* cube;
151 
									for(int l=0;l<6;++l)
152 
										{
153 
											Vec3i pi1(pi0);
154 
											pi1 += N6i[l];
155 
											bool indom = voxel_grid->in_domain(pi1);
156 
											if(!indom ||
157 
												 !comp(inv_cont,float((*voxel_grid)[pi1]),iso))
158 
												{
159 
													float val1 = iso;
160 
													if(indom)
161 
														val1 = (*voxel_grid)[pi1];
162 
 
163 
													if(!interface_cell)
164 
														{
165 
															interface_cell = true;
166 
															Cube** c_ptr;
167 
															cube_hash.create_entry(HashKey3usi(pi0), c_ptr);
168 
															Cube c(pi0);
169 
															cube_table.push_back(c);
170 
															cube = (*c_ptr) = &cube_table.back();
171 
														}
172 
													int face_idx  = dir_to_face(N6i[l]);
173 
													CubeFace& face = cube->faces[face_idx];
174 
													face.used = true;
175 
													float t= (iso-val0)/(val1-val0);
176 
													Vec3f p = Vec3f(pi0)*(1.0f-t)+Vec3f(pi1)*t;
177 
													face.vert = m->create_vertex(p);
178 
// 													face.vert = m->create_vertex(Vec3f(pi0)+
179 
//  																											 N6f[l]/2.0f);
180 
													face.vert->touched = reinterpret_cast<int>(cube);
181 
												}
182 
										}
183 
								}
184 
						}
185 
		}
186 
 
187 
 
188 
		template<class T>
189 
		void Mesher<T>::process_cubes()
190 
		{
191 
			typename list<Cube>::iterator cube_iter = cube_table.begin();
192 
			for(;cube_iter != cube_table.end(); ++cube_iter)
193 
				{
194 
					Cube& cube = *cube_iter;
195 
					for(int i=0;i<6;++i)
196 
						if(cube.faces[i].used)
197 
							{
198 
								Vec3i pos = cube.pos;
199 
								Vec3i dir = N6i[i];
200 
								for(int e=0;e<4;++e)
201 
									{
202 
										Vec3i edge = edges[i][e];
203 
										Vec3i to_edge = cross(edge,dir);
204 
										CubeFace& face = cube.faces[i];
205 
 
206 
										Cube* cube_n;
207 
										Vec3i dir_n;
208 
										Vec3i edge_n;
209 
 
210 
 
211 
										// 											if(cube_hash.find_entry(HashKey3usi(pos + dir + to_edge),
212 
										// 																							cube_n))
213 
										// 												{
214 
										// 													assert(sqr_length(dir + to_edge)==2);
215 
										// 													dir_n = - to_edge;
216 
										// 													edge_n = - edge;
217 
										// 												}
218 
										// 											else if(cube_hash.find_entry(HashKey3usi(pos + to_edge),
219 
										// 																									 cube_n))
220 
										// 												{
221 
										// 													assert(sqr_length(to_edge)==1);
222 
										// 													dir_n = dir;
223 
										// 													edge_n = - edge;
224 
										// 												}
225 
										// 											else
226 
										// 												{
227 
										// 													dir_n = to_edge;
228 
										// 													edge_n = - edge;
229 
										// 													cube_n = &cube;
230 
										// 												}
231 
 
232 
										Cube* dummy;
233 
										if(cube.faces[dir_to_face(to_edge)].used)
234 
											{
235 
												dir_n = to_edge;
236 
												edge_n = - edge;
237 
												cube_n = &cube;
238 
											}
239 
										else if(!cube_hash.find_entry(HashKey3usi(pos+dir+to_edge),
240 
																									dummy))
241 
											{
242 
												assert(sqr_length(to_edge)==1);
243 
												dir_n = dir;
244 
												edge_n = - edge;
245 
												cube_hash.find_entry(HashKey3usi(pos + to_edge),
246 
																						 cube_n);
247 
											}
248 
										else
249 
											{
250 
												dir_n = - to_edge;
251 
												edge_n = - edge;
252 
												cube_n = dummy;
253 
											}
254 
 
255 
 
256 
										int i_n = dir_to_face(dir_n);
257 
										CubeFace& face_n = cube_n->faces[i_n];
258 
										int e_n = dir_to_edge(i_n, edge_n);
259 
										assert(face_n.used);
260 
 
261 
										HalfEdgeIter& he = face.get_he(m,e);
262 
										face.vert->he = he;
263 
										he->vert = face_n.vert;
264 
 
265 
										HalfEdgeIter& he_n = face_n.get_he(m, (e_n+3)%4);
266 
										link(he, he_n);
267 
 
268 
										HalfEdgeIter& he_opp = face_n.get_he(m, e_n);
269 
										face_n.vert->he = he_opp;
270 
										glue(he, he_opp);
271 
									}
272 
							}
273 
 
274 
				}
275 
		}
276 
 
277 
		template<class T>
278 
		void Mesher<T>::create_faces()
279 
		{
280 
			HalfEdgeIter he0 = m->halfedges_begin();
281 
			while(he0 != m->halfedges_end())
282 
				{
283 
					if(he0->face == NULL_FACE_ITER)
284 
						{
285 
							FaceIter fi = m->create_face();
286 
							fi->last = he0;
287 
 
288 
							assert(he0->face == NULL_FACE_ITER);
289 
							assert(he0->prev != NULL_HALFEDGE_ITER);
290 
							assert(he0->next != NULL_HALFEDGE_ITER);
291 
							assert(he0->vert != NULL_VERTEX_ITER);
292 
							assert(he0->vert->he != NULL_HALFEDGE_ITER);
293 
							assert(he0->opp != NULL_HALFEDGE_ITER);
294 
 
295 
							he0->face = fi;
296 
 
297 
							HalfEdgeIter he = he0->next;
298 
							while(he != he0)
299 
								{
300 
									assert(he->face == NULL_FACE_ITER);
301 
									assert(he->prev != NULL_HALFEDGE_ITER);
302 
									assert(he->next != NULL_HALFEDGE_ITER);
303 
									assert(he->vert != NULL_VERTEX_ITER);
304 
									assert(he->vert->he != NULL_HALFEDGE_ITER);
305 
									assert(he->opp != NULL_HALFEDGE_ITER);
306 
									he->face = fi;
307 
									he = he->next;
308 
								}
309 
						}
310 
					++he0;
311 
				}
312 
		}
313 
 
314 
 
315 
		template<class T>
316 
		void Mesher<T>::triangulate(float iso)
317 
		{
318 
#ifndef NDEBUG
319 
			cout << "Initial mesh valid : " << m->is_valid() << endl;
320 
#endif
321 
			Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(voxel_grid);
322 
 
323 
			int work;
324 
			do
325 
				{
326 
					work = 0;
327 
					for(FaceIter fi =m->faces_begin(); fi != m->faces_end(); ++fi)
328 
						{
329 
							vector<VertexIter> verts;
330 
							for(FaceCirculator fc(fi); !fc.end(); ++fc)
331 
								verts.push_back(fc.get_vertex());
332 
 
333 
							if(verts.size() == 3) continue;
334 
 
335 
							vector<pair<int,int> > vpairs;
336 
							const int N = verts.size();
337 
							for(int i=0;i<N-2;++i)
338 
								for(int j=i+2;j<N; ++j)
339 
									if(!is_connected(verts[i], verts[j]))
340 
										vpairs.push_back(pair<int,int>(i,j));
64 jab 341 
							if(vpairs.empty() && N > 3)
342 
							{
343 
									cout << __FILE__ << __LINE__ << " " << N << endl;
344 
									for(int i=0;i<N-2;++i)
345 
											for(int j=i+2;j<N; ++j)
346 
											{
347 
													cout << "i=" << i << " j=" << j << " ";
348 
													cout << &*verts[i] << " " << &*verts[j] << " ";
349 
													cout << is_connected(verts[i], verts[j]) << endl;
350 
											}
351 
							}
352 
 
61 jab 353 
							if(!vpairs.empty())
354 
								{
355 
									float min_val=FLT_MAX;
356 
									int min_k = -1;
357 
									for(int k=0;k<vpairs.size(); ++k)
358 
										{
359 
											int i = vpairs[k].first;
360 
											int j = vpairs[k].second;
361 
											Vec3f centre =
362 
												(verts[i]->pos +
363 
												 verts[j]->pos)/2.0f;
364 
 
365 
											float val;
366 
											if(grid_interp.in_domain(centre))
367 
												val = fabs(grid_interp(centre)-iso);
368 
											else
369 
												val = 0.0f;
370 
											if(val<min_val)
371 
												{
372 
													min_val = val;
373 
													min_k = k;
374 
												}
375 
										}
376 
									assert(min_k != -1);
377 
									int i = vpairs[min_k].first;
378 
									int j = vpairs[min_k].second;
379 
									m->split_face(fi, verts[i], verts[j]);
380 
									++work;
381 
								}
382 
						}
383 
				}
384 
			while(work);
385 
 
386 
#ifndef NDEBUG
387 
			cout << "Valid after triangulation : " << m->is_valid() << endl;
388 
#endif
389 
		}
390 
 
391 
		template<class T>
392 
		void Mesher<T>::remove_duplicates()
393 
		{
394 
			// BUG --- we never do this loop more than one time.
395 
			// does it matter.
396 
			VertexIter vi = m->vertices_begin();
397 
			bool did_work = false;
398 
			do
399 
				{
400 
					did_work = false;
401 
					while(vi != m->vertices_end())
402 
						{
403 
							bool did = false;
404 
							VertexCirculator vc(vi);
405 
							int k=0;
406 
							for(;!vc.end();++vc)
407 
								{
408 
									assert(++k<1000);
409 
									HalfEdgeIter he = vc.get_halfedge();
410 
									VertexIter n = vc.get_vertex();
411 
									if(n->touched == vi->touched && m->collapse_precond(vi,he))
412 
										{
413 
											VertexIter vi2 = vi;
414 
											++vi;
415 
											m->collapse_halfedge(vi2,he,true);
416 
											did = did_work = true;
417 
											break;
418 
										}
419 
								}
420 
							if(!did) ++vi;
421 
						}
422 
 
423 
				}
424 
			while(did_work);
425 
#ifndef NDEBUG
426 
			cout << "Valid after removal of excess vertices : " << m->is_valid() << endl;
427 
#endif
428 
		}
429 
 
430 
		template<class T>
431 
		void Mesher<T>::push_vertices(float iso)
432 
		{
433 
			GradientFilter<RGrid<T> > grad(voxel_grid);
434 
			TrilinFilter<GradientFilter<RGrid<T> > > ginter(&grad);
435 
			TrilinFilter<RGrid<T> > inter(voxel_grid);
436 
 
437 
			cout << "Pushing vertices" << endl;
438 
			for(int i=0;i<4;++i)
439 
				{
440 
					cout << "." << flush;
441 
					for(VertexIter vi= m->vertices_begin();
442 
							vi != m->vertices_end(); ++vi)
443 
						{
444 
							Vec3f p = vi->get_pos();
445 
							if(ginter.in_domain(p))
446 
								{
447 
									Vec3f g = ginter(p);
448 
									float s = sqr_length(g);
449 
									if(s>0.0001)
450 
										{
451 
											float v = inter(p)-iso;
452 
											vi->set_pos(p-g*v/s);
453 
										}
454 
								}
455 
						}
456 
				}
457 
		}
458 
 
459 
	}
460 
 
461 
	template<typename T>
462 
	void fair_polygonize(const RGrid<T>& voxel_grid,
463 
											 Manifold& mani,
464 
											 float iso,
465 
											 bool invert_contour)
466 
	{
467 
		Mesher<T> m(&voxel_grid, &mani, iso, invert_contour);
468 
		m.process_cubes();
469 
		m.create_faces();
64 jab 470 
		m.triangulate(iso);
61 jab 471 
 		m.remove_duplicates();
472 
		//m.push_vertices(iso);
473 
	}
474 
 
475 
	template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,
476 
																							 Manifold&, float, bool);
477 
	template void fair_polygonize<float>(const RGrid<float>&,
478 
																			 Manifold&, float, bool);
479 
}