Rev 170 | Go to most recent revision | Blame | Last modification | View Log | RSS feed
#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);
}