WebSVN – gelsvn – Diff – /trunk/src/HMesh/volume_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 "volume_polygonize.h"
 
namespace HMesh
{
	using namespace Geometry;
	using namespace Util;
	using namespace CGLA;
	using namespace std;
 
	namespace 
	{
 
		template<class T>
		Vec3f push_vertex(const RGrid<T>& grid, const Vec3f& pos0, float iso)
		{
			Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&grid);
			Vec3f corner_min = pos0 - Vec3f(1);
			Vec3f corner_max = pos0 + Vec3f(1);
			Vec3f pos = pos0;
			for(int i=0;i<1000;++i)
				if(inter.in_domain(pos))
					{
						float val = inter(pos)-iso;
						Vec3f grad = inter.grad(pos);
						float sqlen_grad = sqr_length(grad);
 
						// Compute the projection in the gradient direction
						Vec3f d1 = - val * grad/sqlen_grad;
						Vec3f p_new = pos + 0.5 * d1; // "0.5 *" to add damping
 
						// Don't move outside the box.
						if(!(corner_max.all_g(p_new) && 
								 corner_min.all_l(p_new)))
							return pos;
								
						// Move so that the gradient points in the direction
						// of the original position
						Vec3f d2 = (pos0-pos) - 
							grad*dot((pos0-pos),grad)/sqlen_grad;
									
						// Compute new position and distance up to original position.
						p_new += 0.5 * d2; // "0.5 *" to add damping
									
						// Don't move outside the box.
						if(corner_max.all_g(p_new) && 
							 corner_min.all_l(p_new))
							pos = p_new;
						else
							return pos;
 
						// Did we converge?
						if(sqrt(sqr_length(d1)+sqr_length(d2))<.001)
							return pos;
					}
			return pos;
		}
								
		void compute_dual(Manifold& m, Manifold& m2, const vector<Vec3f>& face_positions)
		{
			// make sure every face knows its number
			m.enumerate_faces();
 
			int i=0;
			vector<VertexIter> vertices(m.no_faces());
			for(FaceIter f=m.faces_begin(); f!=m.faces_end(); ++f, ++i)
					vertices[i] = m2.create_vertex(face_positions[f->touched]);
 
 
			vector<HalfEdgeIter> hedges(m.no_halfedges(), NULL_HALFEDGE_ITER);
			int he_num=0;
			for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v)
				if(!is_boundary(v))
				{
						FaceIter f = m2.create_face();
						vector<int> indices;
						for(VertexCirculator vc(v); !vc.end(); ++vc)
						{
								TouchType t = vc.get_halfedge()->touched;
								int h_idx = he_num++;
								vc.get_halfedge()->touched = h_idx;
								HalfEdgeIter h = m2.create_halfedge();
								h->face = f;
								h->touched = t;
								
								int f_idx = vc.get_face()->touched;
								h->vert = vertices[f_idx];
								hedges[h_idx] = h;
 
								indices.push_back(h_idx);
						}
						f->last = hedges[indices[0]];
						int N=indices.size();
						for(int j=0;j<N;++j)
						{
								HalfEdgeIter h1 = hedges[indices[j]];
								HalfEdgeIter h2 = hedges[indices[(j+1)%N]];
								
								link(h2,h1);
								h2->vert->he = h1;
						}
				}
			for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v)
				if(!is_boundary(v))
						for(VertexCirculator vc(v); !vc.end(); ++vc)
								if(!is_boundary(vc.get_vertex()))
								{
										int idx0 = vc.get_halfedge()->touched;
										int idx1 = vc.get_opp_halfedge()->touched;
										glue(hedges[idx0], hedges[idx1]);
								}
 
			for(HalfEdgeIter h = m2.halfedges_begin(); h != m2.halfedges_end(); ++h)
					if(h->opp == NULL_HALFEDGE_ITER)
							glue(h, m2.create_halfedge());
			for(VertexIter v = m2.vertices_begin(); v != m2.vertices_end(); ++v)
					check_boundary_consistency(v);
		}
 
		/*
			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(): 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, 
													bool push_vertices,
													bool compute_face_positions,
													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);
			
			HalfEdgeIter he = cube_n->faces[i_n].he[e_n];
			if(cube_n == &cube)
					he->touched = 1;
			else
					he->touched = 0;
			
			return he;
		}
				
		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, 
																		 bool push_vertices,
																		 bool compute_face_positions,
																		 vector<Vec3f>& face_positions)
		{
			face_positions.clear();
			
			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();
 
								if(compute_face_positions)
								{
										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)
								{
										face.he[e] = m.create_halfedge();
										face.he[e]->face = hmesh_face;
								}
								hmesh_face->last = face.he[0];
							}
 
				}
 
			/* 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
								// Let e be one of the four edges of the face
								for(int e=0;e<4;++e)
									{
											HalfEdgeIter he = face.he[e];
											link(he,face.he[(e+1)%4]);
											link(face.he[(e+3)%4],he);
											if(he->opp == NULL_HALFEDGE_ITER)
											{
													HalfEdgeIter he_n = 
															edge_to_glue(voxel_grid, cube, i, e);
													if(he_n != NULL_HALFEDGE_ITER)
													{
															glue(he, he_n);
															he->touched = he_n->touched;
													}
											}
									}
							}
 
				}
 
			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);
											if(push_vertices)
													pos = push_vertex(voxel_grid, pos, iso);
											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;
										}
							}
				}
		}
	}
	
 
	
	template<typename T>
	void mc_polygonize(const RGrid<T>& voxel_grid, 
											 Manifold& mani_out, 
											 float iso)
	{
			Manifold mani;
			CubeGrid<T> cube_grid(voxel_grid, iso);
			vector<Vec3f> face_positions;
			cube_grid.cuberille_mesh(voxel_grid, mani, false, true, face_positions);
			compute_dual(mani, mani_out, 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, push, false, face_positions);
	}
 
		template void mc_polygonize<short>(const RGrid<short>&,
																			 Manifold&, float);
		template void mc_polygonize<unsigned short>(const RGrid<unsigned short>&,
																							 Manifold&, float);
		template void mc_polygonize<unsigned char>(const RGrid<unsigned char>&,
																							 Manifold&, float);
		template void mc_polygonize<float>(const RGrid<float>&, Manifold&, float);
 
		template void cuberille_polygonize<short>(const RGrid<short>&,
																							Manifold&, float, bool);
 
		template void cuberille_polygonize<unsigned short>(const RGrid<unsigned short>&,
																											 Manifold&, float, bool);
 
		template void cuberille_polygonize<unsigned char>(const RGrid<unsigned char>&,
																											Manifold&, float, bool);
 
		template void cuberille_polygonize<float>(const RGrid<float>&,
																							Manifold&, float, bool);
 
}
	
 
 
	
	
 

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

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/volume_polygonize.cpp @ 512 – Rev 192 → 215