Rev 191 | 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 "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);
}