Subversion Repositories gelsvn

Rev

Rev 169 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 169 Rev 170
Line 7... Line 7...
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 
 
 13 
 
13 
 
 14 
 
14 
#include "Geometry/RGrid.h"
 15 
#include "Geometry/RGrid.h"
15 
#include "Geometry/TrilinFilter.h"
 16 
#include "Geometry/TrilinFilter.h"
16 
#include "Geometry/GradientFilter.h"
 17 
#include "Geometry/GradientFilter.h"
Line 20... Line 21...
20 
 
 21 
 
21 
#include "fair_polygonize.h"
 22 
#include "fair_polygonize.h"
22 
 
 23 
 
23 
namespace HMesh
 24 
namespace HMesh
24 
{
 25 
{
25 
		using namespace Geometry;
 26 
	using namespace Geometry;
26 
		using namespace Util;
 27 
	using namespace Util;
27 
		using namespace CGLA;
 28 
	using namespace CGLA;
28 
		using namespace std;
 29 
	using namespace std;
29 
 
 - 
 
30 
		namespace
- 
 
31 
		{
 - 
 
32 
				/*
 - 
 
33 
					Vectors along edges. For each of six cube faces there are four such vectors
 - 
 
34 
					since the faces of a cube are quads.
 - 
 
35 
				*/
 - 
 
36 
				const Vec3i edges[6][4] = {
 - 
 
37 
						{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},
 - 
 
38 
						{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},
 - 
 
39 
						{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},
 - 
 
40 
						{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},
 - 
 
41 
						{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},
 - 
 
42 
						{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};
 - 
 
43 
 
 30 
 
44 
				/* Convert a direction vector to a cube face index. */
 - 
 
45 
				int dir_to_face(const Vec3i& dir)
 31 
	namespace
46 
				{
 32 
	{
47 
						assert(dot(dir,dir)==1);
 - 
 
48 
 
 33 
 
49 
						if(dir[0] != 0)
 - 
 
50 
								return (dir[0]<0) ? 0 : 1;
 34 
		void compute_dual(Manifold& m, const vector<Vec3f>& face_positions)
51 
						if(dir[1] != 0)
 35 
		{
52 
								return (dir[1]<0) ? 2 : 3;
 36 
			// make sure every face knows its number
53 
						if(dir[2] != 0)
 37 
			m.enumerate_faces();
54 
								return (dir[2]<0) ? 4 : 5;
 - 
 
55 
 
 38 
 
- 
 
 39 
			vector<Vec3f> vertices(m.no_faces());
- 
 
 40 
			vector<int> faces;
- 
 
 41 
			vector<int> indices;
- 
 
 42 
 
- 
 
 43 
			// Create new vertices. Each face becomes a vertex whose positio
- 
 
 44 
			// is the centre of the face
- 
 
 45 
			int i=0;
- 
 
 46 
			for(FaceIter f=m.faces_begin(); f!=m.faces_end(); ++f,++i)
- 
 
 47 
				vertices[i] = face_positions[f->touched];
- 
 
 48 
 
- 
 
 49 
			// Create new faces. Each vertex is a new face with N=valency of
- 
 
 50 
			// edges.
- 
 
 51 
			i=0;
- 
 
 52 
			for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v,++i)
- 
 
 53 
				if(!is_boundary(v))
- 
 
 54 
					{
- 
 
 55 
						int no_edges = 0;
- 
 
 56 
						vector<int> index_tmp;
- 
 
 57
- 
 
 58 
						// Circulate around vertex
- 
 
 59 
						VertexCirculator vc(v);
- 
 
 60 
						for(; !vc.end(); ++vc)
- 
 
 61 
							{
- 
 
 62 
								FaceIter f = vc.get_face();
- 
 
 63 
								if(f != NULL_FACE_ITER)
- 
 
 64 
									{
- 
 
 65 
										index_tmp.push_back(f->touched);
- 
 
 66 
										++no_edges;
- 
 
 67 
									}
- 
 
 68 
							}
- 
 
 69 
						// Push vertex indices for this face onto indices vector
- 
 
 70
- 
 
 71 
						// The circulator moves around the face in a clockwise f
- 
 
 72 
						// so we just reverse the ordering.
- 
 
 73 
						indices.insert(indices.end(), index_tmp.rbegin(), index_tmp.rend());
- 
 
 74 
						faces.push_back(no_edges);
- 
 
 75 
					}
- 
 
 76 
 
- 
 
 77 
			// Clear the manifold before new geometry is inserted.
- 
 
 78 
			m.clear();
- 
 
 79 
 
- 
 
 80 
			// And build
- 
 
 81 
			build_manifold(m, vertices.size(), &vertices[0], faces.size(),
- 
 
 82 
										 &faces[0],&indices[0]);
- 
 
 83 
 
- 
 
 84 
		}
- 
 
 85 
 
- 
 
 86 
		/*
- 
 
 87 
			Vectors along edges. For each of six cube faces there are four such vectors
- 
 
 88 
			since the faces of a cube are quads.
- 
 
 89 
		*/
- 
 
 90 
		const Vec3i edges[6][4] = {
- 
 
 91 
			{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},
- 
 
 92 
			{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},
- 
 
 93 
			{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},
- 
 
 94 
			{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},
- 
 
 95 
			{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},
- 
 
 96 
			{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};
- 
 
 97 
 
- 
 
 98 
		/* Convert a direction vector to a cube face index. */
- 
 
 99 
		int dir_to_face(const Vec3i& dir)
- 
 
 100 
		{
- 
 
 101 
			assert(dot(dir,dir)==1);
- 
 
 102 
 
- 
 
 103 
			if(dir[0] != 0)
- 
 
 104 
				return (dir[0]<0) ? 0 : 1;
- 
 
 105 
			if(dir[1] != 0)
- 
 
 106 
				return (dir[1]<0) ? 2 : 3;
- 
 
 107 
			if(dir[2] != 0)
- 
 
 108 
				return (dir[2]<0) ? 4 : 5;
- 
 
 109 
 
56 
						return -1;
 110 
			return -1;
57 
				}
 111 
		}
58
112
59 
				/* Convert a face index and an edge vector to an edge index. */
 113 
		/* Convert a face index and an edge vector to an edge index. */
60 
				int dir_to_edge(int face, const Vec3i& edge)
 114 
		int dir_to_edge(int face, const Vec3i& edge)
61 
				{
 115 
		{
62 
						for(int i=0;i<4;++i)
 116 
			for(int i=0;i<4;++i)
63 
								if(edges[face][i] == edge) return i;
 117 
				if(edges[face][i] == edge) return i;
64 
						assert(0);
 118 
			assert(0);
65 
						return -1;
 119 
			return -1;
- 
 
 120 
		}
- 
 
 121 
 
- 
 
 122 
		Vec3f cube_corner(const Vec3i& pos, int i, int e)
- 
 
 123 
		{
- 
 
 124 
			Vec3i to_face = N6i[i];
- 
 
 125 
			Vec3i along_edge = edges[i][e];
- 
 
 126 
			Vec3i to_edge = cross(along_edge,to_face);
- 
 
 127 
			return Vec3f(pos) + 0.5*Vec3f(to_face+to_edge+along_edge);
66 
				}
 128 
		}
67
129
68 
				/* Cube faces are created whenever a voxel that is inside is adjacent
 130 
		/* Cube faces are created whenever a voxel that is inside is adjacent
69 
					 to a voxel that is outside. Thus the voxel can be seen as a small
 131 
			 to a voxel that is outside. Thus the voxel can be seen as a small
70 
					 box that either contains material or not. A face is created where
 132 
			 box that either contains material or not. A face is created where
71 
					 full and empty boxes meet. A vertex is placed on every box face.
 133 
			 full and empty boxes meet. A vertex is placed on every box face.
72
134
73 
					 A halfedge is created for every edge of the (quadrilateral) face.
 135 
			 A halfedge is created for every edge of the (quadrilateral) face.
74 
					 This halfedge connects the face to its neighbour. */
 136 
			 This halfedge connects the face to its neighbour. */
75 
				struct CubeFace
 137 
		struct CubeFace
76 
				{
 138 
		{
77 
						bool used;
 139 
			bool used;
78 
						VertexIter vert;
 - 
 
79 
						HalfEdgeIter he[4];
 140 
			HalfEdgeIter he[4];
80
141
81 
						CubeFace();
 142 
			CubeFace();
82 
 
 143 
 
83 
						HalfEdgeIter& get_he(Manifold& m, int i)
 144 
			HalfEdgeIter& get_he(Manifold& m, int i)
84 
								{
 145 
			{
85 
										if(he[i] == NULL_HALFEDGE_ITER)
 146 
				if(he[i] == NULL_HALFEDGE_ITER)
86 
												he[i] = m.create_halfedge();
 147 
					he[i] = m.create_halfedge();
87 
										return he[i];
 148 
				return he[i];
88 
								}
 149 
			}
89 
				};
 150 
		};
90 
 
 151 
 
91 
				CubeFace::CubeFace(): used(false), vert(NULL_VERTEX_ITER)
 152 
		CubeFace::CubeFace(): used(false)
92 
				{
 153 
		{
93 
						for(int i=0;i<4;++i)
 154 
			for(int i=0;i<4;++i)
94 
								he[i] = NULL_HALFEDGE_ITER;
 155 
				he[i] = NULL_HALFEDGE_ITER;
95 
				}
 156 
		}
96
157
97 
 
 158 
 
98 
				/* A cube consists of six faces and a position. The cube centre is
 159 
		/* A cube consists of six faces and a position. The cube centre is
99 
					 an integer (voxel grid) position. */
 160 
			 an integer (voxel grid) position. */
100 
				struct Cube
 161 
		struct Cube
101 
				{
 162 
		{
102 
						Vec3i pos;
 163 
			Vec3i pos;
103 
						CubeFace faces[6];
 164 
			CubeFace faces[6];
104
165
105 
						Cube(const Vec3i& _pos): pos(_pos) {}
 166 
			Cube(const Vec3i& _pos): pos(_pos) {}
106 
				};
 167 
		};
107 
 
 168 
 
108 
 
 169 
 
109 
				/* The mesher class does the actual work. It is not pretty and it
 170 
		/* The mesher class does the actual work. It is not pretty and it
110 
					 is not seen by the user. */
 171 
			 is not seen by the user. */
111 
				template<class T>
 172 
		template<class T>
112 
				class CubeGrid
 173 
		class CubeGrid
113 
				{
 174 
		{
114 
						float iso;
 175 
			float iso;
115 
						HashTable<HashKey3usi,Cube*> cube_hash;  // Hashtable of cube pointers
 176 
			HashTable<HashKey3usi,Cube*> cube_hash;  // Hashtable of cube pointers
116 
						list<Cube> cube_table;                   // Linked list containing
 177 
			list<Cube> cube_table;                   // Linked list containing
117 
                                          					 // the actual cubes.
 178 
			// the actual cubes.
118 
 
 179 
 
119 
				public:
 180 
		public:
120
181
121 
						CubeGrid(const RGrid<T>& _voxel_grid, // input voxel grid
 182 
			CubeGrid(const RGrid<T>& _voxel_grid, // input voxel grid
122 
										 float iso,                   // iso value
 183 
							 float iso);                   // iso value
123 
										 bool inv_cont);              // invert contour?
 - 
 
124 
 
 - 
 
125 
						void find_neighbour_face_disconnect(Cube& cube, int face_idx, int edge_idx,
 - 
 
126 
																								Cube*& cube_n, int& i_n, int& e_n);
 - 
 
127 
 
 - 
 
128 
						void find_neighbour_face_connect(Cube& cube, int face_idx, int edge_idx,
 - 
 
129 
																						 Cube*& cube_n, int& i_n, int& e_n);
 - 
 
130
- 
 
131 
						void dual_cuberille_mesh(const RGrid<T>&, Manifold& m, bool connect=false);
 - 
 
132 
						void cuberille_mesh(const RGrid<T>&, Manifold& m, bool connect=false);
 - 
 
133 
				};
 - 
 
134 
 
 - 
 
135 
				template<class T>
 - 
 
136 
				void CubeGrid<T>::find_neighbour_face_disconnect(Cube& cube, int face_idx, int edge_idx,
 - 
 
137 
																												 Cube*& cube_n,  // Neighbouring cube
 - 
 
138 
																												 int& i_n,  // direction from center to neighbour face
 - 
 
139 
																												 int& e_n) // Edge vector for neighbour face
 - 
 
140 
				{
 - 
 
141 
						/*
 - 
 
142 
							Algorithm: Pick faces to glue
 - 
 
143 
 
 184 
 
144 
							Algorithm which connects cube faces. If both A and B
 - 
 
145 
							are above the isovalue and 0 and C are below, A and B
 - 
 
146 
							will be seperated, i.e. the 0A and AC faces will be
- 
 
147 
							connected. There is a different out-commented algorithm
 - 
 
148 
							at the end of this file which connects AO and 0B
 - 
 
149 
							instead.
 - 
 
150
- 
 
151 
							Figure.
 - 
 
152 
 
 185 
 
153 
							0|B
 - 
 
154 
							-+-
 - 
 
155 
							A|C
 - 
 
156
- 
 
157 
							Explanation:
 - 
 
158
- 
 
159 
							Let the current cube be A. Since we are here, the voxel
 - 
 
160 
							of A is above the isovalue, and the voxel of O is below
 186 
			HalfEdgeIter edge_to_glue(const RGrid<T>& voxel_grid,
161 
							since there is a face AO.
 - 
 
162 
 
 - 
 
163 
							We first check if there is a face between A and C.
- 
 
164 
							if that is the case, the AO face is connected to the AC
187 
																Cube& cube, int face_idx, int edge_idx);
165 
							face.
- 
 
166 
 
 - 
 
167 
							Otherwise, we check if B is also inside the isocontour.
 - 
 
168 
							If that is true there is an OB face, and we connect the
 - 
 
169 
							A0 and the 0B faces.
 - 
 
170
- 
 
171 
							Failing that, there must be a CB face, and we connect
 - 
 
172 
							with that.
 - 
 
173 
						*/
 - 
 
174 
						Vec3i pos = cube.pos;
 - 
 
175 
						Vec3i dir = N6i[face_idx];
 - 
 
176 
 
 - 
 
177 
 
 - 
 
178 
						Vec3i dir_n;   // direction from center to neighbour face
 - 
 
179 
						Vec3i edge_n; // Edge vector for neighbour face
 - 
 
180 
 
 - 
 
181 
						Vec3i edge = edges[face_idx][edge_idx];	       // Vector _along_ edge
 - 
 
182 
						Vec3i to_edge = cross(edge,dir); // vector to edge
 - 
 
183 
 
 - 
 
184
- 
 
185 
						// Begin "Pick faces to glue"
 - 
 
186 
						Cube* dummy;
 - 
 
187 
						if(cube.faces[dir_to_face(to_edge)].used)
 - 
 
188 
						{
 - 
 
189 
								dir_n = to_edge;
 - 
 
190 
								edge_n = - edge;
- 
 
191 
								cube_n = &cube;
- 
 
192 
						}
 - 
 
193 
						else if(!cube_hash.find_entry(HashKey3usi(pos+dir+to_edge),
 - 
 
194 
																					dummy))
 - 
 
195 
						{
 188
196 
								assert(sqr_length(to_edge)==1);
 - 
 
197 
								dir_n = dir;
 - 
 
198 
								edge_n = - edge;
 - 
 
199 
								cube_hash.find_entry(HashKey3usi(pos + to_edge),
189 
			void cuberille_mesh(const RGrid<T>&, Manifold& m, vector<Vec3f>& face_positions);
200 
																		 cube_n);
 - 
 
201 
						}
 - 
 
202 
						else
- 
 
203 
						{
 - 
 
204 
								dir_n = - to_edge;
 - 
 
205 
								edge_n = - edge;
 - 
 
206 
								cube_n = dummy;
 - 
 
207 
						}
 - 
 
208 
 
 - 
 
209 
						i_n = dir_to_face(dir_n);
 - 
 
210 
						e_n = dir_to_edge(i_n, edge_n);
 - 
 
211 
				}
 190 
		};
212 
 
 191 
 
213 
				template<class T>
 192 
		template<class T>
214 
				void CubeGrid<T>::find_neighbour_face_connect(Cube& cube, int face_idx, int edge_idx,
 193 
		HalfEdgeIter CubeGrid<T>::edge_to_glue(const RGrid<T>& voxel_grid,
215 
																											Cube*& cube_n,  // Neighbouring cube
 194 
																					 Cube& cube, int face_idx, int edge_idx)
216 
																											int& i_n,  // direction from center to neighbour face
 - 
 
217 
																											int& e_n) // Edge vector for neighbour face
 - 
 
218 
				{
 195 
		{
219 
						Vec3i pos = cube.pos;
 196 
			/*  Figure.
220 
						Vec3i dir = N6i[face_idx];
 - 
 
221 
 
 197 
 
- 
 
 198 
			0|B
- 
 
 199 
			-+-
- 
 
 200 
			A|C
222 
						Vec3i dir_n;   // direction from center to neighbour face
 201
223 
						Vec3i edge_n; // Edge vector for neighbour face
 202 
			Explanation: ... pending
- 
 
 203 
			*/
224 
 
 204 
 
- 
 
 205 
			Vec3i dir = N6i[face_idx];
225 
						Vec3i edge = edges[face_idx][edge_idx];	       // Vector _along_ edge
 206 
			Vec3i along_edge = edges[face_idx][edge_idx];	 // Vector _along_ edge
226 
						Vec3i to_edge = cross(edge,dir); // vector to edge
 207 
			Vec3i to_edge = cross(along_edge,dir); // vector to edge
- 
 
 208 
 
- 
 
 209 
			Vec3i a_pos = cube.pos;
- 
 
 210 
			Vec3i o_pos = a_pos + dir;
- 
 
 211 
			Vec3i b_pos = o_pos + to_edge;
- 
 
 212 
			Vec3i c_pos = a_pos + to_edge;
- 
 
 213
- 
 
 214 
			if(!voxel_grid.in_domain(b_pos))
- 
 
 215 
				return NULL_HALFEDGE_ITER;
227 
 
 216 
 
- 
 
 217 
			Cube* cube_n=0;
- 
 
 218 
			Vec3i dir_n;   // direction from center to neighbour face
- 
 
 219 
			Vec3i along_edge_n; // Edge vector for neighbour face
228 
						if(cube_hash.find_entry(HashKey3usi(pos + dir + to_edge),
 220 
			if(!cube_hash.find_entry(HashKey3usi(b_pos),	cube_n))
- 
 
 221 
				{
229 
																		cube_n))
 222 
					if(cube.faces[dir_to_face(to_edge)].used)
230 
						{
 223 
						{
231 
								assert(sqr_length(dir + to_edge)==2);
 224 
							dir_n = to_edge;
232 
																dir_n = - to_edge;
 225 
							along_edge_n = - along_edge;
233 
																edge_n = - edge;
 226 
							cube_n = &cube;
234 
						}
 227 
						}
- 
 
 228 
					else
- 
 
 229 
						{
- 
 
 230 
							dir_n = dir;
- 
 
 231 
							along_edge_n = - along_edge;
235 
						else if(cube_hash.find_entry(HashKey3usi(pos + to_edge),
232 
							cube_hash.find_entry(HashKey3usi(c_pos), cube_n);
- 
 
 233 
							assert(cube_n->faces[dir_to_face(dir_n)].used);
236 
																				 cube_n))
 234 
						}
- 
 
 235 
				}
- 
 
 236 
			else
- 
 
 237 
				{
- 
 
 238 
					float a_val = voxel_grid[a_pos];
- 
 
 239 
					float b_val = voxel_grid[b_pos];
- 
 
 240 
					float c_val = voxel_grid[c_pos];
- 
 
 241 
					float o_val = voxel_grid[o_pos];
- 
 
 242 
					bool connect = (a_val*b_val-c_val*o_val)/(a_val+b_val-c_val-o_val) >iso;
- 
 
 243
- 
 
 244 
					if(connect || !cube.faces[dir_to_face(to_edge)].used)
237 
						{
 245 
						{
238 
								assert(sqr_length(to_edge)==1);
 246 
							dir_n = - to_edge;
239 
								dir_n = dir;
 247 
							along_edge_n = - along_edge;
240 
								edge_n = - edge;
 248 
							assert(cube_n->faces[dir_to_face(dir_n)].used);
241 
						}
 249 
						}
242 
						else
 250 
					else
243 
						{
 251 
						{
244 
								dir_n = to_edge;
 252 
							dir_n = to_edge;
245 
								edge_n = - edge;
 253 
							along_edge_n = - along_edge;
246 
								cube_n = &cube;
 254 
							cube_n = &cube;
247 
						}
 255 
						}
248
- 
 
249 
						i_n = dir_to_face(dir_n);
 - 
 
250 
						e_n = dir_to_edge(i_n, edge_n);
 - 
 
251 
				}
 256 
				}
- 
 
 257 
			int i_n = dir_to_face(dir_n);
- 
 
 258 
			int e_n = dir_to_edge(i_n, along_edge_n);
- 
 
 259 
			return cube_n->faces[i_n].he[e_n];
252 
 
 260 
		}
253
261
254 
				template<class T>
 262 
		template<class T>
255 
				inline bool comp(bool d, T a, T b)
 263 
		inline bool comp(bool d, T a, T b)
256 
				{
 264 
		{
257 
						return d? (a<b) : (b<a);
265 
			return d? (a<b) : (b<a);
258 
				}
 266 
		}
259
267
260 
				template<class T>
 268 
		template<class T>
261 
				CubeGrid<T>::CubeGrid(const RGrid<T>& voxel_grid, float _iso, bool inv_cont): iso(_iso)
 269 
		CubeGrid<T>::CubeGrid(const RGrid<T>& voxel_grid, float _iso): iso(_iso)
262 
				{
 270 
		{
263 
						/* Loop over all voxels */
 271 
			/* Loop over all voxels */
264 
						for(int k=0;k<voxel_grid.get_dims()[ 2 ];++k)
 272 
			for(int k=0;k<voxel_grid.get_dims()[ 2 ];++k)
265 
								for(int j=0;j<voxel_grid.get_dims()[ 1 ];++j)
 273 
				for(int j=0;j<voxel_grid.get_dims()[ 1 ];++j)
266 
										for(int i=0;i<voxel_grid.get_dims()[ 0 ];++i)
 274 
					for(int i=0;i<voxel_grid.get_dims()[ 0 ];++i)
267 
										{
 275 
						{
268 
												Vec3i pi0(i,j,k);
 276 
							Vec3i pi0(i,j,k);
269 
												float val0 = voxel_grid[pi0];
 277 
							float val0 = voxel_grid[pi0];
270 
 
 278 
 
271 
												/* If the voxel value is above the isovale (or below if we
 279 
							/* If the voxel value is above the isovalue */
272 
													 are finding the inverse contour) */
 - 
 
273 
												if(comp(inv_cont,val0,iso))
 280 
							if(val0>iso)
274 
												{
 281 
								{
275 
														bool interface_cell = false;
 282 
									bool interface_cell = false;
276 
														vector<VertexIter> vertices;
 283 
									vector<VertexIter> vertices;
277 
														Cube* cube;
 284 
									Cube* cube;
278 
 
 285 
 
279 
														// For each neighbouring voxel.
 286 
									// For each neighbouring voxel.
280 
														for(int l=0;l<6;++l)
 287 
									for(int l=0;l<6;++l)
281 
														{
 288 
										{
282 
																// First make sure the voxel is inside the voxel grid.
 289 
											// First make sure the voxel is inside the voxel grid.
283 
																Vec3i pi1(pi0);
 290 
											Vec3i pi1(pi0);
284 
																pi1 += N6i[l];
 291 
											pi1 += N6i[l];
285 
																bool indom = voxel_grid.in_domain(pi1);
 292 
											bool indom = voxel_grid.in_domain(pi1);
286
293
287 
																/* If the neighbour is **not** in the grid or it
294 
											/* If the neighbour is **not** in the grid or it
288 
																	 is on the other side of the isovalue */
 295 
												 is on the other side of the isovalue */
289 
																if(!indom ||
- 
 
290 
																	 !comp(inv_cont,float(voxel_grid[pi1]),iso))
 296 
											if(indom && voxel_grid[pi1]<=iso)
291 
																{
 297 
												{
292 
																		/* Store the cube in a hash table. This code is
 298 
													/* Store the cube in a hash table. This code is
293 
																			 executed the first time we create a face
 299 
														 executed the first time we create a face
294 
																			 for the cube. In other words only cubes that
300 
														 for the cube. In other words only cubes that
295 
																			 actually have faces are stored in the hash */
 301 
														 actually have faces are stored in the hash */
296 
																		if(!interface_cell)
 302 
													if(!interface_cell)
297 
																		{
 303 
														{
298 
																				interface_cell = true;
 304 
															interface_cell = true;
299 
																				Cube** c_ptr;
 305 
															Cube** c_ptr;
300 
																				cube_hash.create_entry(HashKey3usi(pi0), c_ptr);
 306 
															cube_hash.create_entry(HashKey3usi(pi0), c_ptr);
301 
																				Cube c(pi0);
 307 
															Cube c(pi0);
302 
																				cube_table.push_back(c);
 308 
															cube_table.push_back(c);
303 
																				cube = (*c_ptr) = &cube_table.back();
 309 
															cube = (*c_ptr) = &cube_table.back();
304 
																		}
 - 
 
305 
																		// Now add the face to the cube.
 - 
 
306 
																		int face_idx  = dir_to_face(N6i[l]);
 - 
 
307 
																		CubeFace& face = cube->faces[face_idx];
 - 
 
308 
																		face.used = true;
 - 
 
309 
 
 - 
 
310 
																}
 - 
 
311 
														}
 310 
														}
312 
												}
 311 
													// Now add the face to the cube.
- 
 
 312 
													int face_idx  = dir_to_face(N6i[l]);
- 
 
 313 
													CubeFace& face = cube->faces[face_idx];
313 
										}
314 
													face.used = true;
314 
				}
 - 
 
315 
 
 315 
 
316 
 
 - 
 
317 
				/* This complicated function is used to connect the cube faces.
 - 
 
318 
					 If the faces are considered to be little square pieces of paper, this
 - 
 
319 
					 function glues them along the edges. However, it takes care to glue
 - 
 
320 
					 faces in such a way as to not create edges with more than two adjacent
 - 
 
321 
					 faces. */
 - 
 
322 
				template<class T>
 - 
 
323 
				void CubeGrid<T>::dual_cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m, bool connect)
 - 
 
324 
				{
 - 
 
325 
						/* For each cube in the table of cubes */
 - 
 
326 
						typename list<Cube>::iterator cube_iter = cube_table.begin();
 - 
 
327 
						for(;cube_iter != cube_table.end(); ++cube_iter)
 - 
 
328 
						{
 - 
 
329 
								/* For each of the six faces, if the face is used */
 - 
 
330 
								Cube& cube = *cube_iter;
 - 
 
331 
								Vec3i pi0 = cube.pos;
 - 
 
332 
								float val0 = voxel_grid[pi0];
 - 
 
333 
								for(int face_idx=0; face_idx<6; ++face_idx)
 - 
 
334 
								{
 - 
 
335 
										CubeFace& face = cube.faces[face_idx];
 - 
 
336 
										if(face.used)
 - 
 
337 
										{
 - 
 
338 
												Vec3i pi1(pi0);
 - 
 
339 
												pi1 += N6i[face_idx];
 - 
 
340 
												Vec3f p;
 - 
 
341 
												if(voxel_grid.in_domain(pi1))
 - 
 
342 
												{
 - 
 
343 
														float val1 = voxel_grid[pi1];
 - 
 
344 
														float t= (iso-val0)/(val1-val0);
 - 
 
345 
														p = Vec3f(pi0)*(1.0f-t)+Vec3f(pi1)*t;
 - 
 
346 
												}
 316 
												}
347 
												else
 - 
 
348 
														p = Vec3f(pi0) * 0.5 + Vec3f(pi1) * 0.5;
 - 
 
349 
												face.vert = m.create_vertex(p);
 - 
 
350 
												// We mark the vertex with the cube index.
 - 
 
351 
												face.vert->touched = reinterpret_cast<TouchType>(&cube);
 - 
 
352 
										}
 317 
										}
353 
								}
 318 
								}
354 
						}
 319 
						}
355
320 
		}
356 
						/* For each cube in the table of cubes */
 - 
 
357 
						cube_iter = cube_table.begin();
 - 
 
358 
						for(;cube_iter != cube_table.end(); ++cube_iter)
 - 
 
359 
						{
 - 
 
360 
								/* For each of the six faces, if the face is used */
 - 
 
361 
								Cube& cube = *cube_iter;
 - 
 
362 
								for(int i=0;i<6;++i)
 - 
 
363 
										if(cube.faces[i].used)
 - 
 
364 
										{
 - 
 
365 
												CubeFace& face = cube.faces[i];  // The face
 - 
 
366
- 
 
367 
												// Let e be one of the four edges of the face
 - 
 
368 
												for(int e=0;e<4;++e)
 - 
 
369 
												{
 - 
 
370 
														Cube* cube_n;  // Neighbouring cube
 - 
 
371 
														int i_n, e_n;
 - 
 
372 
 
 321 
 
373 
														// Do the actual gluing ...
 - 
 
374 
														if(!connect)
- 
 
375 
																find_neighbour_face_disconnect(cube, i, e, cube_n, i_n, e_n);
 - 
 
376 
														else
 - 
 
377 
																find_neighbour_face_connect(cube, i, e, cube_n, i_n, e_n);
 - 
 
378 
 
 - 
 
379 
														CubeFace& face_n = cube_n->faces[i_n];
 - 
 
380 
														HalfEdgeIter& he = face.get_he(m,e);
 - 
 
381 
														face.vert->he = he;
 - 
 
382 
														he->vert = face_n.vert;
 - 
 
383 
 
 - 
 
384 
														/* Link this halfedge to the _next_ halfedge on the
- 
 
385 
															 connected face. */
 - 
 
386 
														assert(face_n.used);
 - 
 
387 
														HalfEdgeIter& he_n = face_n.get_he(m, (e_n+3)%4);
 - 
 
388 
														link(he, he_n);
 - 
 
389 
 
 - 
 
390 
														/* Glue this halfedge to the corresponding one on the
 - 
 
391 
															 connected face.*/
 - 
 
392 
														HalfEdgeIter& he_opp = face_n.get_he(m, e_n);
 - 
 
393 
														face_n.vert->he = he_opp;
 - 
 
394 
														glue(he, he_opp);
 - 
 
395 
												}
 - 
 
396 
										}
 - 
 
397 
 
 322 
 
- 
 
 323 
		template<class T>
- 
 
 324 
		void CubeGrid<T>::cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m,
- 
 
 325 
																		 vector<Vec3f>& face_positions)
- 
 
 326 
		{
- 
 
 327 
			face_positions.clear();
- 
 
 328 
			Geometry::TrilinFilter<Geometry::RGrid<T> > inter(&voxel_grid);
- 
 
 329 
			typename list<Cube>::iterator cube_iter = cube_table.begin();
- 
 
 330 
 
- 
 
 331 
			/* For each cube in the table of cubes */
- 
 
 332 
			cube_iter = cube_table.begin();
- 
 
 333 
			for(;cube_iter != cube_table.end(); ++cube_iter)
- 
 
 334 
				{
- 
 
 335 
					/* For each of the six faces, if the face is used */
- 
 
 336 
					Cube& cube = *cube_iter;
- 
 
 337 
					for(int i=0;i<6;++i)
- 
 
 338 
						if(cube.faces[i].used)
- 
 
 339 
							{
- 
 
 340 
								CubeFace& face = cube.faces[i];  // The face
- 
 
 341 
								FaceIter hmesh_face = m.create_face();
- 
 
 342 
								hmesh_face->touched = face_positions.size();
- 
 
 343 
								Vec3i p0(cube.pos);
- 
 
 344 
								Vec3i p1 = p0 + N6i[i];
- 
 
 345 
								float v0 = voxel_grid[p0];
- 
 
 346 
								float v1 = voxel_grid[p1];
- 
 
 347 
								float t = (iso-v0)/(v1-v0);
- 
 
 348 
								face_positions.push_back(Vec3f(p0)*(1.0f-t)+Vec3f(p1)*t);
- 
 
 349
- 
 
 350 
								// Let e be one of the four edges of the face
- 
 
 351 
								for(int e=0;e<4;++e)
- 
 
 352 
									{
- 
 
 353 
										HalfEdgeIter he = face.get_he(m,e);
- 
 
 354 
										link(he,face.get_he(m,(e+1)%4));
- 
 
 355 
										link(face.get_he(m,(e+3)%4),he);
- 
 
 356 
										he->face = hmesh_face;
- 
 
 357 
 
- 
 
 358 
										HalfEdgeIter he_n = edge_to_glue(voxel_grid, cube, i, e);
- 
 
 359 
										if(he_n != NULL_HALFEDGE_ITER)
- 
 
 360 
											glue(he, he_n);
- 
 
 361 
									}
- 
 
 362 
								hmesh_face->last = face.get_he(m,0);
398 
						}
 363 
							}
399 
 
 364 
 
- 
 
 365 
				}
- 
 
 366 
 
400 
						HalfEdgeIter he0 = m.halfedges_begin();
 367 
			for(HalfEdgeIter h = m.halfedges_begin(); h != m.halfedges_end(); ++h)
- 
 
 368 
				{
401 
						while(he0 != m.halfedges_end())
 369 
					if (h->opp == NULL_HALFEDGE_ITER)
402 
						{
 370 
						{
403 
								// If this halfedge is not assigned to a face.
 - 
 
404 
								if(he0->face == NULL_FACE_ITER)
 - 
 
405 
								{
 - 
 
406 
										// Create a face and let the current halfedge belong to it.
 - 
 
407 
										FaceIter fi = m.create_face();
 371 
							HalfEdgeIter ho = m.create_halfedge();
408 
										fi->last = he0;
 - 
 
409 
										he0->face = fi;
 - 
 
410
- 
 
411 
										// Loop over all connected halfedges and let them also belong
- 
 
412 
										// to the face.
 - 
 
413 
										HalfEdgeIter he = he0->next;
 - 
 
414 
										while(he != he0)
 - 
 
415 
										{
 - 
 
416 
												he->face = fi;
 - 
 
417 
												he = he->next;
 - 
 
418 
										}
 - 
 
419 
							}
 - 
 
420 
								++he0;
 372 
							glue(ho,h);
421 
						}
 373 
						}
422 
				}
 374 
				}
423 
 
 375 
 
424 
				template<class T>
 376 
			cube_iter = cube_table.begin();
425 
				void CubeGrid<T>::cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m, bool connect)
 377 
			for(;cube_iter != cube_table.end(); ++cube_iter)
426 
				{
 378 
				{
427 
						typename list<Cube>::iterator cube_iter = cube_table.begin();
 379 
					//For each of the six faces, if the face is used
428 
 
 - 
 
429 
						/* For each cube in the table of cubes */
 380 
					Cube& cube = *cube_iter;
430 
						cube_iter = cube_table.begin();
 381 
					for(int i=0;i<6;++i)
431 
						for(;cube_iter != cube_table.end(); ++cube_iter)
 382 
						if(cube.faces[i].used)
432 
						{
 383 
							{
433 
								/* For each of the six faces, if the face is used */
 384 
								CubeFace& face = cube.faces[i];  // The face
434 
								Cube& cube = *cube_iter;
 385 
								// Let e be one of the four edges of the face
435 
								for(int i=0;i<6;++i)
 386 
								for(int e=0;e<4;++e)
436 
										if(cube.faces[i].used)
 387 
									if(face.he[e]->vert == NULL_VERTEX_ITER)
437 
										{
 388 
										{
438 
												CubeFace& face = cube.faces[i];  // The face
 389 
											HalfEdgeIter last = face.he[e];
439 
												FaceIter hmesh_face = m.create_face();
 390 
											HalfEdgeIter he = last;
440 
												// Let e be one of the four edges of the face
 391 
											while(he->opp->face != NULL_FACE_ITER)
441 
												for(int e=0;e<4;++e)
 - 
 
442 
												{
 392 
												{
443 
														Cube* cube_n;  // Neighbouring cube
 - 
 
444 
														int i_n, e_n;
 393 
													he = he->opp->prev;
445 
 
 - 
 
446 
														// Do the actual gluing ...
 - 
 
447 
														if(!connect)
394 
													if(he == last)
448 
																find_neighbour_face_disconnect(cube, i, e, cube_n, i_n, e_n);
 - 
 
449 
														else
 395 
														break;
450 
																find_neighbour_face_connect(cube, i, e, cube_n, i_n, e_n);
 - 
 
451
- 
 
452 
														CubeFace& face_n = cube_n->faces[i_n];
 - 
 
453 
														HalfEdgeIter he = face.get_he(m,e);
 - 
 
454 
														HalfEdgeIter he_n = face_n.get_he(m,e_n);
 - 
 
455
- 
 
456 
														link(he,face.get_he(m,(e+1)%4));
 - 
 
457 
														link(face.get_he(m,(e+3)%4),he);
 - 
 
458 
														glue(he, he_n);
 - 
 
459 
														he->face = hmesh_face;
 - 
 
460 
												}
 396 
												}
- 
 
 397 
											last = he;
- 
 
 398
- 
 
 399 
											Vec3f pos = cube_corner(cube.pos, i, e);
461 
												hmesh_face->last = face.get_he(m,0);
 400 
											VertexIter vert = m.create_vertex(pos);
- 
 
 401
462 
 
 402 
											do
463 
 
 403 
												{
- 
 
 404 
													he->vert = vert;
- 
 
 405 
													he = he->next->opp;
464 
										}
 406 
												}
- 
 
 407 
											while(he->face != NULL_FACE_ITER && he != last);
465 
 
 408 
 
466 
						}
 - 
 
467 
						cube_iter = cube_table.begin();
 - 
 
468 
						for(;cube_iter != cube_table.end(); ++cube_iter)
 - 
 
469 
						{
 - 
 
470 
								/* For each of the six faces, if the face is used */
 - 
 
471 
								Cube& cube = *cube_iter;
 - 
 
472 
								for(int i=0;i<6;++i)
 - 
 
473 
										if(cube.faces[i].used)
 - 
 
474 
										{
 - 
 
475 
												CubeFace& face = cube.faces[i];  // The face
 - 
 
476 
												// Let e be one of the four edges of the face
 - 
 
477 
												for(int e=0;e<4;++e)
 - 
 
478 
														if(face.he[e]->vert == NULL_VERTEX_ITER)
 409 
											if(he->face == NULL_FACE_ITER)
479 
														{
 410 
												{
480 
																Vec3i to_face = N6i[i];
 411 
													he->vert = vert;
481 
																Vec3i along_edge = edges[i][e];
- 
 
482 
																Vec3i to_edge = cross(along_edge,to_face);
 - 
 
483 
																Vec3f pos = Vec3f(cube.pos) + 0.5*Vec3f(to_face + along_edge + to_edge);
 - 
 
484 
																HalfEdgeIter last = face.he[e];
 - 
 
485 
																VertexIter vi = m.create_vertex(pos);
 - 
 
486 
																vi->he = last->opp;
 412 
													link(he,last->opp);
487 
																HalfEdgeIter he = last;
 - 
 
488 
																do
 413 
												}
489 
																{
 414
490 
																		he->vert = vi;
 - 
 
491 
																		he = he->next->opp;
 415 
											vert->he = last->opp;
492 
																}
 - 
 
493 
																while(he != last);
 - 
 
494 
														}
 - 
 
495 
										}
 416 
										}
496 
						}
 417 
							}
497 
				}
 418 
				}
- 
 
 419 
		}
- 
 
 420 
	}
498 
 
 421
499 
				/* A fairly complex algorithm which triangulates faces. It triangulates
 422 
		/* A fairly complex algorithm which triangulates faces. It triangulates
500 
					 by iteratively splitting faces. It always splits a face by picking the
 423 
			 by iteratively splitting faces. It always splits a face by picking the
501 
					 edge whose midpoint is closest to the isovalue. */
 424 
			 edge whose midpoint is closest to the isovalue. */
502
425
503 
				template<class T>
 426 
		template<class T>
504 
				void triangulate(const RGrid<T>& voxel_grid, Manifold& m, float iso)
 427 
		void triangulate(const RGrid<T>& voxel_grid, Manifold& m, float iso)
505 
				{
 428 
		{
506 
#ifndef NDEBUG
 429 
#ifndef NDEBUG
507 
						cout << "Initial mesh valid : " << m.is_valid() << endl;
 430 
			cout << "Initial mesh valid : " << m.is_valid() << endl;
508 
#endif
 431 
#endif
509 
						Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(&voxel_grid);
 432 
			Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(&voxel_grid);
510 
						int work;
 433 
			int work;
511 
						do
 434 
			do
512 
						{
 435 
				{
513 
								// For every face.
 436 
					// For every face.
514 
								work = 0;
 437 
					work = 0;
515 
								for(FaceIter fi =m.faces_begin(); fi != m.faces_end(); ++fi)
 438 
					for(FaceIter fi =m.faces_begin(); fi != m.faces_end(); ++fi)
516 
								{
 439 
						{
517 
										// Create a vector of vertices.
 440 
							// Create a vector of vertices.
518 
										vector<VertexIter> verts;
 441 
							vector<VertexIter> verts;
519 
										for(FaceCirculator fc(fi); !fc.end(); ++fc)
 442 
							for(FaceCirculator fc(fi); !fc.end(); ++fc)
520 
												verts.push_back(fc.get_vertex());
 443 
								verts.push_back(fc.get_vertex());
521
444
522 
										// If there are just three we are done.
 445 
							// If there are just three we are done.
523 
										if(verts.size() == 3) continue;
 446 
							if(verts.size() == 3) continue;
524
447
525 
										// Find vertex pairs that may be connected.
 448 
							// Find vertex pairs that may be connected.
526 
										vector<pair<int,int> > vpairs;
 449 
							vector<pair<int,int> > vpairs;
527 
										const int N = verts.size();
 450 
							const int N = verts.size();
528 
										for(int i=0;i<N-2;++i)
 451 
							for(int i=0;i<N-2;++i)
529 
												for(int j=i+2;j<N; ++j)
 452 
								for(int j=i+2;j<N; ++j)
530 
														if(!is_connected(verts[i], verts[j]))
 453 
									if(!is_connected(verts[i], verts[j]))
531 
																vpairs.push_back(pair<int,int>(i,j));
 454 
										vpairs.push_back(pair<int,int>(i,j));
532
455
533 
										assert(!vpairs.empty());
 456 
							if(vpairs.empty())
- 
 
 457 
								{
- 
 
 458 
									cout << "Warning: could not triangulate a face."
- 
 
 459 
											 << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;
- 
 
 460 
									continue;
- 
 
 461 
								}
534
462
535 
										/* For all vertex pairs, find the absolute value of iso-x
 463 
							/* For all vertex pairs, find the absolute value of iso-x
536 
											 where x is the interpolated value of the volume at the
464 
								 where x is the interpolated value of the volume at the
537 
											 midpoint of the line segment connecting the two vertices.*/
 465 
								 midpoint of the line segment connecting the two vertices.*/
538
466
539 
										float min_val=FLT_MAX;
 467 
							float min_val=FLT_MAX;
540 
										int min_k = -1;
 468 
							int min_k = -1;
541 
										for(size_t k=0;k<vpairs.size(); ++k)
 469 
							for(size_t k=0;k<vpairs.size(); ++k)
542 
										{
 470 
								{
543 
												int i = vpairs[k].first;
 471 
									int i = vpairs[k].first;
544 
												int j = vpairs[k].second;
 472 
									int j = vpairs[k].second;
545 
												Vec3f centre =
473 
									Vec3f centre =
546 
														(verts[i]->pos +
474 
										(verts[i]->pos +
547 
														 verts[j]->pos)/2.0f;
 475 
										 verts[j]->pos)/2.0f;
548
476
549 
												float val;
 477 
									float val;
550 
												if(grid_interp.in_domain(centre))
 478 
									if(grid_interp.in_domain(centre))
551 
														val = fabs(grid_interp(centre)-iso);
 479 
										val = fabs(grid_interp(centre)-iso);
552 
												else
 480 
									else
553 
														val = 0.0f;
 481 
										val = 0.0f;
554 
												if(val<min_val)
 482 
									if(val<min_val)
555 
												{
 - 
 
556 
														min_val = val;
 - 
 
557 
														min_k = k;
 - 
 
558 
												}
 - 
 
559 
										}
 - 
 
560 
										assert(min_k != -1);
 - 
 
561
- 
 
562 
										{
 483 
										{
563 
												// Split faces along edge whose midpoint is closest to isovalue
 - 
 
564 
												int i = vpairs[min_k].first;
 484 
											min_val = val;
565 
												int j = vpairs[min_k].second;
 485 
											min_k = k;
566 
												m.split_face(fi, verts[i], verts[j]);
 - 
 
567 
										}
 486 
										}
568 
										++work;
 - 
 
569
- 
 
570 
								}
 487 
								}
- 
 
 488 
							assert(min_k != -1);
- 
 
 489
- 
 
 490 
							{
- 
 
 491 
								// Split faces along edge whose midpoint is closest to isovalue
- 
 
 492 
								int i = vpairs[min_k].first;
- 
 
 493 
								int j = vpairs[min_k].second;
- 
 
 494 
								m.split_face(fi, verts[i], verts[j]);
- 
 
 495 
							}
- 
 
 496 
							++work;
- 
 
 497
571 
						}
 498 
						}
572 
						while(work);
 - 
 
573 
 
 - 
 
574 
#ifndef NDEBUG
 - 
 
575 
						cout << "Valid after triangulation : " << m.is_valid() << endl;
 - 
 
576 
#endif
 - 
 
577 
				}
 499 
				}
- 
 
 500 
			while(work);
578 
 
 501 
 
579 
				/* This function collapses every edge whose endpoints belong to the same
- 
 
580 
					 cube if this edge can be removed without violating the manifold condition*/
 - 
 
581 
				void remove_duplicates(Manifold& m)
 - 
 
582 
				{
 - 
 
583 
						// BUG --- we never do this loop more than one time.
 - 
 
584 
						// does it matter.
 - 
 
585 
						VertexIter vi = m.vertices_begin();
 - 
 
586 
						bool did_work = false;
 - 
 
587 
						do
 - 
 
588 
						{
 - 
 
589 
								did_work = false;
 - 
 
590 
								while(vi != m.vertices_end())
 - 
 
591 
								{
 - 
 
592 
										bool did = false;
 - 
 
593 
										VertexCirculator vc(vi);
 - 
 
594 
										int k=0;
 - 
 
595 
										for(;!vc.end();++vc)
 - 
 
596 
										{
 - 
 
597 
												assert(++k<1000);
 - 
 
598 
												HalfEdgeIter he = vc.get_halfedge();
 - 
 
599 
												VertexIter n = vc.get_vertex();
 - 
 
600 
												if(n->touched == vi->touched && m.collapse_precond(he))
 - 
 
601 
												{
 - 
 
602 
														++vi;
 - 
 
603 
														m.collapse_halfedge(he,true);
 - 
 
604 
														did = did_work = true;
 - 
 
605 
														break;
 - 
 
606 
												}
 - 
 
607 
										}
 - 
 
608 
										if(!did) ++vi;
 - 
 
609 
								}
 - 
 
610 
 
 - 
 
611 
						}
 - 
 
612 
						while(did_work);
 - 
 
613 
#ifndef NDEBUG
 502 
#ifndef NDEBUG
614 
						cout << "Valid after removal of excess vertices : " << m.is_valid() << endl;
 503 
			cout << "Valid after triangulation : " << m.is_valid() << endl;
615 
#endif
 504 
#endif
616 
				}
 - 
 
617 
 
 - 
 
618 
				/* A function (not used now) which pushes vertices onto the isosurface. */
 - 
 
619 
				template<class T>
 - 
 
620 
				void push_vertices(const RGrid<T>& voxel_grid, Manifold& m, float iso)
 - 
 
621 
				{
 - 
 
622 
						GradientFilter<RGrid<T> > grad(&voxel_grid);
 - 
 
623 
						TrilinFilter<GradientFilter<RGrid<T> > > ginter(&grad);
 - 
 
624 
						TrilinFilter<RGrid<T> > inter(&voxel_grid);
 - 
 
625 
						for(int i=0;i<4;++i)
 - 
 
626 
						{
 - 
 
627 
								for(VertexIter vi= m.vertices_begin();
 - 
 
628 
										vi != m.vertices_end(); ++vi)
 - 
 
629 
								{
 - 
 
630 
										Vec3f p = vi->pos;
 - 
 
631 
										if(ginter.in_domain(p))
 - 
 
632 
										{
 - 
 
633 
												Vec3f g = ginter(p);
 - 
 
634 
												float s = sqr_length(g);
 - 
 
635 
												if(s>0.01)
 - 
 
636 
												{
 - 
 
637 
														float v = inter(p)-iso;
 - 
 
638 
														vi->pos = (p-g*v/s);
 - 
 
639 
												}
 - 
 
640 
										}
 - 
 
641 
								}
 - 
 
642 
						}
 - 
 
643 
				}
 - 
 
644 
 
 - 
 
645 
		}
 505 
		}
646 
 
 506 
 
- 
 
 507
647 
		template<typename T>
 508 
	template<typename T>
648 
		void fair_polygonize(const RGrid<T>& voxel_grid,
509 
	void fair_polygonize(const RGrid<T>& voxel_grid,
649 
												 Manifold& mani,
510 
											 Manifold& mani,
650 
												 float iso,
 511 
											 float iso)
- 
 
 512 
	{
- 
 
 513 
		CubeGrid<T> cube_grid(voxel_grid, iso);
- 
 
 514 
		vector<Vec3f> face_positions;
- 
 
 515 
		cube_grid.cuberille_mesh(voxel_grid, mani, face_positions);
- 
 
 516 
		cout << "Initial mesh valid : " << mani.is_valid() << endl;
- 
 
 517 
		compute_dual(mani, face_positions);
- 
 
 518 
	}
- 
 
 519
- 
 
 520 
	template<typename T>
- 
 
 521 
	void cuberille_polygonize(const RGrid<T>& voxel_grid,
651 
												 bool invert_contour,
522 
														Manifold& mani,
652 
												 bool connect,
 523 
														float iso,
653 
												 bool push)
 524 
														bool push)
654 
		{
 525 
	{
655 
				CubeGrid<T> cube_grid(voxel_grid, iso, invert_contour);
 526 
		CubeGrid<T> cube_grid(voxel_grid, iso);
- 
 
 527 
		vector<Vec3f> face_positions;
656 
				cube_grid.dual_cuberille_mesh(voxel_grid, mani, connect);
 528 
		cube_grid.cuberille_mesh(voxel_grid, mani, face_positions);
- 
 
 529 
		if(push)
- 
 
 530 
			{
- 
 
 531 
				for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)
- 
 
 532 
					{
657 
 				triangulate(voxel_grid, mani, iso);
 533 
						VertexCirculator vc(v);
- 
 
 534 
						Vec3f p(0);
- 
 
 535 
						int n=0;
658 
 				remove_duplicates(mani);
 536 
						while(!vc.end())
- 
 
 537 
							{
- 
 
 538 
								if(vc.get_face() != NULL_FACE_ITER)
- 
 
 539 
									{
659 
				if(push) push_vertices(voxel_grid, mani, iso);
 540 
										p += face_positions[vc.get_face()->touched];
- 
 
 541 
										++n;
- 
 
 542 
									}
- 
 
 543 
								++vc;
- 
 
 544 
							}
- 
 
 545 
						v->pos = p/float(n);
- 
 
 546 
					}
660 
		}
 547 
			}
- 
 
 548 
	}
- 
 
 549 
 
- 
 
 550 
	template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,
- 
 
 551 
																							 Manifold&, float);
- 
 
 552 
	template void fair_polygonize<float>(const RGrid<float>&,
- 
 
 553 
																			 Manifold&, float);
- 
 
 554 
 
- 
 
 555 
	template void cuberille_polygonize<unsigned char>(const RGrid<unsigned char>&,
- 
 
 556 
																										Manifold&, float, bool);
- 
 
 557 
	template void cuberille_polygonize<float>(const RGrid<float>&,
- 
 
 558 
																						Manifold&, float, bool);
661 
 
 559 
 
662 
		template<typename T>
 - 
 
663 
		void cuberille_polygonize(const RGrid<T>& voxel_grid,
560 
	template void triangulate(const RGrid<unsigned char>& voxel_grid, Manifold& m, float iso);
664 
															Manifold& mani,
- 
 
665 
															float iso,
 - 
 
666 
															bool invert_contour,
- 
 
667 
															bool connect,
 - 
 
668 
															bool push)
 - 
 
669 
		{
 - 
 
670 
				CubeGrid<T> cube_grid(voxel_grid, iso, invert_contour);
 - 
 
671 
				cube_grid.cuberille_mesh(voxel_grid, mani, connect);
 - 
 
672 
				if(push) push_vertices(voxel_grid, mani, iso);
 - 
 
673 
		}
 - 
 
674 
 
 561 
 
675 
		template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,
 - 
 
676 
																								 Manifold&, float, bool,bool,bool);
 - 
 
677 
		template void fair_polygonize<float>(const RGrid<float>&,
 - 
 
678 
																				 Manifold&, float, bool,bool,bool);
 - 
 
679 
 
 - 
 
680 
		template void cuberille_polygonize<unsigned char>(const RGrid<unsigned char>&,
 - 
 
681 
																											Manifold&, float, bool,bool,bool);
 - 
 
682 
		template void cuberille_polygonize<float>(const RGrid<float>&,
 562 
	template void triangulate(const RGrid<float>& voxel_grid, Manifold& m, float iso);
683 
																							Manifold&, float, bool,bool,bool);
 - 
 
684 
}
 563 
}
- 
 
 564
685 
 
 565 
 
686 
 
 566 
 
- 
 
 567
- 
 
 568