WebSVN – gelsvn – Diff – /trunk/src/HMesh/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 "HMesh/build_manifold.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 HMesh
{
	using namespace Geometry;
	using namespace Util;
	using namespace CGLA;
	using namespace std;
 
	namespace 
	{
 
		void compute_dual(Manifold& m, const vector<Vec3f>& face_positions)
		{
			// make sure every face knows its number
			m.enumerate_faces();
 
			vector<Vec3f> vertices(m.no_faces());
			vector<int> faces;
			vector<int> indices;
 
			// Create new vertices. Each face becomes a vertex whose positio
			// is the centre of the face
			int i=0;
			for(FaceIter f=m.faces_begin(); f!=m.faces_end(); ++f,++i)
				vertices[i] = face_positions[f->touched];
 
			// Create new faces. Each vertex is a new face with N=valency of
			// edges.
			i=0;
			for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v,++i)
				if(!is_boundary(v))
					{
						int no_edges = 0;
						vector<int> index_tmp;
							
						// Circulate around vertex
						VertexCirculator vc(v);
						for(; !vc.end(); ++vc)
							{
								FaceIter f = vc.get_face();
								if(f != NULL_FACE_ITER)
									{
										index_tmp.push_back(f->touched);
										++no_edges;
									}
							}
						// Push vertex indices for this face onto indices vector
							
						// The circulator moves around the face in a clockwise f
						// so we just reverse the ordering.
						indices.insert(indices.end(), index_tmp.rbegin(), index_tmp.rend());
						faces.push_back(no_edges);
					}
 
			// Clear the manifold before new geometry is inserted.
			m.clear();
 
			// And build
			build_manifold(m, vertices.size(), &vertices[0], faces.size(),
										 &faces[0],&indices[0]);
 
		}
 
		/*
			Vectors along edges. For each of six cube faces there are four such vectors
			since the faces of a cube are quads.
		*/
		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)}};
 
		/* Convert a direction vector to a cube face index. */
		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;
 
			return -1;
		}
		
		/* Convert a face index and an edge vector to an edge index. */
		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;
		}
 
		Vec3f cube_corner(const Vec3i& pos, int i, int e)
		{
			Vec3i to_face = N6i[i];
			Vec3i along_edge = edges[i][e];	
			Vec3i to_edge = cross(along_edge,to_face);
			return Vec3f(pos) + 0.5*Vec3f(to_face+to_edge+along_edge);
		}
			
		/* Cube faces are created whenever a voxel that is inside is adjacent
			 to a voxel that is outside. Thus the voxel can be seen as a small
			 box that either contains material or not. A face is created where
			 full and empty boxes meet. A vertex is placed on every box face.
				 
			 A halfedge is created for every edge of the (quadrilateral) face.
			 This halfedge connects the face to its neighbour. */
		struct CubeFace
		{
			bool used;
			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)
		{
			for(int i=0;i<4;++i)
				he[i] = NULL_HALFEDGE_ITER;
		}
		
 
		/* A cube consists of six faces and a position. The cube centre is
			 an integer (voxel grid) position. */
		struct Cube
		{
			Vec3i pos;
			CubeFace faces[6];
			
			Cube(const Vec3i& _pos): pos(_pos) {}
		};
 
 
		/* The mesher class does the actual work. It is not pretty and it
			 is not seen by the user. */
		template<class T>
		class CubeGrid
		{
			float iso;
			HashTable<HashKey3usi,Cube*> cube_hash;  // Hashtable of cube pointers
			list<Cube> cube_table;                   // Linked list containing
			// the actual cubes.
 
		public:
			
			CubeGrid(const RGrid<T>& _voxel_grid, // input voxel grid
							 float iso);                   // iso value
 
 
			HalfEdgeIter edge_to_glue(const RGrid<T>& voxel_grid,
																Cube& cube, int face_idx, int edge_idx);
					
			void cuberille_mesh(const RGrid<T>&, Manifold& m, vector<Vec3f>& face_positions);
		};
 
		template<class T>
		HalfEdgeIter CubeGrid<T>::edge_to_glue(const RGrid<T>& voxel_grid,
																					 Cube& cube, int face_idx, int edge_idx)
		{
			/*  Figure.
 
			0|B
			-+-
			A|C
															
			Explanation: ... pending
			*/
 
			Vec3i dir = N6i[face_idx];
			Vec3i along_edge = edges[face_idx][edge_idx];	 // Vector _along_ edge
			Vec3i to_edge = cross(along_edge,dir); // vector to edge
 
			Vec3i a_pos = cube.pos;
			Vec3i o_pos = a_pos + dir;
			Vec3i b_pos = o_pos + to_edge;
			Vec3i c_pos = a_pos + to_edge;
			
			if(!voxel_grid.in_domain(b_pos))
				return NULL_HALFEDGE_ITER;
 
			Cube* cube_n=0;
			Vec3i dir_n;   // direction from center to neighbour face
			Vec3i along_edge_n; // Edge vector for neighbour face
			if(!cube_hash.find_entry(HashKey3usi(b_pos),	cube_n))
				{
					if(cube.faces[dir_to_face(to_edge)].used)
						{
							dir_n = to_edge;
							along_edge_n = - along_edge; 
							cube_n = &cube;  
						}
					else
						{
							dir_n = dir;
							along_edge_n = - along_edge;
							cube_hash.find_entry(HashKey3usi(c_pos), cube_n);
							assert(cube_n->faces[dir_to_face(dir_n)].used);
						}									
				}
			else
				{
					float a_val = voxel_grid[a_pos];
					float b_val = voxel_grid[b_pos];
					float c_val = voxel_grid[c_pos];
					float o_val = voxel_grid[o_pos];
					bool connect = (a_val*b_val-c_val*o_val)/(a_val+b_val-c_val-o_val) >iso;
								
					if(connect || !cube.faces[dir_to_face(to_edge)].used)
						{
							dir_n = - to_edge;
							along_edge_n = - along_edge;
							assert(cube_n->faces[dir_to_face(dir_n)].used);
						}
					else
						{
							dir_n = to_edge;
							along_edge_n = - along_edge; 
							cube_n = &cube;  
						}
				}
			int i_n = dir_to_face(dir_n);
			int e_n = dir_to_edge(i_n, along_edge_n);
			return cube_n->faces[i_n].he[e_n];
		}
				
		template<class T>
		inline bool comp(bool d, T a, T b)
		{
			return d? (a<b) : (b<a); 
		}
		
		template<class T>
		CubeGrid<T>::CubeGrid(const RGrid<T>& voxel_grid, float _iso): iso(_iso)
		{
			/* Loop over all voxels */
			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 the voxel value is above the isovalue */
							if(val0>iso)
								{
									bool interface_cell = false;
									vector<VertexIter> vertices;
									Cube* cube;
 
									// For each neighbouring voxel.
									for(int l=0;l<6;++l)
										{
											// First make sure the voxel is inside the voxel grid.
											Vec3i pi1(pi0);
											pi1 += N6i[l];
											bool indom = voxel_grid.in_domain(pi1);
												
											/* If the neighbour is **not** in the grid or it 
												 is on the other side of the isovalue */
											if(indom && voxel_grid[pi1]<=iso)
												{
													/* Store the cube in a hash table. This code is
														 executed the first time we create a face
														 for the cube. In other words only cubes that 
														 actually have faces are stored in the hash */
													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();
														}
													// Now add the face to the cube.
													int face_idx  = dir_to_face(N6i[l]);
													CubeFace& face = cube->faces[face_idx];
													face.used = true;
 
												}
										}
								}
						}		
		}
 
 
		template<class T>
		void CubeGrid<T>::cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m, 
																		 vector<Vec3f>& face_positions)
		{
			face_positions.clear();
			Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&voxel_grid);
			typename list<Cube>::iterator cube_iter = cube_table.begin();
 
			/* For each cube in the table of cubes */
			cube_iter = cube_table.begin();
			for(;cube_iter != cube_table.end(); ++cube_iter)
				{
					/* For each of the six faces, if the face is used */
					Cube& cube = *cube_iter;
					for(int i=0;i<6;++i)
						if(cube.faces[i].used)
							{
								CubeFace& face = cube.faces[i];  // The face
								FaceIter hmesh_face = m.create_face();
								hmesh_face->touched = face_positions.size();
								Vec3i p0(cube.pos);
								Vec3i p1 = p0 + N6i[i];
								float v0 = voxel_grid[p0];
								float v1 = voxel_grid[p1];
								float t = (iso-v0)/(v1-v0);
								face_positions.push_back(Vec3f(p0)*(1.0f-t)+Vec3f(p1)*t);
								
								// Let e be one of the four edges of the face
								for(int e=0;e<4;++e)
									{
										HalfEdgeIter he = face.get_he(m,e);
										link(he,face.get_he(m,(e+1)%4));
										link(face.get_he(m,(e+3)%4),he);
										he->face = hmesh_face;
 
										HalfEdgeIter he_n = edge_to_glue(voxel_grid, cube, i, e);
										if(he_n != NULL_HALFEDGE_ITER)
											glue(he, he_n);
									}
								hmesh_face->last = face.get_he(m,0);
							}
 
				}
 
			for(HalfEdgeIter h = m.halfedges_begin(); h != m.halfedges_end(); ++h)
				{
					if (h->opp == NULL_HALFEDGE_ITER)
						{
							HalfEdgeIter ho = m.create_halfedge();
							glue(ho,h);
						}
				}
 
			cube_iter = cube_table.begin();
			for(;cube_iter != cube_table.end(); ++cube_iter)
				{
					//For each of the six faces, if the face is used 
					Cube& cube = *cube_iter;
					for(int i=0;i<6;++i)
						if(cube.faces[i].used)
							{
								CubeFace& face = cube.faces[i];  // The face
								// Let e be one of the four edges of the face
								for(int e=0;e<4;++e)
									if(face.he[e]->vert == NULL_VERTEX_ITER)
										{
											HalfEdgeIter last = face.he[e];
											HalfEdgeIter he = last;
											while(he->opp->face != NULL_FACE_ITER)
												{
													he = he->opp->prev;
													if(he == last) 
														break;
												}
											last = he;
																
											Vec3f pos = cube_corner(cube.pos, i, e);
											VertexIter vert = m.create_vertex(pos);
																
											do
												{
													he->vert = vert;
													he = he->next->opp;
												}
											while(he->face != NULL_FACE_ITER && he != last);
 
											if(he->face == NULL_FACE_ITER)
												{
													he->vert = vert;
													link(he,last->opp);
												}
											
											vert->he = last->opp;
										}
							}
				}
		}
	}
	
		/* A fairly complex algorithm which triangulates faces. It triangulates
			 by iteratively splitting faces. It always splits a face by picking the
			 edge whose midpoint is closest to the isovalue. */
					
		template<class T>
		void triangulate(const RGrid<T>& voxel_grid, Manifold& m, 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
				{
					// For every face.
					work = 0;
					for(FaceIter fi =m.faces_begin(); fi != m.faces_end(); ++fi)
						{
							// Create a vector of vertices.
							vector<VertexIter> verts;
							for(FaceCirculator fc(fi); !fc.end(); ++fc)
								verts.push_back(fc.get_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(!is_connected(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 absolute value of iso-x
								 where x is the interpolated value of the volume at the 
								 midpoint of the line segment connecting the two vertices.*/
											 
							float min_val=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;
									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);
										
							{
								// 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(fi, verts[i], verts[j]);
							}
							++work;
										
						}
				}
			while(work);
 
#ifndef NDEBUG
			cout << "Valid after triangulation : " << m.is_valid() << endl;
#endif
		}
 
	
	template<typename T>
	void fair_polygonize(const RGrid<T>& voxel_grid, 
											 Manifold& mani, 
											 float iso)
	{
		CubeGrid<T> cube_grid(voxel_grid, iso);
		vector<Vec3f> face_positions;
		cube_grid.cuberille_mesh(voxel_grid, mani, face_positions);
		cout << "Initial mesh valid : " << mani.is_valid() << endl;
		compute_dual(mani, face_positions);
	}
	
	template<typename T>
	void cuberille_polygonize(const RGrid<T>& voxel_grid, 
														Manifold& mani, 
														float iso,
														bool push)
	{
		CubeGrid<T> cube_grid(voxel_grid, iso);
		vector<Vec3f> face_positions;
		cube_grid.cuberille_mesh(voxel_grid, mani, face_positions);
		if(push)
			{
				for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
					{
						VertexCirculator vc(v);
						Vec3f p(0);
						int n=0;
						while(!vc.end())
							{
								if(vc.get_face() != NULL_FACE_ITER)
									{
										p += face_positions[vc.get_face()->touched];
										++n;
									}
								++vc;
							}
						v->pos = p/float(n);
					}
			}
	}
 
	template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,
																							 Manifold&, float);
	template void fair_polygonize<float>(const RGrid<float>&,
																			 Manifold&, float);
 
	template void cuberille_polygonize<unsigned char>(const RGrid<unsigned char>&,
																										Manifold&, float, bool);
	template void cuberille_polygonize<float>(const RGrid<float>&,
																						Manifold&, float, bool);
 
	template void triangulate(const RGrid<unsigned char>& voxel_grid, Manifold& m, float iso);
 
	template void triangulate(const RGrid<float>& voxel_grid, Manifold& m, float iso);
}
	
 
 
	
	
 

Rev 170	Rev 182
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 "HMesh/build_manifold.h"`	12	`#include "HMesh/build_manifold.h"`
13		13
14		14
15	`#include "Geometry/RGrid.h"`	15	`#include "Geometry/RGrid.h"`
16	`#include "Geometry/TrilinFilter.h"`	16	`#include "Geometry/TrilinFilter.h"`
17	`#include "Geometry/GradientFilter.h"`	17	`#include "Geometry/GradientFilter.h"`
18	`#include "Geometry/Neighbours.h"`	18	`#include "Geometry/Neighbours.h"`
19	`#include "Util/HashTable.h"`	19	`#include "Util/HashTable.h"`
20	`#include "Util/HashKey.h"`	20	`#include "Util/HashKey.h"`
21		21
22	`#include "fair_polygonize.h"`	22	`#include "fair_polygonize.h"`
23		23
24	`namespace HMesh`	24	`namespace HMesh`
25	`{`	25	`{`
26	`using namespace Geometry;`	26	`using namespace Geometry;`
27	`using namespace Util;`	27	`using namespace Util;`
28	`using namespace CGLA;`	28	`using namespace CGLA;`
29	`using namespace std;`	29	`using namespace std;`
30		30
31	`namespace`	31	`namespace`
32	`{`	32	`{`
33		33
34	`void compute_dual(Manifold& m, const vector<Vec3f>& face_positions)`	34	`void compute_dual(Manifold& m, const vector<Vec3f>& face_positions)`
35	`{`	35	`{`
36	`// make sure every face knows its number`	36	`// make sure every face knows its number`
37	`m.enumerate_faces();`	37	`m.enumerate_faces();`
38		38
39	`vector<Vec3f> vertices(m.no_faces());`	39	`vector<Vec3f> vertices(m.no_faces());`
40	`vector<int> faces;`	40	`vector<int> faces;`
41	`vector<int> indices;`	41	`vector<int> indices;`
42		42
43	`// Create new vertices. Each face becomes a vertex whose positio`	43	`// Create new vertices. Each face becomes a vertex whose positio`
44	`// is the centre of the face`	44	`// is the centre of the face`
45	`int i=0;`	45	`int i=0;`
46	`for(FaceIter f=m.faces_begin(); f!=m.faces_end(); ++f,++i)`	46	`for(FaceIter f=m.faces_begin(); f!=m.faces_end(); ++f,++i)`
47	`vertices[i] = face_positions[f->touched];`	47	`vertices[i] = face_positions[f->touched];`
48		48
49	`// Create new faces. Each vertex is a new face with N=valency of`	49	`// Create new faces. Each vertex is a new face with N=valency of`
50	`// edges.`	50	`// edges.`
51	`i=0;`	51	`i=0;`
52	`for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v,++i)`	52	`for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v,++i)`
53	`if(!is_boundary(v))`	53	`if(!is_boundary(v))`
54	`{`	54	`{`
55	`int no_edges = 0;`	55	`int no_edges = 0;`
56	`vector<int> index_tmp;`	56	`vector<int> index_tmp;`
57		57
58	`// Circulate around vertex`	58	`// Circulate around vertex`
59	`VertexCirculator vc(v);`	59	`VertexCirculator vc(v);`
60	`for(; !vc.end(); ++vc)`	60	`for(; !vc.end(); ++vc)`
61	`{`	61	`{`
62	`FaceIter f = vc.get_face();`	62	`FaceIter f = vc.get_face();`
63	`if(f != NULL_FACE_ITER)`	63	`if(f != NULL_FACE_ITER)`
64	`{`	64	`{`
65	`index_tmp.push_back(f->touched);`	65	`index_tmp.push_back(f->touched);`
66	`++no_edges;`	66	`++no_edges;`
67	`}`	67	`}`
68	`}`	68	`}`
69	`// Push vertex indices for this face onto indices vector`	69	`// Push vertex indices for this face onto indices vector`
70		70
71	`// The circulator moves around the face in a clockwise f`	71	`// The circulator moves around the face in a clockwise f`
72	`// so we just reverse the ordering.`	72	`// so we just reverse the ordering.`
73	`indices.insert(indices.end(), index_tmp.rbegin(), index_tmp.rend());`	73	`indices.insert(indices.end(), index_tmp.rbegin(), index_tmp.rend());`
74	`faces.push_back(no_edges);`	74	`faces.push_back(no_edges);`
75	`}`	75	`}`
76		76
77	`// Clear the manifold before new geometry is inserted.`	77	`// Clear the manifold before new geometry is inserted.`
78	`m.clear();`	78	`m.clear();`
79		79
80	`// And build`	80	`// And build`
81	`build_manifold(m, vertices.size(), &vertices[0], faces.size(),`	81	`build_manifold(m, vertices.size(), &vertices[0], faces.size(),`
82	`&faces[0],&indices[0]);`	82	`&faces[0],&indices[0]);`
83		83
84	`}`	84	`}`
85		85
86	`/*`	86	`/*`
87	`Vectors along edges. For each of six cube faces there are four such vectors`	87	`Vectors along edges. For each of six cube faces there are four such vectors`
88	`since the faces of a cube are quads.`	88	`since the faces of a cube are quads.`
89	`*/`	89	`*/`
90	`const Vec3i edges[6][4] = {`	90	`const Vec3i edges[6][4] = {`
91	`{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},`	91	`{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},`
92	`{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},`	92	`{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},`
93	`{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},`	93	`{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},`
94	`{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},`	94	`{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},`
95	`{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},`	95	`{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},`
96	`{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};`	96	`{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};`
97		97
98	`/* Convert a direction vector to a cube face index. */`	98	`/* Convert a direction vector to a cube face index. */`
99	`int dir_to_face(const Vec3i& dir)`	99	`int dir_to_face(const Vec3i& dir)`
100	`{`	100	`{`
101	`assert(dot(dir,dir)==1);`	101	`assert(dot(dir,dir)==1);`
102		102
103	`if(dir[0] != 0)`	103	`if(dir[0] != 0)`
104	`return (dir[0]<0) ? 0 : 1;`	104	`return (dir[0]<0) ? 0 : 1;`
105	`if(dir[1] != 0)`	105	`if(dir[1] != 0)`
106	`return (dir[1]<0) ? 2 : 3;`	106	`return (dir[1]<0) ? 2 : 3;`
107	`if(dir[2] != 0)`	107	`if(dir[2] != 0)`
108	`return (dir[2]<0) ? 4 : 5;`	108	`return (dir[2]<0) ? 4 : 5;`
109		109
110	`return -1;`	110	`return -1;`
111	`}`	111	`}`
112		112
113	`/* Convert a face index and an edge vector to an edge index. */`	113	`/* Convert a face index and an edge vector to an edge index. */`
114	`int dir_to_edge(int face, const Vec3i& edge)`	114	`int dir_to_edge(int face, const Vec3i& edge)`
115	`{`	115	`{`
116	`for(int i=0;i<4;++i)`	116	`for(int i=0;i<4;++i)`
117	`if(edges[face][i] == edge) return i;`	117	`if(edges[face][i] == edge) return i;`
118	`assert(0);`	118	`assert(0);`
119	`return -1;`	119	`return -1;`
120	`}`	120	`}`
121		121
122	`Vec3f cube_corner(const Vec3i& pos, int i, int e)`	122	`Vec3f cube_corner(const Vec3i& pos, int i, int e)`
123	`{`	123	`{`
124	`Vec3i to_face = N6i[i];`	124	`Vec3i to_face = N6i[i];`
125	`Vec3i along_edge = edges[i][e];`	125	`Vec3i along_edge = edges[i][e];`
126	`Vec3i to_edge = cross(along_edge,to_face);`	126	`Vec3i to_edge = cross(along_edge,to_face);`
127	`return Vec3f(pos) + 0.5*Vec3f(to_face+to_edge+along_edge);`	127	`return Vec3f(pos) + 0.5*Vec3f(to_face+to_edge+along_edge);`
128	`}`	128	`}`
129		129
130	`/* Cube faces are created whenever a voxel that is inside is adjacent`	130	`/* Cube faces are created whenever a voxel that is inside is adjacent`
131	`to a voxel that is outside. Thus the voxel can be seen as a small`	131	`to a voxel that is outside. Thus the voxel can be seen as a small`
132	`box that either contains material or not. A face is created where`	132	`box that either contains material or not. A face is created where`
133	`full and empty boxes meet. A vertex is placed on every box face.`	133	`full and empty boxes meet. A vertex is placed on every box face.`
134		134
135	`A halfedge is created for every edge of the (quadrilateral) face.`	135	`A halfedge is created for every edge of the (quadrilateral) face.`
136	`This halfedge connects the face to its neighbour. */`	136	`This halfedge connects the face to its neighbour. */`
137	`struct CubeFace`	137	`struct CubeFace`
138	`{`	138	`{`
139	`bool used;`	139	`bool used;`
140	`HalfEdgeIter he[4];`	140	`HalfEdgeIter he[4];`
141		141
142	`CubeFace();`	142	`CubeFace();`
143		143
144	`HalfEdgeIter& get_he(Manifold& m, int i)`	144	`HalfEdgeIter& get_he(Manifold& m, int i)`
145	`{`	145	`{`
146	`if(he[i] == NULL_HALFEDGE_ITER)`	146	`if(he[i] == NULL_HALFEDGE_ITER)`
147	`he[i] = m.create_halfedge();`	147	`he[i] = m.create_halfedge();`
148	`return he[i];`	148	`return he[i];`
149	`}`	149	`}`
150	`};`	150	`};`
151		151
152	`CubeFace::CubeFace(): used(false)`	152	`CubeFace::CubeFace(): used(false)`
153	`{`	153	`{`
154	`for(int i=0;i<4;++i)`	154	`for(int i=0;i<4;++i)`
155	`he[i] = NULL_HALFEDGE_ITER;`	155	`he[i] = NULL_HALFEDGE_ITER;`
156	`}`	156	`}`
157		157
158		158
159	`/* A cube consists of six faces and a position. The cube centre is`	159	`/* A cube consists of six faces and a position. The cube centre is`
160	`an integer (voxel grid) position. */`	160	`an integer (voxel grid) position. */`
161	`struct Cube`	161	`struct Cube`
162	`{`	162	`{`
163	`Vec3i pos;`	163	`Vec3i pos;`
164	`CubeFace faces[6];`	164	`CubeFace faces[6];`
165		165
166	`Cube(const Vec3i& _pos): pos(_pos) {}`	166	`Cube(const Vec3i& _pos): pos(_pos) {}`
167	`};`	167	`};`
168		168
169		169
170	`/* The mesher class does the actual work. It is not pretty and it`	170	`/* The mesher class does the actual work. It is not pretty and it`
171	`is not seen by the user. */`	171	`is not seen by the user. */`
172	`template<class T>`	172	`template<class T>`
173	`class CubeGrid`	173	`class CubeGrid`
174	`{`	174	`{`
175	`float iso;`	175	`float iso;`
176	`HashTable<HashKey3usi,Cube*> cube_hash; // Hashtable of cube pointers`	176	`HashTable<HashKey3usi,Cube*> cube_hash; // Hashtable of cube pointers`
177	`list<Cube> cube_table; // Linked list containing`	177	`list<Cube> cube_table; // Linked list containing`
178	`// the actual cubes.`	178	`// the actual cubes.`
179		179
180	`public:`	180	`public:`
181		181
182	`CubeGrid(const RGrid<T>& _voxel_grid, // input voxel grid`	182	`CubeGrid(const RGrid<T>& _voxel_grid, // input voxel grid`
183	`float iso); // iso value`	183	`float iso); // iso value`
184		184
185		185
186	`HalfEdgeIter edge_to_glue(const RGrid<T>& voxel_grid,`	186	`HalfEdgeIter edge_to_glue(const RGrid<T>& voxel_grid,`
187	`Cube& cube, int face_idx, int edge_idx);`	187	`Cube& cube, int face_idx, int edge_idx);`
188		188
189	`void cuberille_mesh(const RGrid<T>&, Manifold& m, vector<Vec3f>& face_positions);`	189	`void cuberille_mesh(const RGrid<T>&, Manifold& m, vector<Vec3f>& face_positions);`
190	`};`	190	`};`
191		191
192	`template<class T>`	192	`template<class T>`
193	`HalfEdgeIter CubeGrid<T>::edge_to_glue(const RGrid<T>& voxel_grid,`	193	`HalfEdgeIter CubeGrid<T>::edge_to_glue(const RGrid<T>& voxel_grid,`
194	`Cube& cube, int face_idx, int edge_idx)`	194	`Cube& cube, int face_idx, int edge_idx)`
195	`{`	195	`{`
196	`/* Figure.`	196	`/* Figure.`
197		197
198	`0\|B`	198	`0\|B`
199	`-+-`	199	`-+-`
200	`A\|C`	200	`A\|C`
201		201
202	`Explanation: ... pending`	202	`Explanation: ... pending`
203	`*/`	203	`*/`
204		204
205	`Vec3i dir = N6i[face_idx];`	205	`Vec3i dir = N6i[face_idx];`
206	`Vec3i along_edge = edges[face_idx][edge_idx]; // Vector _along_ edge`	206	`Vec3i along_edge = edges[face_idx][edge_idx]; // Vector _along_ edge`
207	`Vec3i to_edge = cross(along_edge,dir); // vector to edge`	207	`Vec3i to_edge = cross(along_edge,dir); // vector to edge`
208		208
209	`Vec3i a_pos = cube.pos;`	209	`Vec3i a_pos = cube.pos;`
210	`Vec3i o_pos = a_pos + dir;`	210	`Vec3i o_pos = a_pos + dir;`
211	`Vec3i b_pos = o_pos + to_edge;`	211	`Vec3i b_pos = o_pos + to_edge;`
212	`Vec3i c_pos = a_pos + to_edge;`	212	`Vec3i c_pos = a_pos + to_edge;`
213		213
214	`if(!voxel_grid.in_domain(b_pos))`	214	`if(!voxel_grid.in_domain(b_pos))`
215	`return NULL_HALFEDGE_ITER;`	215	`return NULL_HALFEDGE_ITER;`
216		216
217	`Cube* cube_n=0;`	217	`Cube* cube_n=0;`
218	`Vec3i dir_n; // direction from center to neighbour face`	218	`Vec3i dir_n; // direction from center to neighbour face`
219	`Vec3i along_edge_n; // Edge vector for neighbour face`	219	`Vec3i along_edge_n; // Edge vector for neighbour face`
220	`if(!cube_hash.find_entry(HashKey3usi(b_pos), cube_n))`	220	`if(!cube_hash.find_entry(HashKey3usi(b_pos), cube_n))`
221	`{`	221	`{`
222	`if(cube.faces[dir_to_face(to_edge)].used)`	222	`if(cube.faces[dir_to_face(to_edge)].used)`
223	`{`	223	`{`
224	`dir_n = to_edge;`	224	`dir_n = to_edge;`
225	`along_edge_n = - along_edge;`	225	`along_edge_n = - along_edge;`
226	`cube_n = &cube;`	226	`cube_n = &cube;`
227	`}`	227	`}`
228	`else`	228	`else`
229	`{`	229	`{`
230	`dir_n = dir;`	230	`dir_n = dir;`
231	`along_edge_n = - along_edge;`	231	`along_edge_n = - along_edge;`
232	`cube_hash.find_entry(HashKey3usi(c_pos), cube_n);`	232	`cube_hash.find_entry(HashKey3usi(c_pos), cube_n);`
233	`assert(cube_n->faces[dir_to_face(dir_n)].used);`	233	`assert(cube_n->faces[dir_to_face(dir_n)].used);`
234	`}`	234	`}`
235	`}`	235	`}`
236	`else`	236	`else`
237	`{`	237	`{`
238	`float a_val = voxel_grid[a_pos];`	238	`float a_val = voxel_grid[a_pos];`
239	`float b_val = voxel_grid[b_pos];`	239	`float b_val = voxel_grid[b_pos];`
240	`float c_val = voxel_grid[c_pos];`	240	`float c_val = voxel_grid[c_pos];`
241	`float o_val = voxel_grid[o_pos];`	241	`float o_val = voxel_grid[o_pos];`
242	`bool connect = (a_valb_val-c_valo_val)/(a_val+b_val-c_val-o_val) >iso;`	242	`bool connect = (a_valb_val-c_valo_val)/(a_val+b_val-c_val-o_val) >iso;`
243		243
244	`if(connect \|\| !cube.faces[dir_to_face(to_edge)].used)`	244	`if(connect \|\| !cube.faces[dir_to_face(to_edge)].used)`
245	`{`	245	`{`
246	`dir_n = - to_edge;`	246	`dir_n = - to_edge;`
247	`along_edge_n = - along_edge;`	247	`along_edge_n = - along_edge;`
248	`assert(cube_n->faces[dir_to_face(dir_n)].used);`	248	`assert(cube_n->faces[dir_to_face(dir_n)].used);`
249	`}`	249	`}`
250	`else`	250	`else`
251	`{`	251	`{`
252	`dir_n = to_edge;`	252	`dir_n = to_edge;`
253	`along_edge_n = - along_edge;`	253	`along_edge_n = - along_edge;`
254	`cube_n = &cube;`	254	`cube_n = &cube;`
255	`}`	255	`}`
256	`}`	256	`}`
257	`int i_n = dir_to_face(dir_n);`	257	`int i_n = dir_to_face(dir_n);`
258	`int e_n = dir_to_edge(i_n, along_edge_n);`	258	`int e_n = dir_to_edge(i_n, along_edge_n);`
259	`return cube_n->faces[i_n].he[e_n];`	259	`return cube_n->faces[i_n].he[e_n];`
260	`}`	260	`}`
261		261
262	`template<class T>`	262	`template<class T>`
263	`inline bool comp(bool d, T a, T b)`	263	`inline bool comp(bool d, T a, T b)`
264	`{`	264	`{`
265	`return d? (a<b) : (b<a);`	265	`return d? (a<b) : (b<a);`
266	`}`	266	`}`
267		267
268	`template<class T>`	268	`template<class T>`
269	`CubeGrid<T>::CubeGrid(const RGrid<T>& voxel_grid, float _iso): iso(_iso)`	269	`CubeGrid<T>::CubeGrid(const RGrid<T>& voxel_grid, float _iso): iso(_iso)`
270	`{`	270	`{`
271	`/* Loop over all voxels */`	271	`/* Loop over all voxels */`
272	`for(int k=0;k<voxel_grid.get_dims()[ 2 ];++k)`	272	`for(int k=0;k<voxel_grid.get_dims()[ 2 ];++k)`
273	`for(int j=0;j<voxel_grid.get_dims()[ 1 ];++j)`	273	`for(int j=0;j<voxel_grid.get_dims()[ 1 ];++j)`
274	`for(int i=0;i<voxel_grid.get_dims()[ 0 ];++i)`	274	`for(int i=0;i<voxel_grid.get_dims()[ 0 ];++i)`
275	`{`	275	`{`
276	`Vec3i pi0(i,j,k);`	276	`Vec3i pi0(i,j,k);`
277	`float val0 = voxel_grid[pi0];`	277	`float val0 = voxel_grid[pi0];`
278		278
279	`/* If the voxel value is above the isovalue */`	279	`/* If the voxel value is above the isovalue */`
280	`if(val0>iso)`	280	`if(val0>iso)`
281	`{`	281	`{`
282	`bool interface_cell = false;`	282	`bool interface_cell = false;`
283	`vector<VertexIter> vertices;`	283	`vector<VertexIter> vertices;`
284	`Cube* cube;`	284	`Cube* cube;`
285		285
286	`// For each neighbouring voxel.`	286	`// For each neighbouring voxel.`
287	`for(int l=0;l<6;++l)`	287	`for(int l=0;l<6;++l)`
288	`{`	288	`{`
289	`// First make sure the voxel is inside the voxel grid.`	289	`// First make sure the voxel is inside the voxel grid.`
290	`Vec3i pi1(pi0);`	290	`Vec3i pi1(pi0);`
291	`pi1 += N6i[l];`	291	`pi1 += N6i[l];`
292	`bool indom = voxel_grid.in_domain(pi1);`	292	`bool indom = voxel_grid.in_domain(pi1);`
293		293
294	`/* If the neighbour is not in the grid or it`	294	`/* If the neighbour is not in the grid or it`
295	`is on the other side of the isovalue */`	295	`is on the other side of the isovalue */`
296	`if(indom && voxel_grid[pi1]<=iso)`	296	`if(indom && voxel_grid[pi1]<=iso)`
297	`{`	297	`{`
298	`/* Store the cube in a hash table. This code is`	298	`/* Store the cube in a hash table. This code is`
299	`executed the first time we create a face`	299	`executed the first time we create a face`
300	`for the cube. In other words only cubes that`	300	`for the cube. In other words only cubes that`
301	`actually have faces are stored in the hash */`	301	`actually have faces are stored in the hash */`
302	`if(!interface_cell)`	302	`if(!interface_cell)`
303	`{`	303	`{`
304	`interface_cell = true;`	304	`interface_cell = true;`
305	`Cube** c_ptr;`	305	`Cube** c_ptr;`
306	`cube_hash.create_entry(HashKey3usi(pi0), c_ptr);`	306	`cube_hash.create_entry(HashKey3usi(pi0), c_ptr);`
307	`Cube c(pi0);`	307	`Cube c(pi0);`
308	`cube_table.push_back(c);`	308	`cube_table.push_back(c);`
309	`cube = (*c_ptr) = &cube_table.back();`	309	`cube = (*c_ptr) = &cube_table.back();`
310	`}`	310	`}`
311	`// Now add the face to the cube.`	311	`// Now add the face to the cube.`
312	`int face_idx = dir_to_face(N6i[l]);`	312	`int face_idx = dir_to_face(N6i[l]);`
313	`CubeFace& face = cube->faces[face_idx];`	313	`CubeFace& face = cube->faces[face_idx];`
314	`face.used = true;`	314	`face.used = true;`
315		315
316	`}`	316	`}`
317	`}`	317	`}`
318	`}`	318	`}`
319	`}`	319	`}`
320	`}`	320	`}`
321		321
322		322
323	`template<class T>`	323	`template<class T>`
324	`void CubeGrid<T>::cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m,`	324	`void CubeGrid<T>::cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m,`
325	`vector<Vec3f>& face_positions)`	325	`vector<Vec3f>& face_positions)`
326	`{`	326	`{`
327	`face_positions.clear();`	327	`face_positions.clear();`
328	`Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&voxel_grid);`	328	`Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&voxel_grid);`
329	`typename list<Cube>::iterator cube_iter = cube_table.begin();`	329	`typename list<Cube>::iterator cube_iter = cube_table.begin();`
330		330
331	`/* For each cube in the table of cubes */`	331	`/* For each cube in the table of cubes */`
332	`cube_iter = cube_table.begin();`	332	`cube_iter = cube_table.begin();`
333	`for(;cube_iter != cube_table.end(); ++cube_iter)`	333	`for(;cube_iter != cube_table.end(); ++cube_iter)`
334	`{`	334	`{`
335	`/* For each of the six faces, if the face is used */`	335	`/* For each of the six faces, if the face is used */`
336	`Cube& cube = *cube_iter;`	336	`Cube& cube = *cube_iter;`
337	`for(int i=0;i<6;++i)`	337	`for(int i=0;i<6;++i)`
338	`if(cube.faces[i].used)`	338	`if(cube.faces[i].used)`
339	`{`	339	`{`
340	`CubeFace& face = cube.faces[i]; // The face`	340	`CubeFace& face = cube.faces[i]; // The face`
341	`FaceIter hmesh_face = m.create_face();`	341	`FaceIter hmesh_face = m.create_face();`
342	`hmesh_face->touched = face_positions.size();`	342	`hmesh_face->touched = face_positions.size();`
343	`Vec3i p0(cube.pos);`	343	`Vec3i p0(cube.pos);`
344	`Vec3i p1 = p0 + N6i[i];`	344	`Vec3i p1 = p0 + N6i[i];`
345	`float v0 = voxel_grid[p0];`	345	`float v0 = voxel_grid[p0];`
346	`float v1 = voxel_grid[p1];`	346	`float v1 = voxel_grid[p1];`
347	`float t = (iso-v0)/(v1-v0);`	347	`float t = (iso-v0)/(v1-v0);`
348	`face_positions.push_back(Vec3f(p0)(1.0f-t)+Vec3f(p1)t);`	348	`face_positions.push_back(Vec3f(p0)(1.0f-t)+Vec3f(p1)t);`
349		349
350	`// Let e be one of the four edges of the face`	350	`// Let e be one of the four edges of the face`
351	`for(int e=0;e<4;++e)`	351	`for(int e=0;e<4;++e)`
352	`{`	352	`{`
353	`HalfEdgeIter he = face.get_he(m,e);`	353	`HalfEdgeIter he = face.get_he(m,e);`
354	`link(he,face.get_he(m,(e+1)%4));`	354	`link(he,face.get_he(m,(e+1)%4));`
355	`link(face.get_he(m,(e+3)%4),he);`	355	`link(face.get_he(m,(e+3)%4),he);`
356	`he->face = hmesh_face;`	356	`he->face = hmesh_face;`
357		357
358	`HalfEdgeIter he_n = edge_to_glue(voxel_grid, cube, i, e);`	358	`HalfEdgeIter he_n = edge_to_glue(voxel_grid, cube, i, e);`
359	`if(he_n != NULL_HALFEDGE_ITER)`	359	`if(he_n != NULL_HALFEDGE_ITER)`
360	`glue(he, he_n);`	360	`glue(he, he_n);`
361	`}`	361	`}`
362	`hmesh_face->last = face.get_he(m,0);`	362	`hmesh_face->last = face.get_he(m,0);`
363	`}`	363	`}`
364		364
365	`}`	365	`}`
366		366
367	`for(HalfEdgeIter h = m.halfedges_begin(); h != m.halfedges_end(); ++h)`	367	`for(HalfEdgeIter h = m.halfedges_begin(); h != m.halfedges_end(); ++h)`
368	`{`	368	`{`
369	`if (h->opp == NULL_HALFEDGE_ITER)`	369	`if (h->opp == NULL_HALFEDGE_ITER)`
370	`{`	370	`{`
371	`HalfEdgeIter ho = m.create_halfedge();`	371	`HalfEdgeIter ho = m.create_halfedge();`
372	`glue(ho,h);`	372	`glue(ho,h);`
373	`}`	373	`}`
374	`}`	374	`}`
375		375
376	`cube_iter = cube_table.begin();`	376	`cube_iter = cube_table.begin();`
377	`for(;cube_iter != cube_table.end(); ++cube_iter)`	377	`for(;cube_iter != cube_table.end(); ++cube_iter)`
378	`{`	378	`{`
379	`//For each of the six faces, if the face is used`	379	`//For each of the six faces, if the face is used`
380	`Cube& cube = *cube_iter;`	380	`Cube& cube = *cube_iter;`
381	`for(int i=0;i<6;++i)`	381	`for(int i=0;i<6;++i)`
382	`if(cube.faces[i].used)`	382	`if(cube.faces[i].used)`
383	`{`	383	`{`
384	`CubeFace& face = cube.faces[i]; // The face`	384	`CubeFace& face = cube.faces[i]; // The face`
385	`// Let e be one of the four edges of the face`	385	`// Let e be one of the four edges of the face`
386	`for(int e=0;e<4;++e)`	386	`for(int e=0;e<4;++e)`
387	`if(face.he[e]->vert == NULL_VERTEX_ITER)`	387	`if(face.he[e]->vert == NULL_VERTEX_ITER)`
388	`{`	388	`{`
389	`HalfEdgeIter last = face.he[e];`	389	`HalfEdgeIter last = face.he[e];`
390	`HalfEdgeIter he = last;`	390	`HalfEdgeIter he = last;`
391	`while(he->opp->face != NULL_FACE_ITER)`	391	`while(he->opp->face != NULL_FACE_ITER)`
392	`{`	392	`{`
393	`he = he->opp->prev;`	393	`he = he->opp->prev;`
394	`if(he == last)`	394	`if(he == last)`
395	`break;`	395	`break;`
396	`}`	396	`}`
397	`last = he;`	397	`last = he;`
398		398
399	`Vec3f pos = cube_corner(cube.pos, i, e);`	399	`Vec3f pos = cube_corner(cube.pos, i, e);`
400	`VertexIter vert = m.create_vertex(pos);`	400	`VertexIter vert = m.create_vertex(pos);`
401		401
402	`do`	402	`do`
403	`{`	403	`{`
404	`he->vert = vert;`	404	`he->vert = vert;`
405	`he = he->next->opp;`	405	`he = he->next->opp;`
406	`}`	406	`}`
407	`while(he->face != NULL_FACE_ITER && he != last);`	407	`while(he->face != NULL_FACE_ITER && he != last);`
408		408
409	`if(he->face == NULL_FACE_ITER)`	409	`if(he->face == NULL_FACE_ITER)`
410	`{`	410	`{`
411	`he->vert = vert;`	411	`he->vert = vert;`
412	`link(he,last->opp);`	412	`link(he,last->opp);`
413	`}`	413	`}`
414		414
415	`vert->he = last->opp;`	415	`vert->he = last->opp;`
416	`}`	416	`}`
417	`}`	417	`}`
418	`}`	418	`}`
419	`}`	419	`}`
420	`}`	420	`}`
421		421
422	`/* A fairly complex algorithm which triangulates faces. It triangulates`	422	`/* A fairly complex algorithm which triangulates faces. It triangulates`
423	`by iteratively splitting faces. It always splits a face by picking the`	423	`by iteratively splitting faces. It always splits a face by picking the`
424	`edge whose midpoint is closest to the isovalue. */`	424	`edge whose midpoint is closest to the isovalue. */`
425		425
426	`template<class T>`	426	`template<class T>`
427	`void triangulate(const RGrid<T>& voxel_grid, Manifold& m, float iso)`	427	`void triangulate(const RGrid<T>& voxel_grid, Manifold& m, float iso)`
428	`{`	428	`{`
429	`#ifndef NDEBUG`	429	`#ifndef NDEBUG`
430	`cout << "Initial mesh valid : " << m.is_valid() << endl;`	430	`cout << "Initial mesh valid : " << m.is_valid() << endl;`
431	`#endif`	431	`#endif`
432	`Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(&voxel_grid);`	432	`Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(&voxel_grid);`
433	`int work;`	433	`int work;`
434	`do`	434	`do`
435	`{`	435	`{`
436	`// For every face.`	436	`// For every face.`
437	`work = 0;`	437	`work = 0;`
438	`for(FaceIter fi =m.faces_begin(); fi != m.faces_end(); ++fi)`	438	`for(FaceIter fi =m.faces_begin(); fi != m.faces_end(); ++fi)`
439	`{`	439	`{`
440	`// Create a vector of vertices.`	440	`// Create a vector of vertices.`
441	`vector<VertexIter> verts;`	441	`vector<VertexIter> verts;`
442	`for(FaceCirculator fc(fi); !fc.end(); ++fc)`	442	`for(FaceCirculator fc(fi); !fc.end(); ++fc)`
443	`verts.push_back(fc.get_vertex());`	443	`verts.push_back(fc.get_vertex());`
444		444
445	`// If there are just three we are done.`	445	`// If there are just three we are done.`
446	`if(verts.size() == 3) continue;`	446	`if(verts.size() == 3) continue;`
447		447
448	`// Find vertex pairs that may be connected.`	448	`// Find vertex pairs that may be connected.`
449	`vector<pair<int,int> > vpairs;`	449	`vector<pair<int,int> > vpairs;`
450	`const int N = verts.size();`	450	`const int N = verts.size();`
451	`for(int i=0;i<N-2;++i)`	451	`for(int i=0;i<N-2;++i)`
452	`for(int j=i+2;j<N; ++j)`	452	`for(int j=i+2;j<N; ++j)`
453	`if(!is_connected(verts[i], verts[j]))`	453	`if(!is_connected(verts[i], verts[j]))`
454	`vpairs.push_back(pair<int,int>(i,j));`	454	`vpairs.push_back(pair<int,int>(i,j));`
455		455
456	`if(vpairs.empty())`	456	`if(vpairs.empty())`
457	`{`	457	`{`
458	`cout << "Warning: could not triangulate a face."`	458	`cout << "Warning: could not triangulate a face."`
459	`<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;`	459	`<< "Probably a vertex appears more than one time in other vertex's one-ring" << endl;`
460	`continue;`	460	`continue;`
461	`}`	461	`}`
462		462
463	`/* For all vertex pairs, find the absolute value of iso-x`	463	`/* For all vertex pairs, find the absolute value of iso-x`
464	`where x is the interpolated value of the volume at the`	464	`where x is the interpolated value of the volume at the`
465	`midpoint of the line segment connecting the two vertices.*/`	465	`midpoint of the line segment connecting the two vertices.*/`
466		466
467	`float min_val=FLT_MAX;`	467	`float min_val=FLT_MAX;`
468	`int min_k = -1;`	468	`int min_k = -1;`
469	`for(size_t k=0;k<vpairs.size(); ++k)`	469	`for(size_t k=0;k<vpairs.size(); ++k)`
470	`{`	470	`{`
471	`int i = vpairs[k].first;`	471	`int i = vpairs[k].first;`
472	`int j = vpairs[k].second;`	472	`int j = vpairs[k].second;`
473	`Vec3f centre =`	473	`Vec3f centre =`
474	`(verts[i]->pos +`	474	`(verts[i]->pos +`
475	`verts[j]->pos)/2.0f;`	475	`verts[j]->pos)/2.0f;`
476		476
477	`float val;`	477	`float val;`
478	`if(grid_interp.in_domain(centre))`	478	`if(grid_interp.in_domain(centre))`
479	`val = fabs(grid_interp(centre)-iso);`	479	`val = fabs(grid_interp(centre)-iso);`
480	`else`	480	`else`
481	`val = 0.0f;`	481	`val = 0.0f;`
482	`if(val<min_val)`	482	`if(val<min_val)`
483	`{`	483	`{`
484	`min_val = val;`	484	`min_val = val;`
485	`min_k = k;`	485	`min_k = k;`
486	`}`	486	`}`
487	`}`	487	`}`
488	`assert(min_k != -1);`	488	`assert(min_k != -1);`
489		489
490	`{`	490	`{`
491	`// Split faces along edge whose midpoint is closest to isovalue`	491	`// Split faces along edge whose midpoint is closest to isovalue`
492	`int i = vpairs[min_k].first;`	492	`int i = vpairs[min_k].first;`
493	`int j = vpairs[min_k].second;`	493	`int j = vpairs[min_k].second;`
494	`m.split_face(fi, verts[i], verts[j]);`	494	`m.split_face(fi, verts[i], verts[j]);`
495	`}`	495	`}`
496	`++work;`	496	`++work;`
497		497
498	`}`	498	`}`
499	`}`	499	`}`
500	`while(work);`	500	`while(work);`

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/fair_polygonize.cpp @ 182 – Rev 170 → 182