WebSVN – gelsvn – Diff – /trunk/src/Geometry/fair_polygonize.cpp


#include <cfloat>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <list>
 
#include "CGLA/Mat2x3f.h"
#include "CGLA/Mat2x2f.h"
#include "HMesh/Manifold.h"
#include "HMesh/VertexCirculator.h"
#include "HMesh/FaceCirculator.h"
#include "HMeshUtil/triangulate.h"
#include "HMeshUtil/build_manifold.h"
#include "HMeshUtil/x3d_save.h"
#include "HMeshUtil/obj_save.h"
#include "HMeshUtil/mesh_optimization.h"
#include "HMeshUtil/smooth.h"
 
#include "Geometry/RGrid.h"
#include "Geometry/TrilinFilter.h"
#include "Geometry/GradientFilter.h"
#include "Geometry/Neighbours.h"
#include "Util/HashTable.h"
#include "Util/HashKey.h"
 
#include "fair_polygonize.h"
 
namespace Geometry
{
	using namespace HMesh;
	using namespace HMeshUtil;
	using namespace CMP;
 
	using namespace CGLA;
	using namespace std;
 
	namespace 
	{
 
		const Vec3i edges[6][4] = {
			{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},
			{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},
			{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},
			{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},
			{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},
			{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};
 
		int dir_to_face(const Vec3i& dir)
		{
			assert(dot(dir,dir)==1);
 
			if(dir[0] != 0)
				return (dir[0]<0) ? 0 : 1;
			if(dir[1] != 0)
				return (dir[1]<0) ? 2 : 3;
			if(dir[2] != 0)
				return (dir[2]<0) ? 4 : 5;
		}
		
 
		int dir_to_edge(int face, const Vec3i& edge)
		{
			for(int i=0;i<4;++i)
				if(edges[face][i] == edge) return i;
			assert(0);
			return -1;
		}
			
		struct CubeFace
		{
			bool used;
			VertexIter vert;
			HalfEdgeIter he[4];
				
			CubeFace();
 
			HalfEdgeIter& get_he(Manifold* m, int i)
			{
				if(he[i] == NULL_HALFEDGE_ITER)
					he[i] = m->create_halfedge();
				return he[i];
			}
		};
 
		CubeFace::CubeFace(): used(false), vert(NULL_VERTEX_ITER)
		{
			for(int i=0;i<4;++i)
				he[i] = NULL_HALFEDGE_ITER;
		}
		
 
		struct Cube
		{
			Vec3i pos;
			CubeFace faces[6];
			
			Cube(const Vec3i& _pos): pos(_pos) {}
		};
 
 
 
		template<class T>
		class Mesher
		{
			const RGrid<T>* const voxel_grid;
			Manifold* m;		
			HashTable<HashKey3usi,Cube*> cube_hash;
			list<Cube> cube_table;
			std::vector<VertexIter> mul_vertices;
 
		public:
			
			Mesher(const RGrid<T>* _voxel_grid, 
						 Manifold* _m, 
						 float iso, bool inv_cont);
 
			void process_cubes();
 
			void create_faces();
 
			void triangulate(float iso);
 
			void remove_duplicates();
 
			void push_vertices(float iso);
 
		};
	
		template<class T>
		inline bool comp(bool d, T a, T b)
		{
			return d? (a<b) : (b<a); 
		}
		
		template<class T>
		Mesher<T>::Mesher(const RGrid<T>* _voxel_grid, 
											Manifold* _m, float iso, bool inv_cont): 
			voxel_grid(_voxel_grid),
			m(_m)
		{
			for(int k=0;k<voxel_grid->get_dims()[ 2 ];++k)
				for(int j=0;j<voxel_grid->get_dims()[ 1 ];++j)
					for(int i=0;i<voxel_grid->get_dims()[ 0 ];++i)
						{
							Vec3i pi0(i,j,k);
							float val0 = (*voxel_grid)[pi0];
							if(comp(inv_cont,val0,iso))
								{
									bool interface_cell = false;
									vector<VertexIter> vertices;
									Cube* cube;
									for(int l=0;l<6;++l)
										{
											Vec3i pi1(pi0);
											pi1 += N6i[l];
											bool indom = voxel_grid->in_domain(pi1);
											if(!indom || 
												 !comp(inv_cont,float((*voxel_grid)[pi1]),iso))
												{
													float val1 = iso;
													if(indom)
														val1 = (*voxel_grid)[pi1];
													
													if(!interface_cell)
														{
															interface_cell = true;
															Cube** c_ptr;
															cube_hash.create_entry(HashKey3usi(pi0), c_ptr);
															Cube c(pi0);
															cube_table.push_back(c);
															cube = (*c_ptr) = &cube_table.back();
														}
													int face_idx  = dir_to_face(N6i[l]);
													CubeFace& face = cube->faces[face_idx];
													face.used = true;
													float t= (iso-val0)/(val1-val0);
													Vec3f p = Vec3f(pi0)*(1.0f-t)+Vec3f(pi1)*t;
													face.vert = m->create_vertex(p);
// 													face.vert = m->create_vertex(Vec3f(pi0)+
//  																											 N6f[l]/2.0f);
													face.vert->touched = reinterpret_cast<int>(cube);
												}
										}
								}
						}		
		}
 
 
		template<class T>
		void Mesher<T>::process_cubes()
		{
			typename list<Cube>::iterator cube_iter = cube_table.begin();
			for(;cube_iter != cube_table.end(); ++cube_iter)
				{
					Cube& cube = *cube_iter;
					for(int i=0;i<6;++i)
						if(cube.faces[i].used)
							{
								Vec3i pos = cube.pos;
								Vec3i dir = N6i[i];
								for(int e=0;e<4;++e)
									{
										Vec3i edge = edges[i][e];
										Vec3i to_edge = cross(edge,dir);
										CubeFace& face = cube.faces[i];
 
										Cube* cube_n;
										Vec3i dir_n;
										Vec3i edge_n;
										
											
										// 											if(cube_hash.find_entry(HashKey3usi(pos + dir + to_edge), 
										// 																							cube_n))
										// 												{
										// 													assert(sqr_length(dir + to_edge)==2);
										// 													dir_n = - to_edge;
										// 													edge_n = - edge;
										// 												}
										// 											else if(cube_hash.find_entry(HashKey3usi(pos + to_edge), 
										// 																									 cube_n))
										// 												{
										// 													assert(sqr_length(to_edge)==1);
										// 													dir_n = dir;
										// 													edge_n = - edge;
										// 												}
										// 											else
										// 												{
										// 													dir_n = to_edge;
										// 													edge_n = - edge;
										// 													cube_n = &cube;
										// 												}
 
										Cube* dummy;
										if(cube.faces[dir_to_face(to_edge)].used)
											{
												dir_n = to_edge;
												edge_n = - edge;
												cube_n = &cube;
											}
										else if(!cube_hash.find_entry(HashKey3usi(pos+dir+to_edge),
																									dummy))
											{
												assert(sqr_length(to_edge)==1);
												dir_n = dir;
												edge_n = - edge;
												cube_hash.find_entry(HashKey3usi(pos + to_edge), 
																						 cube_n);
											}
										else 
											{
												dir_n = - to_edge;
												edge_n = - edge;
												cube_n = dummy;
											}
 
 
										int i_n = dir_to_face(dir_n);
										CubeFace& face_n = cube_n->faces[i_n];
										int e_n = dir_to_edge(i_n, edge_n);
										assert(face_n.used);
										
										HalfEdgeIter& he = face.get_he(m,e);
										face.vert->he = he;
										he->vert = face_n.vert;
 
										HalfEdgeIter& he_n = face_n.get_he(m, (e_n+3)%4);
										link(he, he_n);
 
										HalfEdgeIter& he_opp = face_n.get_he(m, e_n);
										face_n.vert->he = he_opp;
										glue(he, he_opp);
									}
							}
 
				}
		}
 
		template<class T>
		void Mesher<T>::create_faces()
		{
			HalfEdgeIter he0 = m->halfedges_begin();
			while(he0 != m->halfedges_end())
				{
					if(he0->face == NULL_FACE_ITER)
						{
							FaceIter fi = m->create_face();
							fi->last = he0;
 
							assert(he0->face == NULL_FACE_ITER);
							assert(he0->prev != NULL_HALFEDGE_ITER);
							assert(he0->next != NULL_HALFEDGE_ITER);
							assert(he0->vert != NULL_VERTEX_ITER);
							assert(he0->vert->he != NULL_HALFEDGE_ITER);
							assert(he0->opp != NULL_HALFEDGE_ITER);
							
							he0->face = fi;
							
							HalfEdgeIter he = he0->next;
							while(he != he0)
								{
									assert(he->face == NULL_FACE_ITER);
									assert(he->prev != NULL_HALFEDGE_ITER);
									assert(he->next != NULL_HALFEDGE_ITER);
									assert(he->vert != NULL_VERTEX_ITER);
									assert(he->vert->he != NULL_HALFEDGE_ITER);
									assert(he->opp != NULL_HALFEDGE_ITER);
									he->face = fi;
									he = he->next;
								}
						}
					++he0;
				}
		}
 
 
		template<class T>
		void Mesher<T>::triangulate(float iso)
		{
#ifndef NDEBUG
			cout << "Initial mesh valid : " << m->is_valid() << endl;
#endif
			Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(voxel_grid);
 
			int work;
			do
				{
					work = 0;
					for(FaceIter fi =m->faces_begin(); fi != m->faces_end(); ++fi)
						{
							vector<VertexIter> verts;
							for(FaceCirculator fc(fi); !fc.end(); ++fc)
								verts.push_back(fc.get_vertex());
								
							if(verts.size() == 3) continue;
								
							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(!is_connected(verts[i], verts[j]))
										vpairs.push_back(pair<int,int>(i,j));
							if(vpairs.empty() && N > 3)
							{
									cout << __FILE__ << __LINE__ << " " << N << endl;
									for(int i=0;i<N-2;++i)
											for(int j=i+2;j<N; ++j)
											{
													cout << "i=" << i << " j=" << j << " ";
													cout << &*verts[i] << " " << &*verts[j] << " ";
													cout << is_connected(verts[i], verts[j]) << endl;
											}
							}
 
							if(!vpairs.empty())
								{
									float min_val=FLT_MAX;
									int min_k = -1;
									for(int k=0;k<vpairs.size(); ++k)
										{
											int i = vpairs[k].first;
											int j = vpairs[k].second;
											Vec3f centre = 
												(verts[i]->pos + 
												 verts[j]->pos)/2.0f;
 
											float val;
											if(grid_interp.in_domain(centre))
												val = fabs(grid_interp(centre)-iso);
											else
												val = 0.0f;
											if(val<min_val)
												{
													min_val = val;
													min_k = k;
												}
										}
									assert(min_k != -1);
									int i = vpairs[min_k].first;
									int j = vpairs[min_k].second;
									m->split_face(fi, verts[i], verts[j]);
									++work;
								}
						}
				}
			while(work);
 
#ifndef NDEBUG
			cout << "Valid after triangulation : " << m->is_valid() << endl;
#endif
		}
 
		template<class T>
		void Mesher<T>::remove_duplicates()
		{
			// BUG --- we never do this loop more than one time.
			// does it matter.
			VertexIter vi = m->vertices_begin();
			bool did_work = false;
			do
				{
					did_work = false;
					while(vi != m->vertices_end())
						{
							bool did = false;
							VertexCirculator vc(vi);
							int k=0;
							for(;!vc.end();++vc)
								{
									assert(++k<1000);
									HalfEdgeIter he = vc.get_halfedge();
									VertexIter n = vc.get_vertex();
									if(n->touched == vi->touched && m->collapse_precond(vi,he))
										{
											VertexIter vi2 = vi;
											++vi;
											m->collapse_halfedge(vi2,he,true);
											did = did_work = true;
											break;
										}
								}
							if(!did) ++vi;
						}
 
				}
			while(did_work);
#ifndef NDEBUG
			cout << "Valid after removal of excess vertices : " << m->is_valid() << endl;
#endif
		}
		
		template<class T>
		void Mesher<T>::push_vertices(float iso)
		{
			GradientFilter<RGrid<T> > grad(voxel_grid);
			TrilinFilter<GradientFilter<RGrid<T> > > ginter(&grad);
			TrilinFilter<RGrid<T> > inter(voxel_grid);
		
			cout << "Pushing vertices" << endl;
			for(int i=0;i<4;++i)
				{
					cout << "." << flush;
					for(VertexIter vi= m->vertices_begin();
							vi != m->vertices_end(); ++vi)
						{
							Vec3f p = vi->get_pos();
							if(ginter.in_domain(p))
								{
									Vec3f g = ginter(p);
									float s = sqr_length(g);
									if(s>0.0001)
										{
											float v = inter(p)-iso;
											vi->set_pos(p-g*v/s);
										}
								}
						}
				}
		}
 
	}
 
	template<typename T>
	void fair_polygonize(const RGrid<T>& voxel_grid, 
											 Manifold& mani, 
											 float iso,
											 bool invert_contour)
	{
		Mesher<T> m(&voxel_grid, &mani, iso, invert_contour);
		m.process_cubes();
		m.create_faces();
 
		m.triangulate(iso);
 		m.remove_duplicates();
		//m.push_vertices(iso);
	}
 
	template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,
																							 Manifold&, float, bool);
	template void fair_polygonize<float>(const RGrid<float>&,
																			 Manifold&, float, bool);
}
 

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

Subversion Repositories gelsvn

(root)/trunk/src/Geometry/fair_polygonize.cpp @ 78 – Rev 61 → 64