dual contouring concept and term translations

Question

I am working on a 3D sandbox mmo with destructible terrain, using Java and LWJGL. I looked at most options to achieve this and found that dual contouring would probably be the best option for me. I've read about every paper i could find on the internet about how to implement such a thing and, since i have a very basic experience in game/3d/complex programming, I've had quite a hard time. but I think I've found the solution i was looking for and would like to check with some more experienced programmer if my solution seems right and if I really successfully wrapped my mind around some of those pesky terms that were haunting my nights for the last few months;

Edit: I will try and add some more precise detail as to what i already do here in bold

so, first of all, I've created an object called voxel when I say object, i refer the the user side's way of handling the data received from the server. I already compress the data as to have it fit a byte or a short and use bitwise operations. Voxels are loaded once from the server then stored to memory unless the server returns a change to that voxel (not sending entire chunks for a small change), a voxel knows it's type and it's "integrity" (in this case, how much that stone voxel has been mined for example. I think it refers the the volume i often read about?). Those voxels "will" be created from a pre generated world map for a persistent world. Yes, the world will have overhangs and most of all, tunnels players can mine.

then i got a chunk, composed of at most 32x32x32 voxels. A chunk knows it's world position and has a 3D array of voxels passed from the server when first loaded to memory, then the client/server only updates the voxels that have changed since the last update. A chunk also knows it's neighbors (if such neighbors have been loaded into memory, or creates a default chunk of 32x32x32 air voxels)

To find the surface, inside a chunk, i check every voxel from the bottom,left,near (I think it is the right order to fit with the way OpenGL read it's coordinates?) to the top,right,far. each time, i take the current voxels and all it's 7 neighbors in a cube shape and check whether their is a type change from solid to air in any of them, indicating a surface running through that cube. If there's a type change, i compare every edges to see if there's a type change and if so, place the entersection of that edge using and average of both corners integrity. doing so for all 12 edges allows me to average all the planes to a single point inside that cube and create a vertex and store it in an array.

The point where I'm confused is this: what should i do next? I think I should check along one edge that displays a type change to find the 3 other cubes that has the same edge and find their vertices. Since i know what direction the ground is (comparing the type change along the edge axis) I can deduce the front and back of the quad. taking 3 of those 4 vertices should allow me to get a normal from their dot product and from there, go around that normal CCW and send those vertices to the GPU for the rendering phase. I guess the real question here is: 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?

In theory, in my mind, it makes sense! but I find it so much simpler than all those papers I've found, that i fear I've missed the boat completely! I get lost when it comes to Isosurface (I guess here, the formula they're talking about would be the average of integrity along any type change?) and hermite data (here, should i read voxel collection instead?) and most maths are slightly beyond my understanding yet (working on it tho and I progressed quite a bit) please, can you help me out here and guide me a bit?

Answer 1

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"?