Subversion Repositories gelsvn

#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 HMesh

{

		using namespace Geometry;

		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;

						return -1;

				}

				/* Convert a face index and an edge vector to an edge index. */

				int dir_to_edge(int face, const Vec3i& edge)

				{

						for(int i=0;i<4;++i)

								if(edges[face][i] == edge) return i;

						assert(0);

						return -1;

				}

				/* 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 CubeGrid

				{

						float iso;

						HashTable<HashKey3usi,Cube*> cube_hash;  // Hashtable of cube pointers

						list<Cube> cube_table;                   // Linked list containing

                                          					 // the actual cubes.

				public:

						CubeGrid(const RGrid<T>& _voxel_grid, // input voxel grid

										 float iso,                   // iso value

										 bool inv_cont);              // invert contour?

						void find_neighbour_face_disconnect(Cube& cube, int face_idx, int edge_idx,

																								Cube*& cube_n, int& i_n, int& e_n);

						void find_neighbour_face_connect(Cube& cube, int face_idx, int edge_idx,

																						 Cube*& cube_n, int& i_n, int& e_n);

						void dual_cuberille_mesh(const RGrid<T>&, Manifold& m, bool connect=false);

						void cuberille_mesh(const RGrid<T>&, Manifold& m, bool connect=false);

				};

				template<class T>

				void CubeGrid<T>::find_neighbour_face_disconnect(Cube& cube, int face_idx, int edge_idx,

																												 Cube*& cube_n,  // Neighbouring cube

																												 int& i_n,  // direction from center to neighbour face

																												 int& e_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 faces 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.

						*/

						Vec3i pos = cube.pos;

						Vec3i dir = N6i[face_idx];

						Vec3i dir_n;   // direction from center to neighbour face

						Vec3i edge_n; // Edge vector for neighbour face

						Vec3i edge = edges[face_idx][edge_idx];	       // Vector _along_ edge

						Vec3i to_edge = cross(edge,dir); // vector to edge

						// 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;

						}

						i_n = dir_to_face(dir_n);

						e_n = dir_to_edge(i_n, edge_n);

				}

				template<class T>

				void CubeGrid<T>::find_neighbour_face_connect(Cube& cube, int face_idx, int edge_idx,

																											Cube*& cube_n,  // Neighbouring cube

																											int& i_n,  // direction from center to neighbour face

																											int& e_n) // Edge vector for neighbour face

				{

						Vec3i pos = cube.pos;

						Vec3i dir = N6i[face_idx];

						Vec3i dir_n;   // direction from center to neighbour face

						Vec3i edge_n; // Edge vector for neighbour face

						Vec3i edge = edges[face_idx][edge_idx];	       // Vector _along_ edge

						Vec3i to_edge = cross(edge,dir); // vector to edge

						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;

						}

						i_n = dir_to_face(dir_n);

						e_n = dir_to_edge(i_n, edge_n);

				}

				template<class T>

				inline bool comp(bool d, T a, T b)

				{

						return d? (a<b) : (b<a); 

				}

				template<class T>

				CubeGrid<T>::CubeGrid(const RGrid<T>& voxel_grid, float _iso, bool inv_cont): iso(_iso)

				{

						/* Loop over all voxels */

						for(int k=0;k<voxel_grid.get_dims()[ 2 ];++k)

								for(int j=0;j<voxel_grid.get_dims()[ 1 ];++j)

										for(int i=0;i<voxel_grid.get_dims()[ 0 ];++i)

										{

												Vec3i pi0(i,j,k);

												float val0 = voxel_grid[pi0];

												/* If the voxel value is above the 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))

																{

																		/* 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;

																}

														}

												}

										}		

				}

				/* This complicated function is used to connect the cube faces.

					 If the faces are considered 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 than two adjacent

					 faces. */

				template<class T>

				void CubeGrid<T>::dual_cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m, bool connect)

				{

						/* 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;

								Vec3i pi0 = cube.pos;

								float val0 = voxel_grid[pi0];

								for(int face_idx=0; face_idx<6; ++face_idx)

								{

										CubeFace& face = cube.faces[face_idx];

										if(face.used)

										{

												Vec3i pi1(pi0);

												pi1 += N6i[face_idx];

												Vec3f p;

												if(voxel_grid.in_domain(pi1))

												{

														float val1 = voxel_grid[pi1];

														float t= (iso-val0)/(val1-val0);

														p = Vec3f(pi0)*(1.0f-t)+Vec3f(pi1)*t;

												}

												else

														p = Vec3f(pi0) * 0.5 + Vec3f(pi1) * 0.5;

												face.vert = m.create_vertex(p);

												// We mark the vertex with the cube index.

												face.vert->touched = reinterpret_cast<TouchType>(&cube);

										}

								}

						}

						/* For each cube in the table of cubes */

						cube_iter = cube_table.begin();

						for(;cube_iter != cube_table.end(); ++cube_iter)

						{

								/* For each of the six faces, if the face is used */

								Cube& cube = *cube_iter;

								for(int i=0;i<6;++i)

										if(cube.faces[i].used)

										{

												CubeFace& face = cube.faces[i];  // The face

												// Let e be one of the four edges of the face

												for(int e=0;e<4;++e)

												{

														Cube* cube_n;  // Neighbouring cube

														int i_n, e_n;

														// Do the actual gluing ...

														if(!connect) 

																find_neighbour_face_disconnect(cube, i, e, cube_n, i_n, e_n);

														else

																find_neighbour_face_connect(cube, i, e, cube_n, i_n, e_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. */

														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);

												}

										}

						}

						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;

						}

				}

				template<class T>

				void CubeGrid<T>::cuberille_mesh(const RGrid<T>& voxel_grid, Manifold& m, bool connect)

				{

						typename list<Cube>::iterator cube_iter = cube_table.begin();

						/* For each cube in the table of cubes */

						cube_iter = cube_table.begin();

						for(;cube_iter != cube_table.end(); ++cube_iter)

						{

								/* For each of the six faces, if the face is used */

								Cube& cube = *cube_iter;

								for(int i=0;i<6;++i)

										if(cube.faces[i].used)

										{

												CubeFace& face = cube.faces[i];  // The face

												FaceIter hmesh_face = m.create_face();

												// Let e be one of the four edges of the face

												for(int e=0;e<4;++e)

												{

														Cube* cube_n;  // Neighbouring cube

														int i_n, e_n;

														// Do the actual gluing ...

														if(!connect) 

																find_neighbour_face_disconnect(cube, i, e, cube_n, i_n, e_n);

														else

																find_neighbour_face_connect(cube, i, e, cube_n, i_n, e_n);

														CubeFace& face_n = cube_n->faces[i_n];

														HalfEdgeIter he = face.get_he(m,e);

														HalfEdgeIter he_n = face_n.get_he(m,e_n);

														link(he,face.get_he(m,(e+1)%4));

														link(face.get_he(m,(e+3)%4),he);

														glue(he, he_n);

														he->face = hmesh_face;

												}

												hmesh_face->last = face.get_he(m,0);

										}

						}

						cube_iter = cube_table.begin();

						for(;cube_iter != cube_table.end(); ++cube_iter)

						{

								/* For each of the six faces, if the face is used */

								Cube& cube = *cube_iter;

								for(int i=0;i<6;++i)

										if(cube.faces[i].used)

										{

												CubeFace& face = cube.faces[i];  // The face

												// Let e be one of the four edges of the face

												for(int e=0;e<4;++e)

														if(face.he[e]->vert == NULL_VERTEX_ITER)

														{

																Vec3i to_face = N6i[i];

																Vec3i along_edge = edges[i][e];	

																Vec3i to_edge = cross(along_edge,to_face);

																Vec3f pos = Vec3f(cube.pos) + 0.5*Vec3f(to_face + along_edge + to_edge);

																HalfEdgeIter last = face.he[e];

																VertexIter vi = m.create_vertex(pos);

																vi->he = last->opp;

																HalfEdgeIter he = last;

																do

																{

																		he->vert = vi;

																		he = he->next->opp;

																}

																while(he != last);

														}

										}

						}

				}

				/* 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 triangulate(const RGrid<T>& voxel_grid, Manifold& m, 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(size_t 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*/

				void remove_duplicates(Manifold& m)

				{

						// 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(he))

												{

														++vi;

														m.collapse_halfedge(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 push_vertices(const RGrid<T>& voxel_grid, Manifold& m, float iso)

				{

						GradientFilter<RGrid<T> > grad(&voxel_grid);

						TrilinFilter<GradientFilter<RGrid<T> > > ginter(&grad);

						TrilinFilter<RGrid<T> > inter(&voxel_grid);

						for(int i=0;i<4;++i)

						{

								for(VertexIter vi= m.vertices_begin();

										vi != m.vertices_end(); ++vi)

								{

										Vec3f p = vi->pos;

										if(ginter.in_domain(p))

										{

												Vec3f g = ginter(p);

												float s = sqr_length(g);

												if(s>0.01)

												{

														float v = inter(p)-iso;

														vi->pos = (p-g*v/s);

												}

										}

								}

						}

				}

		}

		template<typename T>

		void fair_polygonize(const RGrid<T>& voxel_grid, 

												 Manifold& mani, 

												 float iso,

												 bool invert_contour, 

												 bool connect,

												 bool push)

		{

				CubeGrid<T> cube_grid(voxel_grid, iso, invert_contour);

				cube_grid.dual_cuberille_mesh(voxel_grid, mani, connect);

 				triangulate(voxel_grid, mani, iso);

 				remove_duplicates(mani);

				if(push) push_vertices(voxel_grid, mani, iso);

		}

		template<typename T>

		void cuberille_polygonize(const RGrid<T>& voxel_grid, 

															Manifold& mani, 

															float iso,

															bool invert_contour, 

															bool connect,

															bool push)

		{

				CubeGrid<T> cube_grid(voxel_grid, iso, invert_contour);

				cube_grid.cuberille_mesh(voxel_grid, mani, connect);

				if(push) push_vertices(voxel_grid, mani, iso);

		}

		template void fair_polygonize<unsigned char>(const RGrid<unsigned char>&,

																								 Manifold&, float, bool,bool,bool);

		template void fair_polygonize<float>(const RGrid<float>&,

																				 Manifold&, float, bool,bool,bool);

		template void cuberille_polygonize<unsigned char>(const RGrid<unsigned char>&,

																											Manifold&, float, bool,bool,bool);

		template void cuberille_polygonize<float>(const RGrid<float>&,

																							Manifold&, float, bool,bool,bool);

}