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 "HMesh/build_manifold.h"

#include "Geometry/RGrid.h"

#include "Geometry/TrilinFilter.h"

#include "Geometry/GradientFilter.h"

#include "Geometry/Neighbours.h"

#include "Util/HashTable.h"

#include "Util/HashKey.h"

#include "fair_polygonize.h"

namespace HMesh

{

	using namespace Geometry;

	using namespace Util;

	using namespace CGLA;

	using namespace std;

	namespace 

	{

		void compute_dual(Manifold& m, const vector<Vec3f>& face_positions)

		{

			// make sure every face knows its number

			m.enumerate_faces();

			vector<Vec3f> vertices(m.no_faces());

			vector<int> faces;

			vector<int> indices;

			// Create new vertices. Each face becomes a vertex whose positio

			// is the centre of the face

			int i=0;

			for(FaceIter f=m.faces_begin(); f!=m.faces_end(); ++f,++i)

				vertices[i] = face_positions[f->touched];

			// Create new faces. Each vertex is a new face with N=valency of

			// edges.

			i=0;

			for(VertexIter v=m.vertices_begin(); v!= m.vertices_end(); ++v,++i)

				if(!is_boundary(v))

					{

						int no_edges = 0;

						vector<int> index_tmp;

						// Circulate around vertex

						VertexCirculator vc(v);

						for(; !vc.end(); ++vc)

							{

								FaceIter f = vc.get_face();

								if(f != NULL_FACE_ITER)

									{

										index_tmp.push_back(f->touched);

										++no_edges;

									}

							}

						// Push vertex indices for this face onto indices vector

						// The circulator moves around the face in a clockwise f

						// so we just reverse the ordering.

						indices.insert(indices.end(), index_tmp.rbegin(), index_tmp.rend());

						faces.push_back(no_edges);

					}

			// Clear the manifold before new geometry is inserted.

			m.clear();

			// And build

			build_manifold(m, vertices.size(), &vertices[0], faces.size(),

										 &faces[0],&indices[0]);

		}

		/*

			Vectors along edges. For each of six cube faces there are four such vectors

			since the faces of a cube are quads.

		*/

		const Vec3i edges[6][4] = {

			{Vec3i(0,-1,0),Vec3i(0,0,1),Vec3i(0,1,0),Vec3i(0,0,-1)},

			{Vec3i(0,1,0),Vec3i(0,0,1),Vec3i(0,-1,0),Vec3i(0,0,-1)},

			{Vec3i(0,0,-1),Vec3i(1,0,0),Vec3i(0,0,1),Vec3i(-1,0,0)},

			{Vec3i(0,0,1),Vec3i(1,0,0),Vec3i(0,0,-1),Vec3i(-1,0,0)},

			{Vec3i(-1,0,0),Vec3i(0,1,0),Vec3i(1,0,0),Vec3i(0,-1,0)},

			{Vec3i(1,0,0),Vec3i(0,1,0),Vec3i(-1,0,0),Vec3i(0,-1,0)}};

		/* Convert a direction vector to a cube face index. */

		int dir_to_face(const Vec3i& dir)

		{

			assert(dot(dir,dir)==1);

			if(dir[0] != 0)

				return (dir[0]<0) ? 0 : 1;

			if(dir[1] != 0)

				return (dir[1]<0) ? 2 : 3;

			if(dir[2] != 0)

				return (dir[2]<0) ? 4 : 5;

			return -1;

		}

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

		int dir_to_edge(int face, const Vec3i& edge)

		{

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

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

			assert(0);

			return -1;

		}

		Vec3f cube_corner(const Vec3i& pos, int i, int e)

		{

			Vec3i to_face = N6i[i];

			Vec3i along_edge = edges[i][e];	

			Vec3i to_edge = cross(along_edge,to_face);

			return Vec3f(pos) + 0.5*Vec3f(to_face+to_edge+along_edge);

		}

		/* Cube faces are created whenever a voxel that is inside is adjacent

			 to a voxel that is outside. Thus the voxel can be seen as a small

			 box that either contains material or not. A face is created where

			 full and empty boxes meet. A vertex is placed on every box face.

			 A halfedge is created for every edge of the (quadrilateral) face.

			 This halfedge connects the face to its neighbour. */

		struct CubeFace

		{

			bool used;

			HalfEdgeIter he[4];

			CubeFace();

			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)

		{

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

				he[i] = NULL_HALFEDGE_ITER;

		}

		/* A cube consists of six faces and a position. The cube centre is

			 an integer (voxel grid) position. */

		struct Cube

		{

			Vec3i pos;

			CubeFace faces[6];

			Cube(const Vec3i& _pos): pos(_pos) {}

		};

		/* The mesher class does the actual work. It is not pretty and it

			 is not seen by the user. */

		template<class T>

		class CubeGrid

		{

			float iso;

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

			list<Cube> cube_table;                   // Linked list containing

			// the actual cubes.

		public:

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

							 float iso);                   // iso value

			HalfEdgeIter edge_to_glue(const RGrid<T>& voxel_grid,

																Cube& cube, int face_idx, int edge_idx);

			void cuberille_mesh(const RGrid<T>&, Manifold& m, vector<Vec3f>& face_positions);

		};

		template<class T>

		HalfEdgeIter CubeGrid<T>::edge_to_glue(const RGrid<T>& voxel_grid,

																					 Cube& cube, int face_idx, int edge_idx)

		{

			/*  Figure.

			0|B

			-+-

			A|C

			Explanation: ... pending

			*/

			Vec3i dir = N6i[face_idx];

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

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

			Vec3i a_pos = cube.pos;

			Vec3i o_pos = a_pos + dir;

			Vec3i b_pos = o_pos + to_edge;

			Vec3i c_pos = a_pos + to_edge;

			if(!voxel_grid.in_domain(b_pos))

				return NULL_HALFEDGE_ITER;

			Cube* cube_n=0;

			Vec3i dir_n;   // direction from center to neighbour face

			Vec3i along_edge_n; // Edge vector for neighbour face

			if(!cube_hash.find_entry(HashKey3usi(b_pos),	cube_n))

				{

					if(cube.faces[dir_to_face(to_edge)].used)

						{

							dir_n = to_edge;

							along_edge_n = - along_edge; 

							cube_n = &cube;  

						}

					else

						{

							dir_n = dir;

							along_edge_n = - along_edge;

							cube_hash.find_entry(HashKey3usi(c_pos), cube_n);

							assert(cube_n->faces[dir_to_face(dir_n)].used);

						}									

				}

			else

				{

					float a_val = voxel_grid[a_pos];

					float b_val = voxel_grid[b_pos];

					float c_val = voxel_grid[c_pos];

					float o_val = voxel_grid[o_pos];

					bool connect = (a_val*b_val-c_val*o_val)/(a_val+b_val-c_val-o_val) >iso;

					if(connect || !cube.faces[dir_to_face(to_edge)].used)

						{

							dir_n = - to_edge;

							along_edge_n = - along_edge;

							assert(cube_n->faces[dir_to_face(dir_n)].used);

						}

					else

						{

							dir_n = to_edge;

							along_edge_n = - along_edge; 

							cube_n = &cube;  

						}

				}

			int i_n = dir_to_face(dir_n);

			int e_n = dir_to_edge(i_n, along_edge_n);

			return cube_n->faces[i_n].he[e_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): iso(_iso)

		{

			/* Loop over all voxels */

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

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

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

						{

							Vec3i pi0(i,j,k);

							float val0 = voxel_grid[pi0];

							/* If the voxel value is above the isovalue */

							if(val0>iso)

								{

									bool interface_cell = false;

									vector<VertexIter> vertices;

									Cube* cube;

									// For each neighbouring voxel.

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

										{

											// First make sure the voxel is inside the voxel grid.

											Vec3i pi1(pi0);

											pi1 += N6i[l];

											bool indom = voxel_grid.in_domain(pi1);

											/* If the neighbour is **not** in the grid or it 

												 is on the other side of the isovalue */

											if(indom && voxel_grid[pi1]<=iso)

												{

													/* Store the cube in a hash table. This code is

														 executed the first time we create a face

														 for the cube. In other words only cubes that 

														 actually have faces are stored in the hash */

													if(!interface_cell)

														{

															interface_cell = true;

															Cube** c_ptr;

															cube_hash.create_entry(HashKey3usi(pi0), c_ptr);

															Cube c(pi0);

															cube_table.push_back(c);

															cube = (*c_ptr) = &cube_table.back();

														}

													// Now add the face to the cube.

													int face_idx  = dir_to_face(N6i[l]);

													CubeFace& face = cube->faces[face_idx];

													face.used = true;

												}

										}

								}

						}		

		}

		template<class T>

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

																		 vector<Vec3f>& face_positions)

		{

			face_positions.clear();

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

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

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

			cube_iter = cube_table.begin();

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

				{

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

					Cube& cube = *cube_iter;

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

						if(cube.faces[i].used)

							{

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

								FaceIter hmesh_face = m.create_face();

								hmesh_face->touched = face_positions.size();

								Vec3i p0(cube.pos);

								Vec3i p1 = p0 + N6i[i];

								float v0 = voxel_grid[p0];

								float v1 = voxel_grid[p1];

								float t = (iso-v0)/(v1-v0);

								face_positions.push_back(Vec3f(p0)*(1.0f-t)+Vec3f(p1)*t);

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

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

									{

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

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

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

										he->face = hmesh_face;

										HalfEdgeIter he_n = edge_to_glue(voxel_grid, cube, i, e);

										if(he_n != NULL_HALFEDGE_ITER)

											glue(he, he_n);

									}

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

							}

				}

			for(HalfEdgeIter h = m.halfedges_begin(); h != m.halfedges_end(); ++h)

				{

					if (h->opp == NULL_HALFEDGE_ITER)

						{

							HalfEdgeIter ho = m.create_halfedge();

							glue(ho,h);

						}

				}

			cube_iter = cube_table.begin();

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

				{

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

					Cube& cube = *cube_iter;

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

						if(cube.faces[i].used)

							{

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

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

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

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

										{

											HalfEdgeIter last = face.he[e];

											HalfEdgeIter he = last;

											while(he->opp->face != NULL_FACE_ITER)

												{

													he = he->opp->prev;

													if(he == last) 

														break;

												}

											last = he;

											Vec3f pos = cube_corner(cube.pos, i, e);

											VertexIter vert = m.create_vertex(pos);

											do

												{

													he->vert = vert;

													he = he->next->opp;

												}

											while(he->face != NULL_FACE_ITER && he != last);

											if(he->face == NULL_FACE_ITER)

												{

													he->vert = vert;

													link(he,last->opp);

												}

											vert->he = last->opp;

										}

							}

				}

		}

	}

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

							if(vpairs.empty())

								{

									cout << "Warning: could not triangulate a face." 

											 << "Probably a vertex appears more than one time in other vertex's one-ring" << endl;

									continue;

								}

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

		}

	template<typename T>

	void fair_polygonize(const RGrid<T>& voxel_grid, 

											 Manifold& mani, 

											 float iso)

	{

		CubeGrid<T> cube_grid(voxel_grid, iso);

		vector<Vec3f> face_positions;

		cube_grid.cuberille_mesh(voxel_grid, mani, face_positions);

		cout << "Initial mesh valid : " << mani.is_valid() << endl;

		compute_dual(mani, face_positions);

	}

	template<typename T>

	void cuberille_polygonize(const RGrid<T>& voxel_grid, 

														Manifold& mani, 

														float iso,

														bool push)

	{

		CubeGrid<T> cube_grid(voxel_grid, iso);

		vector<Vec3f> face_positions;

		cube_grid.cuberille_mesh(voxel_grid, mani, face_positions);

		if(push)

			{

				for(VertexIter v=mani.vertices_begin(); v != mani.vertices_end(); ++v)

					{

						VertexCirculator vc(v);

						Vec3f p(0);

						int n=0;

						while(!vc.end())

							{

								if(vc.get_face() != NULL_FACE_ITER)

									{

										p += face_positions[vc.get_face()->touched];

										++n;

									}

								++vc;

							}

						v->pos = p/float(n);

					}

			}

	}

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

																							 Manifold&, float);

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

																			 Manifold&, float);

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

																										Manifold&, float, bool);

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

																						Manifold&, float, bool);

	template void triangulate(const RGrid<unsigned char>& voxel_grid, Manifold& m, float iso);

	template void triangulate(const RGrid<float>& voxel_grid, Manifold& m, float iso);

}