Rev 61 | 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 "HMeshUtil/triangulate.h"
#include "HMeshUtil/build_manifold.h"
#include "HMeshUtil/x3d_save.h"
#include "HMeshUtil/obj_save.h"
#include "HMeshUtil/mesh_optimization.h"
#include "HMeshUtil/smooth.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 HMeshUtil;
using namespace CMP;
using namespace CGLA;
using namespace std;
namespace
{
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)}};
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;
}
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;
}
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;
}
struct Cube
{
Vec3i pos;
CubeFace faces[6];
Cube(const Vec3i& _pos): pos(_pos) {}
};
template<class T>
class Mesher
{
const RGrid<T>* const voxel_grid;
Manifold* m;
HashTable<HashKey3usi,Cube*> cube_hash;
list<Cube> cube_table;
std::vector<VertexIter> mul_vertices;
public:
Mesher(const RGrid<T>* _voxel_grid,
Manifold* _m,
float iso, bool inv_cont);
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)
{
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(comp(inv_cont,val0,iso))
{
bool interface_cell = false;
vector<VertexIter> vertices;
Cube* cube;
for(int l=0;l<6;++l)
{
Vec3i pi1(pi0);
pi1 += N6i[l];
bool indom = voxel_grid->in_domain(pi1);
if(!indom ||
!comp(inv_cont,float((*voxel_grid)[pi1]),iso))
{
float val1 = iso;
if(indom)
val1 = (*voxel_grid)[pi1];
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();
}
int face_idx = dir_to_face(N6i[l]);
CubeFace& face = cube->faces[face_idx];
face.used = true;
float t= (iso-val0)/(val1-val0);
Vec3f p = Vec3f(pi0)*(1.0f-t)+Vec3f(pi1)*t;
face.vert = m->create_vertex(p);
// face.vert = m->create_vertex(Vec3f(pi0)+
// N6f[l]/2.0f);
face.vert->touched = reinterpret_cast<int>(cube);
}
}
}
}
}
template<class T>
void Mesher<T>::process_cubes()
{
typename list<Cube>::iterator cube_iter = cube_table.begin();
for(;cube_iter != cube_table.end(); ++cube_iter)
{
Cube& cube = *cube_iter;
for(int i=0;i<6;++i)
if(cube.faces[i].used)
{
Vec3i pos = cube.pos;
Vec3i dir = N6i[i];
for(int e=0;e<4;++e)
{
Vec3i edge = edges[i][e];
Vec3i to_edge = cross(edge,dir);
CubeFace& face = cube.faces[i];
Cube* cube_n;
Vec3i dir_n;
Vec3i edge_n;
// 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;
// }
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;
}
int i_n = dir_to_face(dir_n);
CubeFace& face_n = cube_n->faces[i_n];
int e_n = dir_to_edge(i_n, edge_n);
assert(face_n.used);
HalfEdgeIter& he = face.get_he(m,e);
face.vert->he = he;
he->vert = face_n.vert;
HalfEdgeIter& he_n = face_n.get_he(m, (e_n+3)%4);
link(he, he_n);
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()
{
HalfEdgeIter he0 = m->halfedges_begin();
while(he0 != m->halfedges_end())
{
if(he0->face == NULL_FACE_ITER)
{
FaceIter fi = m->create_face();
fi->last = he0;
assert(he0->face == NULL_FACE_ITER);
assert(he0->prev != NULL_HALFEDGE_ITER);
assert(he0->next != NULL_HALFEDGE_ITER);
assert(he0->vert != NULL_VERTEX_ITER);
assert(he0->vert->he != NULL_HALFEDGE_ITER);
assert(he0->opp != NULL_HALFEDGE_ITER);
he0->face = fi;
HalfEdgeIter he = he0->next;
while(he != he0)
{
assert(he->face == NULL_FACE_ITER);
assert(he->prev != NULL_HALFEDGE_ITER);
assert(he->next != NULL_HALFEDGE_ITER);
assert(he->vert != NULL_VERTEX_ITER);
assert(he->vert->he != NULL_HALFEDGE_ITER);
assert(he->opp != NULL_HALFEDGE_ITER);
he->face = fi;
he = he->next;
}
}
++he0;
}
}
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
{
work = 0;
for(FaceIter fi =m->faces_begin(); fi != m->faces_end(); ++fi)
{
vector<VertexIter> verts;
for(FaceCirculator fc(fi); !fc.end(); ++fc)
verts.push_back(fc.get_vertex());
if(verts.size() == 3) continue;
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() && N > 3)
{
cout << __FILE__ << __LINE__ << " " << N << endl;
for(int i=0;i<N-2;++i)
for(int j=i+2;j<N; ++j)
{
cout << "i=" << i << " j=" << j << " ";
cout << &*verts[i] << " " << &*verts[j] << " ";
cout << is_connected(verts[i], verts[j]) << endl;
}
}
if(!vpairs.empty())
{
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);
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<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
}
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);
}