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 "volume_polygonize.h"

namespace HMesh

{

	using namespace Geometry;

	using namespace Util;

	using namespace CGLA;

	using namespace std;

	namespace 

	{

		template<class T>

		Vec3f push_vertex(const RGrid<T>& grid, const Vec3f& pos0, float iso)

		{

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

			Vec3f corner_min = pos0 - Vec3f(1);

			Vec3f corner_max = pos0 + Vec3f(1);

			Vec3f pos = pos0;

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

				if(inter.in_domain(pos))

					{

						float val = inter(pos)-iso;

						Vec3f grad = inter.grad(pos);

						float sqlen_grad = sqr_length(grad);

						// Compute the projection in the gradient direction

						Vec3f d1 = - val * grad/sqlen_grad;

						Vec3f p_new = pos + 0.5 * d1; // "0.5 *" to add damping

						// Don't move outside the box.

						if(!(corner_max.all_g(p_new) && 

								 corner_min.all_l(p_new)))

							return pos;

						// Move so that the gradient points in the direction

						// of the original position

						Vec3f d2 = (pos0-pos) - 

							grad*dot((pos0-pos),grad)/sqlen_grad;

						// Compute new position and distance up to original position.

						p_new += 0.5 * d2; // "0.5 *" to add damping

						// Don't move outside the box.

						if(corner_max.all_g(p_new) && 

							 corner_min.all_l(p_new))

							pos = p_new;

						else

							return pos;

						// Did we converge?

						if(sqrt(sqr_length(d1)+sqr_length(d2))<.001)

							return pos;

					}

			return pos;

		}

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

		{

			// make sure every face knows its number

			m.enumerate_faces();

			int i=0;

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

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

					vertices[i] = m2.create_vertex(face_positions[f->touched]);

			vector<HalfEdgeIter> hedges(m.no_halfedges(), NULL_HALFEDGE_ITER);

			int he_num=0;

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

				if(!is_boundary(v))

				{

						FaceIter f = m2.create_face();

						vector<int> indices;

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

						{

								TouchType t = vc.get_halfedge()->touched;

								int h_idx = he_num++;

								vc.get_halfedge()->touched = h_idx;

								HalfEdgeIter h = m2.create_halfedge();

								h->face = f;

								h->touched = t;

								int f_idx = vc.get_face()->touched;

								h->vert = vertices[f_idx];

								hedges[h_idx] = h;

								indices.push_back(h_idx);

						}

						f->last = hedges[indices[0]];

						int N=indices.size();

						for(int j=0;j<N;++j)

						{

								HalfEdgeIter h1 = hedges[indices[j]];

								HalfEdgeIter h2 = hedges[indices[(j+1)%N]];

								link(h2,h1);

								h2->vert->he = h1;

						}

				}

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

				if(!is_boundary(v))

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

								if(!is_boundary(vc.get_vertex()))

								{

										int idx0 = vc.get_halfedge()->touched;

										int idx1 = vc.get_opp_halfedge()->touched;

										glue(hedges[idx0], hedges[idx1]);

								}

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

					if(h->opp == NULL_HALFEDGE_ITER)

							glue(h, m2.create_halfedge());

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

					check_boundary_consistency(v);

		}

		/*

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

			since the faces of a cube are quads.

		*/

		const Vec3i edges[6][4] = {

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

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

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

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

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

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

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

		int dir_to_face(const Vec3i& dir)

		{

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

			if(dir[0] != 0)

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

			if(dir[1] != 0)

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

			if(dir[2] != 0)

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

			return -1;

		}

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

		int dir_to_edge(int face, const Vec3i& edge)

		{

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

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

			assert(0);

			return -1;

		}

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

		{

			Vec3i to_face = N6i[i];

			Vec3i along_edge = edges[i][e];	

			Vec3i to_edge = cross(along_edge,to_face);

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

		}

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

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

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

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

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

			 This halfedge connects the face to its neighbour. */

		struct CubeFace

		{

				bool used;

				HalfEdgeIter he[4];

				CubeFace(): used(false)

						{

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

										he[i] = NULL_HALFEDGE_ITER;

						}

		};

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

			 an integer (voxel grid) position. */

		struct Cube

		{

			Vec3i pos;

			CubeFace faces[6];

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

		};

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

			 is not seen by the user. */

		template<class T>

		class CubeGrid

		{

			float iso;

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

			list<Cube> cube_table;                   // Linked list containing

			// the actual cubes.

		public:

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

							 float iso);                   // iso value

			HalfEdgeIter edge_to_glue(const RGrid<T>& voxel_grid,

																Cube& cube, int face_idx, int edge_idx);

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

													bool push_vertices,

													bool compute_face_positions,

													vector<Vec3f>& face_positions);

		};

		template<class T>

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

																					 Cube& cube, int face_idx, int edge_idx)

		{

			/*  Figure.

			0|B

			-+-

			A|C

			Explanation: ... pending

			*/

			Vec3i dir = N6i[face_idx];

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

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

			Vec3i a_pos = cube.pos;

			Vec3i o_pos = a_pos + dir;

			Vec3i b_pos = o_pos + to_edge;

			Vec3i c_pos = a_pos + to_edge;

			if(!voxel_grid.in_domain(b_pos))

				return NULL_HALFEDGE_ITER;

			Cube* cube_n=0;

			Vec3i dir_n;   // direction from center to neighbour face

			Vec3i along_edge_n; // Edge vector for neighbour face

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

				{

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

						{

							dir_n = to_edge;

							along_edge_n = - along_edge; 

							cube_n = &cube;  

						}

					else

						{

							dir_n = dir;

							along_edge_n = - along_edge;

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

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

						}									

				}

			else

				{

					float a_val = voxel_grid[a_pos];

					float b_val = voxel_grid[b_pos];

					float c_val = voxel_grid[c_pos];

					float o_val = voxel_grid[o_pos];

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

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

						{

							dir_n = - to_edge;

							along_edge_n = - along_edge;

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

						}

					else

						{

							dir_n = to_edge;

							along_edge_n = - along_edge; 

							cube_n = &cube;  

						}

				}

			int i_n = dir_to_face(dir_n);

			int e_n = dir_to_edge(i_n, along_edge_n);

			HalfEdgeIter he = cube_n->faces[i_n].he[e_n];

			if(cube_n == &cube)

					he->touched = 1;

			else

					he->touched = 0;

			return he;

		}

		template<class T>

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

		{

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

		}

		template<class T>

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

		{

			/* Loop over all voxels */

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

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

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

						{

							Vec3i pi0(i,j,k);

							float val0 = voxel_grid[pi0];

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

							if(val0>iso)

								{

									bool interface_cell = false;

									vector<VertexIter> vertices;

									Cube* cube;

									// For each neighbouring voxel.

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

										{

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

											Vec3i pi1(pi0);

											pi1 += N6i[l];

											bool indom = voxel_grid.in_domain(pi1);

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

												 is on the other side of the isovalue */

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

												{

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

														 executed the first time we create a face

														 for the cube. In other words only cubes that 

														 actually have faces are stored in the hash */

													if(!interface_cell)

														{

															interface_cell = true;

															Cube** c_ptr;

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

															Cube c(pi0);

															cube_table.push_back(c);

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

														}

													// Now add the face to the cube.

													int face_idx  = dir_to_face(N6i[l]);

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

													face.used = true;

												}

										}

								}

						}		

		}

		template<class T>

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

																		 bool push_vertices,

																		 bool compute_face_positions,

																		 vector<Vec3f>& face_positions)

		{

			face_positions.clear();

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

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

			cube_iter = cube_table.begin();

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

				{

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

					Cube& cube = *cube_iter;

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

						if(cube.faces[i].used)

							{

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

								FaceIter hmesh_face = m.create_face();

								hmesh_face->touched = face_positions.size();

								if(compute_face_positions)

								{

										Vec3i p0(cube.pos);

										Vec3i p1 = p0 + N6i[i];

										float v0 = voxel_grid[p0];

										float v1 = voxel_grid[p1];

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

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

								}

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

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

								{

										face.he[e] = m.create_halfedge();

										face.he[e]->face = hmesh_face;

								}

								hmesh_face->last = face.he[0];

							}

				}

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

			cube_iter = cube_table.begin();

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

				{

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

					Cube& cube = *cube_iter;

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

						if(cube.faces[i].used)

							{

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

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

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

									{

											HalfEdgeIter he = face.he[e];

											link(he,face.he[(e+1)%4]);

											link(face.he[(e+3)%4],he);

											if(he->opp == NULL_HALFEDGE_ITER)

											{

													HalfEdgeIter he_n = 

															edge_to_glue(voxel_grid, cube, i, e);

													if(he_n != NULL_HALFEDGE_ITER)

													{

															glue(he, he_n);

															he->touched = he_n->touched;

													}

											}

									}

							}

				}

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

				{

					if (h->opp == NULL_HALFEDGE_ITER)

						{

							HalfEdgeIter ho = m.create_halfedge();

							glue(ho,h);

						}

				}

			cube_iter = cube_table.begin();

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

				{

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

					Cube& cube = *cube_iter;

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

						if(cube.faces[i].used)

							{

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

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

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

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

										{

											HalfEdgeIter last = face.he[e];

											HalfEdgeIter he = last;

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

												{

													he = he->opp->prev;

													if(he == last) 

														break;

												}

											last = he;

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

											if(push_vertices)

													pos = push_vertex(voxel_grid, pos, iso);

											VertexIter vert = m.create_vertex(pos);

											do

												{

													he->vert = vert;

													he = he->next->opp;

												}

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

											if(he->face == NULL_FACE_ITER)

												{

													he->vert = vert;

													link(he,last->opp);

												}

											vert->he = last->opp;

										}

							}

				}

		}

	}

	template<typename T>

	void mc_polygonize(const RGrid<T>& voxel_grid, 

											 Manifold& mani_out, 

											 float iso)

	{

			Manifold mani;

			CubeGrid<T> cube_grid(voxel_grid, iso);

			vector<Vec3f> face_positions;

			cube_grid.cuberille_mesh(voxel_grid, mani, false, true, face_positions);

			compute_dual(mani, mani_out, face_positions);

	}

	template<typename T>

	void cuberille_polygonize(const RGrid<T>& voxel_grid, 

														Manifold& mani, 

														float iso,

														bool push)

	{

		CubeGrid<T> cube_grid(voxel_grid, iso);

		vector<Vec3f> face_positions;

		cube_grid.cuberille_mesh(voxel_grid, mani, push, false, face_positions);

	}

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

																							 Manifold&, float);

		template void mc_polygonize<float>(const RGrid<float>&, Manifold&, float);

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

																											Manifold&, float, bool);

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

																											Manifold&, float, bool);

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

																							Manifold&, float, bool);

}