Rev 90 | Blame | Compare with Previous | 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 "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.