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

Subversion Repositories gelsvn

(root)/trunk/src/HMesh/volume_polygonize.cpp @ 201 – Rev 191 → 192