WebSVN – gelsvn – Diff – /trunk/src/HMeshUtil/fair_polygonize.cpp


#include <cfloat>
#include <algorithm>
#include <fstream>
#include <iostream>
#include <list>
 
#include "CGLA/Mat2x3f.h"
#include "CGLA/Mat2x2f.h"
#include "HMesh/Manifold.h"
#include "HMesh/VertexCirculator.h"
#include "HMesh/FaceCirculator.h"
 
 
#include "Geometry/RGrid.h"
#include "Geometry/TrilinFilter.h"
#include "Geometry/GradientFilter.h"
#include "Geometry/Neighbours.h"
#include "Util/HashTable.h"
#include "Util/HashKey.h"
 
#include "fair_polygonize.h"
 
namespace Geometry
{
		using namespace HMesh;
		using namespace Util;
		using namespace CGLA;
		using namespace std;
 
		namespace 
		{
				/*
					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;
				}
		
				/* 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;
				}
			
				/* 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;
						VertexIter vert;
						HalfEdgeIter he[4];
				
						CubeFace();
 
						HalfEdgeIter& get_he(Manifold* m, int i)
								{
										if(he[i] == NULL_HALFEDGE_ITER)
												he[i] = m->create_halfedge();
										return he[i];
								}
				};
 
				CubeFace::CubeFace(): used(false), vert(NULL_VERTEX_ITER)
				{
						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 Mesher
				{
						const RGrid<T>* const voxel_grid;        // The voxel grid we work on
						Manifold* m;		                         // Output manifold
						HashTable<HashKey3usi,Cube*> cube_hash;  // Hashtable of cube pointers
						list<Cube> cube_table;                   // Linked list containing
						// the actual cubes.
 
				public:
			
						Mesher(const RGrid<T>* _voxel_grid, // input voxel grid
									 Manifold* _m,                // output HMesh
									 float iso,                   // iso value
									 bool inv_cont);              // invert contour?
 
						void process_cubes();
 
						void create_faces();
 
						void triangulate(float iso);
 
						void remove_duplicates();
 
						void push_vertices(float iso);
 
				};
	
				template<class T>
				inline bool comp(bool d, T a, T b)
				{
						return d? (a<b) : (b<a); 
				}
		
				template<class T>
				Mesher<T>::Mesher(const RGrid<T>* _voxel_grid, 
													Manifold* _m, float iso, bool inv_cont): 
						voxel_grid(_voxel_grid),
						m(_m)
				{
						/* 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 isovale (or below if we
													 are finding the inverse contour) */
												if(comp(inv_cont,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 || 
																	 !comp(inv_cont,float((*voxel_grid)[pi1]),iso))
																{
																		/* Create a cube face for the interface between
																			 our voxel and this neighbour. */
 
																		float val1 = iso;
																		if(indom)
																				val1 = (*voxel_grid)[pi1];
																		
																		/* 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;
 
																		// Finally, we place a vertex corresponding to the
																		// centre of the cube face. The vertex is 
																		// positioned on the line between the two voxels
																		// (which pierces the new face) at the point where
																		// the interpolated value = isovalue
																		float t= (iso-val0)/(val1-val0);
																		Vec3f p = Vec3f(pi0)*(1.0f-t)+Vec3f(pi1)*t;
																		face.vert = m->create_vertex(p);
																		
																		// We mark the vertex with the cube index.
																		face.vert->touched = reinterpret_cast<int>(cube);
																}
														}
												}
										}		
				}
 
 
				/* This complicated function is used to connect the cube faces.
					 If the faces are cosidered to be little square pieces of paper, this
					 function glues them along the edges. However, it takes care to glue
					 faces in such a way as to not create edges with more the two adjacent
					 faces. */
				template<class T>
				void Mesher<T>::process_cubes()
				{
						/* For each cube in the table of cubes */
						typename list<Cube>::iterator 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)
										{
												Vec3i pos = cube.pos;
												Vec3i dir = N6i[i];
 
												// Let e be one of the four edges of the face
												for(int e=0;e<4;++e)
												{
														Vec3i edge = edges[i][e];	       // Vector _along_ edge
														Vec3i to_edge = cross(edge,dir); // vector to edge
														CubeFace& face = cube.faces[i];  // The face
 
														Cube* cube_n;  // Neighbouring cube
														Vec3i dir_n;   // direction from center to neighbour face
														Vec3i edge_n; // Edge vector for neighbour face
 
 
														/*
															Algorithm: Pick faces to glue
 
															Algorithm which connects cube faces. If both A and B
															are above the isovalue and 0 and C are below, A and B
															will be seperated, i.e. the 0A and AC face will be 
															connected. There is a different out-commented algorithm
															at the end of this file which connects AO and 0B
															instead.
															
															Figure.
 
															0|B
															-+-
															A|C
															
															Explanation:
														
															Let the current cube be A. Since we are here, the voxel
															of A is above the isovalue, and the voxel of O is below
															since there is a face AO.
 
															We first check if there is a face between A and C. 
															if that is the case, the AO face is connected to the AC 
															face. 
 
															Otherwise, we check if B is also inside the isocontour.
															If that is true there is an OB face, and we connect the
															A0 and the 0B faces.
															
															Failing that, there must be a CB face, and we connect
															with that.
														*/
										
														// Begin "Pick faces to glue"
														Cube* dummy;
														if(cube.faces[dir_to_face(to_edge)].used)
														{
																dir_n = to_edge;
																edge_n = - edge; 
																cube_n = &cube;  
														}
														else if(!cube_hash.find_entry(HashKey3usi(pos+dir+to_edge),
																													dummy))
														{
																assert(sqr_length(to_edge)==1);
																dir_n = dir;
																edge_n = - edge;
																cube_hash.find_entry(HashKey3usi(pos + to_edge), 
																										 cube_n);
														}
														else 
														{
																dir_n = - to_edge;
																edge_n = - edge;
																cube_n = dummy;
														}
														// End "Pick faces to glue"
 
														// Do the actual gluing ...
														
														/* Let the vertex of this face's halfedge be the one on
															 the connected face. */
														int i_n = dir_to_face(dir_n);
														CubeFace& face_n = cube_n->faces[i_n];
														HalfEdgeIter& he = face.get_he(m,e);
														face.vert->he = he;
														he->vert = face_n.vert;
 
														/* Link this halfedge to the _next_ halfedge on the 
															 connected face. */
														int e_n = dir_to_edge(i_n, edge_n);
														assert(face_n.used);
														HalfEdgeIter& he_n = face_n.get_he(m, (e_n+3)%4);
														link(he, he_n);
 
														/* Glue this halfedge to the corresponding one on the
															 connected face.*/
														HalfEdgeIter& he_opp = face_n.get_he(m, e_n);
														face_n.vert->he = he_opp;
														glue(he, he_opp);
												}
										}
 
						}
				}
 
				template<class T>
				void Mesher<T>::create_faces()
				{
						// For each halfedge
						HalfEdgeIter he0 = m->halfedges_begin();
						while(he0 != m->halfedges_end())
						{
								// If this halfedge is not assigned to a face.
								if(he0->face == NULL_FACE_ITER)
								{
										// Create a face and let the current halfedge belong to it.
										FaceIter fi = m->create_face();
										fi->last = he0;
										he0->face = fi;
 
										// Loop over all connected halfedges and let them also belong 
										// to the face.
										HalfEdgeIter he = he0->next;
										while(he != he0)
										{
												he->face = fi;
												he = he->next;
										}
								}
								++he0;
						}
				}
 
 
				/* A fairly complex algorithm which triangulates faces. It triangulates
					 by iteratively splitting faces. It always splits a face by picking the
					 edge whose midpoint is closest to the isovalue. */
					
				template<class T>
				void Mesher<T>::triangulate(float iso)
				{
#ifndef NDEBUG
						cout << "Initial mesh valid : " << m->is_valid() << endl;
#endif
						Geometry::TrilinFilter<Geometry::RGrid<T> > grid_interp(voxel_grid);
						int work;
						do
						{
								// For every face.
								work = 0;
								for(FaceIter fi =m->faces_begin(); fi != m->faces_end(); ++fi)
								{
										// Create a vector of vertices.
										vector<VertexIter> verts;
										for(FaceCirculator fc(fi); !fc.end(); ++fc)
												verts.push_back(fc.get_vertex());
								
										// If there are just three we are done.
										if(verts.size() == 3) continue;
										
										// Find vertex pairs that may be connected.
										vector<pair<int,int> > vpairs;
										const int N = verts.size();
										for(int i=0;i<N-2;++i)
												for(int j=i+2;j<N; ++j)
														if(!is_connected(verts[i], verts[j]))
																vpairs.push_back(pair<int,int>(i,j));
										
										assert(!vpairs.empty());
									
										/* For all vertex pairs, find the absolute value of iso-x
											 where x is the interpolated value of the volume at the 
											 midpoint of the line segment connecting the two vertices.*/
											 
										float min_val=FLT_MAX;
										int min_k = -1;
										for(int k=0;k<vpairs.size(); ++k)
										{
												int i = vpairs[k].first;
												int j = vpairs[k].second;
												Vec3f centre = 
														(verts[i]->pos + 
														 verts[j]->pos)/2.0f;
												
												float val;
												if(grid_interp.in_domain(centre))
														val = fabs(grid_interp(centre)-iso);
												else
														val = 0.0f;
												if(val<min_val)
												{
														min_val = val;
														min_k = k;
												}
										}
										assert(min_k != -1);
 
										// Split faces along edge whose midpoint is closest to isovalue
										int i = vpairs[min_k].first;
										int j = vpairs[min_k].second;
										m->split_face(fi, verts[i], verts[j]);
										++work;
										
								}
						}
						while(work);
 
#ifndef NDEBUG
						cout << "Valid after triangulation : " << m->is_valid() << endl;
#endif
				}
 
				/* This function collapses every edge whose endpoints belong to the same 
					 cube if this edge can be removed without violating the manifold condition*/
				template<class T>
				void Mesher<T>::remove_duplicates()
				{
						// BUG --- we never do this loop more than one time.
						// does it matter.
						VertexIter vi = m->vertices_begin();
						bool did_work = false;
						do
						{
								did_work = false;
								while(vi != m->vertices_end())
								{
										bool did = false;
										VertexCirculator vc(vi);
										int k=0;
										for(;!vc.end();++vc)
										{
												assert(++k<1000);
												HalfEdgeIter he = vc.get_halfedge();
												VertexIter n = vc.get_vertex();
												if(n->touched == vi->touched && m->collapse_precond(vi,he))
												{
														VertexIter vi2 = vi;
														++vi;
														m->collapse_halfedge(vi2,he,true);
														did = did_work = true;
														break;
												}
										}
										if(!did) ++vi;
								}
 
						}
						while(did_work);
#ifndef NDEBUG
						cout << "Valid after removal of excess vertices : " << m->is_valid() << endl;
#endif
				}
		
				/* A function (not used now) which pushes vertices onto the isosurface. */
				template<class T>
				void Mesher<T>::push_vertices(float iso)
				{
						GradientFilter<RGrid<T> > grad(voxel_grid);
						TrilinFilter<GradientFilter<RGrid<T> > > ginter(&grad);
						TrilinFilter<RGrid<T> > inter(voxel_grid);
		
						cout << "Pushing vertices" << endl;
						for(int i=0;i<4;++i)
						{
								cout << "." << flush;
								for(VertexIter vi= m->vertices_begin();
										vi != m->vertices_end(); ++vi)
								{
										Vec3f p = vi->get_pos();
										if(ginter.in_domain(p))
										{
												Vec3f g = ginter(p);
												float s = sqr_length(g);
												if(s>0.0001)
												{
														float v = inter(p)-iso;
														vi->set_pos(p-g*v/s);
												}
										}
								}
						}
				}
 
		}
 
		template<typename T>
		void fair_polygonize(const RGrid<T>& voxel_grid, 
												 Manifold& mani, 
												 float iso,
												 bool invert_contour)
		{
				Mesher<T> m(&voxel_grid, &mani, iso, invert_contour);
				m.process_cubes();
				m.create_faces();
				m.triangulate(iso);
				m.remove_duplicates();
				//m.push_vertices(iso);
		}
 
		template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,
																								 Manifold&, float, bool);
		template void fair_polygonize<float>(const RGrid<float>&,
																				 Manifold&, float, bool);
}
 
 
// --- Alternative "Pick faces to glue" code.
//                          if(cube_hash.find_entry(HashKey3usi(pos + dir + to_edge),
//                                        cube_n))
// 														{
// 																assert(sqr_length(dir + to_edge)==2);
// 																dir_n = - to_edge;
// 																edge_n = - edge;
// 														}
// 														else if(cube_hash.find_entry(HashKey3usi(pos + to_edge), 
// 																												 cube_n))
// 														{
// 																assert(sqr_length(to_edge)==1);
// 																dir_n = dir;
// 																edge_n = - edge;
// 														}
// 														else
// 														{
// 																dir_n = to_edge;
// 																edge_n = - edge;
// 																											cube_n = &cube;
// 														}
 
 
Generated by GNU Enscript 1.6.6.
 
 
 

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

Subversion Repositories gelsvn

(root)/trunk/src/HMeshUtil/fair_polygonize.cpp @ 94 – Rev 85 → 89