Subversion Repositories gelsvn

Rev

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