//------------------------------------------------------------------------------
template<class T> inline
int DualMC<T>::getCellCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso) const {
// determine for each cube corner if it is outside or inside
int code = 0;
if(data[gA(cx,cy,cz)] >= iso)
code |= 1;
if(data[gA(cx+1,cy,cz)] >= iso)
code |= 2;
if(data[gA(cx,cy+1,cz)] >= iso)
code |= 4;
if(data[gA(cx+1,cy+1,cz)] >= iso)
code |= 8;
if(data[gA(cx,cy,cz+1)] >= iso)
code |= 16;
if(data[gA(cx+1,cy,cz+1)] >= iso)
code |= 32;
if(data[gA(cx,cy+1,cz+1)] >= iso)
code |= 64;
if(data[gA(cx+1,cy+1,cz+1)] >= iso)
code |= 128;
return code;
}
//------------------------------------------------------------------------------
template<class T> inline
int DualMC<T>::getDualPointCode(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso, DMCEdgeCode const edge) const {
int cubeCode = getCellCode(cx, cy, cz, iso);
// is manifold dual marching cubes desired?
if(generateManifold) {
// The Manifold Dual Marching Cubes approach from Rephael Wenger as described in
// chapter 3.3.5 of his book "Isosurfaces: Geometry, Topology, and Algorithms"
// is implemente here.
// If a problematic C16 or C19 configuration shares the ambiguous face
// with another C16 or C19 configuration we simply invert the cube code
// before looking up dual points. Doing this for these pairs ensures
// manifold meshes.
// But this removes the dualism to marching cubes.
// check if we have a potentially problematic configuration
uint8_t const direction = problematicConfigs[uint8_t(cubeCode)];
// If the direction code is in {0,...,5} we have a C16 or C19 configuration.
if(direction != 255) {
// We have to check the neighboring cube, which shares the ambiguous
// face. For this we decode the direction. This could also be done
// with another lookup table.
// copy current cube coordinates into an array.
int32_t neighborCoords[] = {cx,cy,cz};
// get the dimension of the non-zero coordinate axis
unsigned int const component = direction >> 1;
// get the sign of the direction
int32_t delta = (direction & 1) == 1 ? 1 : -1;
// modify the correspong cube coordinate
neighborCoords[component] += delta;
// have we left the volume in this direction?
if(neighborCoords[component] >= 0 && neighborCoords[component] < (dims[component]-1)) {
// get the cube configuration of the relevant neighbor
int neighborCubeCode = getCellCode(neighborCoords[0], neighborCoords[1], neighborCoords[2], iso);
// Look up the neighbor configuration ambiguous face direction.
// If the direction is valid we have a C16 or C19 neighbor.
// As C16 and C19 have exactly one ambiguous face this face is
// guaranteed to be shared for the pair.
if(problematicConfigs[uint8_t(neighborCubeCode)] != 255) {
// replace the cube configuration with its inverse.
cubeCode ^= 0xff;
}
}
}
}
for(int i = 0; i < 4; ++i)
if(dualPointsList[cubeCode][i] & edge) {
return dualPointsList[cubeCode][i];
}
return 0;
}
//------------------------------------------------------------------------------
template<class T> inline
void DualMC<T>::calculateDualPoint(int32_t const cx, int32_t const cy, int32_t const cz, VolumeDataType const iso, int const pointCode, Vertex & v) const {
// initialize the point with lower voxel coordinates
v.x = cx;
v.y = cy;
v.z = cz;
// compute the dual point as the mean of the face vertices belonging to the
// original marching cubes face
Vertex p;
p.x=0;
p.y=0;
p.z=0;
int points = 0;
// sum edge intersection vertices using the point code
if(pointCode & EDGE0) {
p.x += ((float)iso - (float)data[gA(cx,cy,cz)])/((float)data[gA(cx+1,cy,cz)]-(float)data[gA(cx,cy,cz)]);
points++;
}
if(pointCode & EDGE1) {
p.x += 1.0f;
p.z += ((float)iso - (float)data[gA(cx+1,cy,cz)])/((float)data[gA(cx+1,cy,cz+1)]-(float)data[gA(cx+1,cy,cz)]);
points++;
}
if(pointCode & EDGE2) {
p.x += ((float)iso - (float)data[gA(cx,cy,cz+1)])/((float)data[gA(cx+1,cy,cz+1)]-(float)data[gA(cx,cy,cz+1)]);
p.z += 1.0f;
points++;
}
if(pointCode & EDGE3) {
p.z += ((float)iso - (float)data[gA(cx,cy,cz)])/((float)data[gA(cx,cy,cz+1)]-(float)data[gA(cx,cy,cz)]);
points++;
}
if(pointCode & EDGE4) {
p.x += ((float)iso - (float)data[gA(cx,cy+1,cz)])/((float)data[gA(cx+1,cy+1,cz)]-(float)data[gA(cx,cy+1,cz)]);
p.y += 1.0f;
points++;
}
if(pointCode & EDGE5) {
p.x += 1.0f;
p.z += ((float)iso - (float)data[gA(cx+1,cy+1,cz)])/((float)data[gA(cx+1,cy+1,cz+1)]-(float)data[gA(cx+1,cy+1,cz)]);
p.y += 1.0f;
points++;
}
if(pointCode & EDGE6) {
p.x += ((float)iso - (float)data[gA(cx,cy+1,cz+1)])/((float)data[gA(cx+1,cy+1,cz+1)]-(float)data[gA(cx,cy+1,cz+1)]);
p.z += 1.0f;
p.y += 1.0f;
points++;
}
if(pointCode & EDGE7) {
p.z += ((float)iso - (float)data[gA(cx,cy+1,cz)])/((float)data[gA(cx,cy+1,cz+1)]-(float)data[gA(cx,cy+1,cz)]);
p.y += 1.0f;
points++;
}
if(pointCode & EDGE8) {
p.y += ((float)iso - (float)data[gA(cx,cy,cz)])/((float)data[gA(cx,cy+1,cz)]-(float)data[gA(cx,cy,cz)]);
points++;
}
if(pointCode & EDGE9) {
p.x += 1.0f;
p.y += ((float)iso - (float)data[gA(cx+1,cy,cz)])/((float)data[gA(cx+1,cy+1,cz)]-(float)data[gA(cx+1,cy,cz)]);
points++;
}
if(pointCode & EDGE10) {
p.x += 1.0f;
p.y += ((float)iso - (float)data[gA(cx+1,cy,cz+1)])/((float)data[gA(cx+1,cy+1,cz+1)]-(float)data[gA(cx+1,cy,cz+1)]);
p.z += 1.0f;
points++;
}
if(pointCode & EDGE11) {
p.z += 1.0f;
p.y += ((float)iso - (float)data[gA(cx,cy,cz+1)])/((float)data[gA(cx,cy+1,cz+1)]-(float)data[gA(cx,cy,cz+1)]);
points++;
}
// divide by number of accumulated points
float invPoints = 1.0f / (float)points;
p.x*= invPoints;
p.y*= invPoints;
p.z*= invPoints;
// offset point by voxel coordinates
v.x += p.x;
v.y += p.y;
v.z += p.z;
}
//------------------------------------------------------------------------------
template<class T> inline
QuadIndexType DualMC<T>::getSharedDualPointIndex(
int32_t const cx, int32_t const cy, int32_t const cz,
VolumeDataType const iso, DMCEdgeCode const edge,
std::vector<Vertex> & vertices
) {
// create a key for the dual point from its linearized cell ID and point code
DualPointKey key;
key.linearizedCellID = gA(cx,cy,cz);
key.pointCode = getDualPointCode(cx,cy,cz,iso,edge);
// have we already computed the dual point?
auto iterator = pointToIndex.find(key);
if(iterator != pointToIndex.end()) {
// just return the dual point index
return iterator->second;
} else {
// create new vertex and vertex id
QuadIndexType newVertexId = vertices.size();
vertices.emplace_back();
calculateDualPoint(cx,cy,cz,iso,key.pointCode, vertices.back());
// insert vertex ID into map and also return it
pointToIndex[key] = newVertexId;
return newVertexId;
}
}
//------------------------------------------------------------------------------
template<class T> inline
void DualMC<T>::build(
VolumeDataType const * data,
int32_t const dimX, int32_t const dimY, int32_t const dimZ,
VolumeDataType const iso,
bool const generateManifold,
bool const generateSoup,
std::vector<Vertex> & vertices,
std::vector<Quad> & quads
) {
// set members
this->dims[0] = dimX;
this->dims[1] = dimY;
this->dims[2] = dimZ;
this->data = data;
this->generateManifold = generateManifold;
// clear vertices and quad indices
vertices.clear();
quads.clear();
// generate quad soup or shared vertices quad list
if(generateSoup) {
buildQuadSoup(iso,vertices,quads);
} else {
buildSharedVerticesQuads(iso,vertices,quads);
}
}
//------------------------------------------------------------------------------
template<class T> inline
void DualMC<T>::buildQuadSoup(
VolumeDataType const iso,
std::vector<Vertex> & vertices,
std::vector<Quad> & quads
) {
int32_t const reducedX = dims[0] - 2;
int32_t const reducedY = dims[1] - 2;
int32_t const reducedZ = dims[2] - 2;
Vertex vertex0;
Vertex vertex1;
Vertex vertex2;
Vertex vertex3;
int pointCode;
// iterate voxels
for(int32_t z = 0; z < reducedZ; ++z)
for(int32_t y = 0; y < reducedY; ++y)
for(int32_t x = 0; x < reducedX; ++x) {
// construct quad for x edge
if(z > 0 && y > 0) {
// is edge intersected?
bool const entering = data[gA(x,y,z)] < iso && data[gA(x+1,y,z)] >= iso;
bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x+1,y,z)] < iso;
if(entering || exiting){
// generate quad
pointCode = getDualPointCode(x,y,z,iso,EDGE0);
calculateDualPoint(x,y,z,iso,pointCode, vertex0);
pointCode = getDualPointCode(x,y,z-1,iso,EDGE2);
calculateDualPoint(x,y,z-1,iso,pointCode, vertex1);
pointCode = getDualPointCode(x,y-1,z-1,iso,EDGE6);
calculateDualPoint(x,y-1,z-1,iso,pointCode, vertex2);
pointCode = getDualPointCode(x,y-1,z,iso,EDGE4);
calculateDualPoint(x,y-1,z,iso,pointCode, vertex3);
if(entering) {
vertices.emplace_back(vertex0);
vertices.emplace_back(vertex1);
vertices.emplace_back(vertex2);
vertices.emplace_back(vertex3);
} else {
vertices.emplace_back(vertex0);
vertices.emplace_back(vertex3);
vertices.emplace_back(vertex2);
vertices.emplace_back(vertex1);
}
}
}
// construct quad for y edge
if(z > 0 && x > 0) {
// is edge intersected?
bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y+1,z)] >= iso;
bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y+1,z)] < iso;
if(entering || exiting){
// generate quad
pointCode = getDualPointCode(x,y,z,iso,EDGE8);
calculateDualPoint(x,y,z,iso,pointCode, vertex0);
pointCode = getDualPointCode(x,y,z-1,iso,EDGE11);
calculateDualPoint(x,y,z-1,iso,pointCode, vertex1);
pointCode = getDualPointCode(x-1,y,z-1,iso,EDGE10);
calculateDualPoint(x-1,y,z-1,iso,pointCode, vertex2);
pointCode = getDualPointCode(x-1,y,z,iso,EDGE9);
calculateDualPoint(x-1,y,z,iso,pointCode, vertex3);
if(exiting) {
vertices.emplace_back(vertex0);
vertices.emplace_back(vertex1);
vertices.emplace_back(vertex2);
vertices.emplace_back(vertex3);
} else {
vertices.emplace_back(vertex0);
vertices.emplace_back(vertex3);
vertices.emplace_back(vertex2);
vertices.emplace_back(vertex1);
}
}
}
// construct quad for z edge
if(x > 0 && y > 0) {
// is edge intersected?
bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y,z+1)] >= iso;
bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y,z+1)] < iso;
if(entering || exiting){
// generate quad
pointCode = getDualPointCode(x,y,z,iso,EDGE3);
calculateDualPoint(x,y,z,iso,pointCode, vertex0);
pointCode = getDualPointCode(x-1,y,z,iso,EDGE1);
calculateDualPoint(x-1,y,z,iso,pointCode, vertex1);
pointCode = getDualPointCode(x-1,y-1,z,iso,EDGE5);
calculateDualPoint(x-1,y-1,z,iso,pointCode, vertex2);
pointCode = getDualPointCode(x,y-1,z,iso,EDGE7);
calculateDualPoint(x,y-1,z,iso,pointCode, vertex3);
if(exiting) {
vertices.emplace_back(vertex0);
vertices.emplace_back(vertex1);
vertices.emplace_back(vertex2);
vertices.emplace_back(vertex3);
} else {
vertices.emplace_back(vertex0);
vertices.emplace_back(vertex3);
vertices.emplace_back(vertex2);
vertices.emplace_back(vertex1);
}
}
}
}
// generate triangle soup quads
size_t const numQuads = vertices.size() / 4;
quads.reserve(numQuads);
for (size_t i = 0; i < numQuads; ++i) {
quads.emplace_back(i * 4, i * 4 + 1, i * 4 + 2, i * 4 + 3);
}
}
//------------------------------------------------------------------------------
template<class T> inline
void DualMC<T>::buildSharedVerticesQuads(
VolumeDataType const iso,
std::vector<Vertex> & vertices,
std::vector<Quad> & quads
) {
int32_t const reducedX = dims[0] - 2;
int32_t const reducedY = dims[1] - 2;
int32_t const reducedZ = dims[2] - 2;
QuadIndexType i0,i1,i2,i3;
pointToIndex.clear();
// iterate voxels
for(int32_t z = 0; z < reducedZ; ++z)
for(int32_t y = 0; y < reducedY; ++y)
for(int32_t x = 0; x < reducedX; ++x) {
// construct quads for x edge
if(z > 0 && y > 0) {
bool const entering = data[gA(x,y,z)] < iso && data[gA(x+1,y,z)] >= iso;
bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x+1,y,z)] < iso;
if(entering || exiting){
// generate quad
i0 = getSharedDualPointIndex(x,y,z,iso,EDGE0,vertices);
i1 = getSharedDualPointIndex(x,y,z-1,iso,EDGE2,vertices);
i2 = getSharedDualPointIndex(x,y-1,z-1,iso,EDGE6,vertices);
i3 = getSharedDualPointIndex(x,y-1,z,iso,EDGE4,vertices);
if(entering) {
quads.emplace_back(i0,i1,i2,i3);
} else {
quads.emplace_back(i0,i3,i2,i1);
}
}
}
// construct quads for y edge
if(z > 0 && x > 0) {
bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y+1,z)] >= iso;
bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y+1,z)] < iso;
if(entering || exiting){
// generate quad
i0 = getSharedDualPointIndex(x,y,z,iso,EDGE8,vertices);
i1 = getSharedDualPointIndex(x,y,z-1,iso,EDGE11,vertices);
i2 = getSharedDualPointIndex(x-1,y,z-1,iso,EDGE10,vertices);
i3 = getSharedDualPointIndex(x-1,y,z,iso,EDGE9,vertices);
if(exiting) {
quads.emplace_back(i0,i1,i2,i3);
} else {
quads.emplace_back(i0,i3,i2,i1);
}
}
}
// construct quads for z edge
if(x > 0 && y > 0) {
bool const entering = data[gA(x,y,z)] < iso && data[gA(x,y,z+1)] >= iso;
bool const exiting = data[gA(x,y,z)] >= iso && data[gA(x,y,z+1)] < iso;
if(entering || exiting){
// generate quad
i0 = getSharedDualPointIndex(x,y,z,iso,EDGE3,vertices);
i1 = getSharedDualPointIndex(x-1,y,z,iso,EDGE1,vertices);
i2 = getSharedDualPointIndex(x-1,y-1,z,iso,EDGE5,vertices);
i3 = getSharedDualPointIndex(x,y-1,z,iso,EDGE7,vertices);
if(exiting) {
quads.emplace_back(i0,i1,i2,i3);
} else {
quads.emplace_back(i0,i3,i2,i1);
}
}
}
}
}
