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

Subversion Repositories gelsvn

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