0) The size of a chunk doesn't necessarily have to be a power-of-two (e.g. VoxelFarm uses 40^3 - still allows to use 16-bit indices). Speed gains associated with turning multiplication into shifts are negligible compared to choosing optimal chunk granularity for streaming/updates.
Since I didn't understand your question, I'll describe how my voxel engine currently works.
I populate chunks from Signed Distance Fields.
Each chunk contains raw voxel data (8-bit material IDs) and active (or intersecting the surface, zero-crossing) edges (Hermite data - see "Fast and Adaptive Polygon Conversion By Means Of Sparse Volumes"(2011)). An active edge exists only where a material transition takes place.
For each cell in a chunk (cell grid is dual to the voxel grid, i.e. cell corners are voxel centers and vice versa - for a nice pic, see page 7 in "A real-time virtual sculpting application by using an optimized hash-based octree"(2014)) I store up to 3 intersection points with normals on active edges (Hermite data).
As the cubic cell sweeps along the grid, only 3 edges need to be examined.
I calculate the voxel material at the current point (the maximal corner of the cell) and compare it to the previously computed neighbors (along axes X,Y and Z). It's better to explain with a picture, but I'll just dump the code:
// Calculates the so-called 'Hermite data':
// intersection points together with their normals
// (where the surface crosses with the edges of the dual grid).
// NOTE: as the cube sweeps along the grid, only 3 edges need to be examined.
const V3F voxelCenter = GetVoxelCenter( it.x, it.y, it.z ) + chunkOffset;
const SDF::Sample sample = _SDF->SampleAt( voxelCenter ); // sample the distance field
const MaterialID voxelMaterial = GetMaterialAt( sample, _settings );
volume->Set( it.x, it.y, it.z, voxelMaterial ); // store this voxel
// Examine previously computed neighbors and calculate intersections on zero-crossing edges.
if( volume->Get( it.x-1, it.y, it.z ) != voxelMaterial )
{
activeEdges.Add( ComputeIntersection( _SDF, voxelCenter, AXIS_X_ ) );
}
if( volume->Get( it.x, it.y-1, it.z ) != voxelMaterial )
{
activeEdges.Add( ComputeIntersection( _SDF, voxelCenter, AXIS_Y_ ) );
}
if( volume->Get( it.x, it.y, it.z-1 ) != voxelMaterial )
{
activeEdges.Add( ComputeIntersection( _SDF, voxelCenter, AXIS_Z_ ) );
}
Decompression is symmetric and can be performed simultaneously with countouring:
/// Experiments show that 3-byte normals don't have enough precision.
/// 8 bytes - GOOD
class GridEdge
{
UShort4 packed;
public:
GridEdge()
{
packed.v = 0; // initialize with invalid values
}
bool HasIntersection() const {
return !!packed.v;
}
void Encode( const V3F& normal, const float f01 )
{
mxASSERT(V3_IsNormalized(normal));
mxASSERT(f01 >= 0.0f && f01 <= 1.0f);
packed.x = TQuantize<16>::EncodeSNorm( normal.x );
packed.y = TQuantize<16>::EncodeSNorm( normal.y );
packed.z = TQuantize<16>::EncodeSNorm( normal.z );
packed.w = TQuantize<16>::EncodeUNorm( f01 );
}
const NormalDistance Decode() const
{
NormalDistance result;
result.N.x = TQuantize<16>::DecodeSNorm( packed.x );
result.N.y = TQuantize<16>::DecodeSNorm( packed.y );
result.N.z = TQuantize<16>::DecodeSNorm( packed.z );
result.d = TQuantize<16>::DecodeUNorm( packed.w );
//NOTE: normalization is necessary,
//otherwise the extracted surface will look noisy and perturbed (e.g. annoying notches)
float length;
result.N = V3_Normalized( result.N, length );
return result;
}
};
/// exact intersection points with normals (on zero-crossing edges only)
typedef DynamicArray< GridEdge > HermiteDataT;
/// DC cell grid is dual to the voxel grid (its corners are voxel centers)
struct DualCell {
EdgeIndex edges[NUM_AXES_]; //!< indices into Hermite data; -1 == no intersection
VertexIndex vertexId;
public:
DualCell() {
for( int i = 0; i < mxCOUNT_OF(edges); i++ ) {
edges[i] = NO_INTERSECTION;
}
vertexId = NO_VERTEX;
}
};
for( int iVoxelZ = StartZ; iVoxelZ <= MaxZ; iVoxelZ++ )
{
for( int iVoxelY = StartY; iVoxelY <= MaxY; iVoxelY++ )
{
for( int iVoxelX = StartX; iVoxelX <= MaxX; iVoxelX++ )
{
const Int3 it = { iVoxelX, iVoxelY, iVoxelZ };
const V3F voxelCenter = GetVoxelCenter( it.x, it.y, it.z );
const MaterialID voxelMaterial = volume.Get( it.x, it.y, it.z );
const AABB24 cellAabb = { voxelCenter - VOXEL_SIZE, voxelCenter };
// compute a vertex position for the current cell
// and emit quads for the previously visited cell
QEF_Solver::Input solverInput;
solverInput.bounds = cellAabb;
// Cube edge enumeration (edges are split into 3 groups by axes X,Y,Z, numbered using right-hand rule):
// NOTE: voxel centers are at cube corners:
// \WE'RE HERE!
// ______2_____\|/
// /| /| Corner vertices are numbered as follows:
// 5 |11 6 | 6___________7
// / | / |10 /| / Z
// |------3-----| | / | /| | Y
// | |_____1__|___| / | / | | /
// | / | / 4------------5 | |/
// 8| 4 9| 7 | 2________|__3 O-------X
// | / | / | / | /
// |/___________|/ | / | /
// 0 | / |/
// 0/___________1
//
DualCell & cell = dualCellGrid.Get( it.x, it.y, it.z );
// Edge 2:
if( volume.Get( it.x-1, it.y, it.z ) != voxelMaterial )
{
UnpackActiveEdge_And_AddIntersectionIfExists( voxelCenter, cell, AXIS_X_, hermiteData, solverInput );
}
// Edge 6:
if( volume.Get( it.x, it.y-1, it.z ) != voxelMaterial )
{
UnpackActiveEdge_And_AddIntersectionIfExists( voxelCenter, cell, AXIS_Y_, hermiteData, solverInput );
}
// Edge 10:
if( volume.Get( it.x, it.y, it.z-1 ) != voxelMaterial )
{
UnpackActiveEdge_And_AddIntersectionIfExists( voxelCenter, cell, AXIS_Z_, hermiteData, solverInput );
}
// Edges 1 and 7:
{
const DualCell& lowerCell = dualCellGrid.Get( it.x, it.y, it.z-1 );
const V3F lowerCellCenter = GetVoxelCenter( it.x, it.y, it.z-1 );
// 1
AddIntersectionIfExists( lowerCellCenter, lowerCell, AXIS_X_, hermiteData(), solverInput );
// 7
AddIntersectionIfExists( lowerCellCenter, lowerCell, AXIS_Y_, hermiteData(), solverInput );
}
// Edges 3 and 9:
{
const DualCell& frontCell = dualCellGrid.Get( it.x, it.y-1, it.z );
const V3F frontCellCenter = GetVoxelCenter( it.x, it.y-1, it.z );
// 3
AddIntersectionIfExists( frontCellCenter, frontCell, AXIS_X_, hermiteData(), solverInput );
// 9
AddIntersectionIfExists( frontCellCenter, frontCell, AXIS_Z_, hermiteData(), solverInput );
}
// Edges 5 and 11:
{
const DualCell& leftCell = dualCellGrid.Get( it.x-1, it.y, it.z );
const V3F leftCellCenter = GetVoxelCenter( it.x-1, it.y, it.z );
// 5
AddIntersectionIfExists( leftCellCenter, leftCell, AXIS_Y_, hermiteData(), solverInput );
// 11
AddIntersectionIfExists( leftCellCenter, leftCell, AXIS_Z_, hermiteData(), solverInput );
}
// Edge 8:
{
const DualCell& neighbor = dualCellGrid.Get( it.x-1, it.y-1, it.z );
const V3F neighborCenter = GetVoxelCenter( it.x-1, it.y-1, it.z );
// 8
AddIntersectionIfExists( neighborCenter, neighbor, AXIS_Z_, hermiteData(), solverInput );
}
// Edge 4:
{
const DualCell& neighbor = dualCellGrid.Get( it.x-1, it.y, it.z-1 );
const V3F neighborCenter = GetVoxelCenter( it.x-1, it.y, it.z-1 );
// 4
AddIntersectionIfExists( neighborCenter, neighbor, AXIS_Y_, hermiteData(), solverInput );
}
// Edge 0:
{
const DualCell& neighbor = dualCellGrid.Get( it.x, it.y-1, it.z-1 );
const V3F neighborCenter = GetVoxelCenter( it.x, it.y-1, it.z-1 );
// 0
AddIntersectionIfExists( neighborCenter, neighbor, AXIS_X_, hermiteData(), solverInput );
}
if( solverInput.numPoints > 0 )
{
QEF_Solver::Output solverOutput;
solver.Solve( solverInput, solverOutput );
cell.vertexId = mesh.AddVertex( solverOutput.position );
}//if( solverInput.numPoints > 0 )
// Emit quads for each intersecting edge.
const MaterialID v000 = volume.Get( it.x-1, it.y-1, it.z-1 );
const MaterialID v100 = volume.Get( it.x, it.y-1, it.z-1 );
const MaterialID v010 = volume.Get( it.x-1, it.y, it.z-1 );
const MaterialID v001 = volume.Get( it.x-1, it.y-1, it.z );
const bool originInside = (v000 != EMPTY_SPACE); // if the voxel material is not AIR
if( v000 != v100 )
{
// Connect the 4 cells sharing the edge along the X axis:
//
// .-----. Z
// / 0 /| | Y
// /-----/0| | /
// / 1 /|/| |/
// .-----|1|3| O-------X*
// | 1 |/|/
// |-----|2/
// | 2 |/
// .-----.
//
const DualCell& cell1 = dualCellGrid.Get( it.x, it.y-1, it.z );
const DualCell& cell2 = dualCellGrid.Get( it.x, it.y-1, it.z-1 );
const DualCell& cell3 = dualCellGrid.Get( it.x, it.y, it.z-1 );
if( originInside && CCW_FrontFacing ) {
CreateQuad( mesh, cell.vertexId, cell1.vertexId, cell2.vertexId, cell3.vertexId );
} else {
CreateQuad( mesh, cell3.vertexId, cell2.vertexId, cell1.vertexId, cell.vertexId );
}
}
if( v000 != v010 )
{
// Connect the 4 cells sharing the edge along the Y axis:
//
// .-----.-----. Z
// / 3 / 0 /| | Y*
// .-----+-----|0| | /
// | 3 | 0 |/| |/
// |-----X-----|1/ O-------X
// | 2 | 1 |/
// .-----.-----.
//
const DualCell& cell1 = dualCellGrid.Get( it.x, it.y, it.z-1 );
const DualCell& cell2 = dualCellGrid.Get( it.x-1, it.y, it.z-1 );
const DualCell& cell3 = dualCellGrid.Get( it.x-1, it.y, it.z );
if( originInside && CCW_FrontFacing ) {
CreateQuad( mesh, cell.vertexId, cell1.vertexId, cell2.vertexId, cell3.vertexId );
} else {
CreateQuad( mesh, cell3.vertexId, cell2.vertexId, cell1.vertexId, cell.vertexId );
}
}
if( v000 != v001 )
{
// Connect the 4 cells sharing the edge along the Z axis:
//
// Z*
// .-----.-----. | Y
// / 1 / 0 /| | /
// .-----.-----.0| |/
// / 2 / 3 /|/ O-------X
// .-----+-----|3/
// | 2 | 3 |/
// .-----.-----.
const DualCell& cell1 = dualCellGrid.Get( it.x-1, it.y, it.z );
const DualCell& cell2 = dualCellGrid.Get( it.x-1, it.y-1, it.z );
const DualCell& cell3 = dualCellGrid.Get( it.x, it.y-1, it.z );
if( originInside && CCW_FrontFacing ) {
CreateQuad( mesh, cell.vertexId, cell1.vertexId, cell2.vertexId, cell3.vertexId );
} else {
CreateQuad( mesh, cell3.vertexId, cell2.vertexId, cell1.vertexId, cell.vertexId );
}
}
}//X
}//Y
}//Z
How do I get the right order of vertices to form triangles that OpenGL can draw correctly from that large collection of vertices generated while finding the surface?
You simply create a quad for each active edge, connecting the four cells which share this edge (read Surface Nets / Dual Contouring). The orientation of the quad is easily derived from the active edge.
Then you split each quad to generate triangles for rendering:
void AddQuad(
const VertexIndex vertex0,
const VertexIndex vertex1,
const VertexIndex vertex2,
const VertexIndex vertex3
)
{
mxASSERT(vertex0 != NO_VERTEX);
mxASSERT(vertex1 != NO_VERTEX);
mxASSERT(vertex2 != NO_VERTEX);
mxASSERT(vertex3 != NO_VERTEX);
#if 0
const MeshTriangle tri1 = { vertex0, vertex1, vertex2 };
this->triangles.Add( tri1 );
const MeshTriangle tri2 = { vertex0, vertex2, vertex3 };
this->triangles.Add( tri2 );
#else
// Better results are obtained if we triangulate each quad by splitting it along the shortest diagonal.
// https://en.wikipedia.org/wiki/Distance_geometry_problem#Cayley.E2.80.93Menger_determinants
const float diagonal02 = DistanceBetween( vertices[vertex0].xyz, vertices[vertex2].xyz );
const float diagonal13 = DistanceBetween( vertices[vertex1].xyz, vertices[vertex3].xyz );
if( diagonal02 < diagonal13 ) {
const MeshTriangle tri1 = { vertex0, vertex1, vertex2 };
this->triangles.Add( tri1 );
const MeshTriangle tri2 = { vertex0, vertex2, vertex3 };
this->triangles.Add( tri2 );
} else {
const MeshTriangle tri1 = { vertex1, vertex3, vertex0 };
this->triangles.Add( tri1 );
const MeshTriangle tri2 = { vertex1, vertex2, vertex3 };
this->triangles.Add( tri2 );
}
#endif
}
HTH.
And what the hell is "the average of integrity along any type change"?
