# Check if neighbors are active in spatially subdivided voxelization

I have a voxelization that I need to write out to a binary STL file. These voxels are inside a higher-level 3D grid, with each cell in the grid being made up of voxels. Individual cells may or may not have any active voxels in them.

I've optimized the file size significantly by not writing out triangles for interior faces - faces between active voxels.

However, the method I used is very ugly, and I think could be improved. I'm not sure if there's an applicable design pattern or just a better way of thinking about it.

• It's ugly
• It's repetitive

Regardless, I want to learn something from this code and see how I can do better in the future.

Note spacing might be a bit odd, using clang-format and an 80-character limit.

writeVoxelSTL():

void writeVoxelAsSTL(
std::ofstream& fileStream,
unsigned& numTris,
const unsigned xCell,
const unsigned yCell,
const unsigned zCell,
const unsigned xVoxel,
const unsigned yVoxel,
const unsigned zVoxel) const
{
// Now need to check if each of its sides need to be written out
bool left = isNeighborActive(
xCell, yCell, zCell, static_cast<int>(xVoxel) - 1, yVoxel, zVoxel);
bool right =
isNeighborActive(xCell, yCell, zCell, xVoxel + 1, yVoxel, zVoxel);
bool bottom = isNeighborActive(
xCell, yCell, zCell, xVoxel, static_cast<int>(yVoxel) - 1, zVoxel);
bool top = isNeighborActive(xCell, yCell, zCell, xVoxel, yVoxel + 1, zVoxel);
bool front = isNeighborActive(
xCell, yCell, zCell, xVoxel, yVoxel, static_cast<int>(zVoxel) - 1);
bool back = isNeighborActive(xCell, yCell, zCell, xVoxel, yVoxel, zVoxel + 1);

Point minPoint = highLevelGridMinpoint;

// Center of the voxel
Point center;
center.x = minPoint.x + (((xCell * cellDim.x) + xVoxel) * voxelWidth);
center.y = minPoint.y + (((yCell * cellDim.y) + yVoxel) * voxelWidth);
center.z = minPoint.z + (((zCell * cellDim.z) + zVoxel) * voxelWidth);

// From this center, we can get the 8 corners of the voxel and thus the 12
// triangles

// They are defined like this:
//              7-------6
//             /|      /|
//            4-+-----5 |
//            | |     | |   y
//            | 3-----+-2   | z
//            |/      |/    |/
//            0-------1     +--x

Point p0(center.x - halfVoxel, center.y - halfVoxel, center.z - halfVoxel);
Point p1(center.x + halfVoxel, center.y - halfVoxel, center.z - halfVoxel);
Point p2(center.x + halfVoxel, center.y - halfVoxel, center.z + halfVoxel);
Point p3(center.x - halfVoxel, center.y - halfVoxel, center.z + halfVoxel);
Point p4(center.x - halfVoxel, center.y + halfVoxel, center.z - halfVoxel);
Point p5(center.x + halfVoxel, center.y + halfVoxel, center.z - halfVoxel);
Point p6(center.x + halfVoxel, center.y + halfVoxel, center.z + halfVoxel);
Point p7(center.x - halfVoxel, center.y + halfVoxel, center.z + halfVoxel);

// Now for the 12 triangles made of these points

// Left side
if (!left)
{
writeTriangleSTL(fileStream, p0, p3, p4);
writeTriangleSTL(fileStream, p3, p4, p7);
numTris += 2;
}

// Right Side
if (!right)
{
writeTriangleSTL(fileStream, p1, p2, p5);
writeTriangleSTL(fileStream, p2, p5, p6);
numTris += 2;
}

// Bottom Side
if (!bottom)
{
writeTriangleSTL(fileStream, p0, p1, p2);
writeTriangleSTL(fileStream, p0, p2, p3);
numTris += 2;
}

// Top Side
if (!top)
{
writeTriangleSTL(fileStream, p4, p5, p6);
writeTriangleSTL(fileStream, p4, p6, p7);
numTris += 2;
}

// Front side
if (!front)
{
writeTriangleSTL(fileStream, p0, p1, p4);
writeTriangleSTL(fileStream, p1, p4, p5);
numTris += 2;
}

// Back side
if (!back)
{
writeTriangleSTL(fileStream, p3, p2, p7);
writeTriangleSTL(fileStream, p2, p7, p6);
numTris += 2;
}
}


isNeighborActive():

Note: cellDim is the number of voxels in each direction for a cell. gridDim is the number of cells in each direction.

bool isNeighborActive(
long xCell, long yCell, long zCell, long xVoxel, long yVoxel, long zVoxel)
const
{
// Get the "real" coordinates
if (xVoxel < 0)
{
xCell--;
xVoxel = cellDim.x - 1;
if (xCell < 0)
{
return false;
}
}
else if (xVoxel > (cellDim.x - 1))
{
xCell++;
xVoxel = 0;
if (xCell > (gridDim.x - 1))
{
return false;
}
}

if (yVoxel < 0)
{
yCell--;
yVoxel = cellDim.y - 1;
if (yCell < 0)
{
return false;
}
}
else if (yVoxel > (cellDim.y - 1))
{
yCell++;
yVoxel = 0;
if (yCell > (gridDim.y - 1))
{
return false;
}
}

if (zVoxel < 0)
{
zCell--;
zVoxel = cellDim.z - 1;
if (zCell < 0)
{
return false;
}
}
else if (zVoxel > (cellDim.z - 1))
{
zCell++;
zVoxel = 0;
if (zCell > (gridDim.z - 1))
{
return false;
}
}

const cellVoxelization& cell = modelVoxelization.at(xCell, yCell, zCell);

if (cell.getSize() == 0)
{
return false;
}

bool active = cell.at(xVoxel, yVoxel, zVoxel).active;

return active;
}


I considered using some kind of enum for the directions, but then that just shifts the ugliness down to the next method which has have a switch, sadly.

writeVoxelAsSTL takes unsigned arguments, but it seems like the native type is actually a long. Could it take longs instead? For example, unsigned(2^31) looks like it would not be a valid value since it would overflow into -1 upon conversion.

You could take an ostream& argument instead so the code isn't coupled to files.

You could factor out the bools and directions into arrays, and rewrite much of the code into loops.

Further improvement would depend on how this is being used. How is the grid represented? Can you combine neighboring voxel faces? For example, would a 100x1x1 stack require 804 triangles, or 12? Perhaps you could split this into multiple phases, one that extracts relevant planes, and one that writes planes to a stream. Maybe even one that first extracts a sparse list of relevant voxels depending on how you're visiting them.

If this needs to be performant, surface extraction from voxels is a very parallelizable algorithm. Consider just a 2x2x2 box: 38 out of 64 points computed are going to be shared. The center point will be computed 8 times but never used.

• The cells and voxel positions are unsigned for sure. I had to adjust isNeighbor to use some kind of signed argument so I would have a way of (for example) getting the previous neighbor of a voxel in the 0th position, which in this case would now be the -1th position in that axis. Jan 13 '20 at 19:08
• 2: I actually do want fstream here, this will only be used for writing out to a file. Jan 13 '20 at 19:08
• 4: Combining voxel faces would be ideal, but seems exceptionally challenging. We would only want to combine flush faces, and the only way to tell if the faces are flush would be to step outwards to the 26 voxels in each direction to determine their activation as well. IE, if my front face is active, AND my neighbors front face is active, but the neighbors FRONT neighbor isn't active, then we can make them flush. But what about the voxel to the right of my right neighbor, can I combine that as well? there may be an algorithm for this but I do not know it. Jan 13 '20 at 19:11
• One idea could be to collect only surface voxels. For example only voxels for which their +XYZ neighbor's occupancy state is different from their state, a total of 3 tests per voxel. It could be stored as a relative position and a bitflag indicating adjacency. Then each entry outputs [1-3] planes instead of [0-6], and you never need to check for duplicate planes. Regarding looping, I was imagining something like this: godbolt.org/z/PzkJDP
– butt
Jan 13 '20 at 19:40