OpenCL Naive Surface nets

BasicDensityCalc.cl

// Fills a grid of density values, the grid is vertex based, so for a 64^3 cells the grid must be 65^3 to include the extra vertex for the last cell
// This code expects to be combined with other code for a "DensityFunc" before being compiled.

#define VERTEX_SIZE 65

int vertex_index(const int4 pos)
{ 
    int x = clamp(pos.x, 0, VERTEX_SIZE - 1);
    int y = clamp(pos.y, 0, VERTEX_SIZE - 1);
    int z = clamp(pos.z, 0, VERTEX_SIZE - 1);
    return (z * VERTEX_SIZE * VERTEX_SIZE + y * VERTEX_SIZE + x);
}

// Signature for density function
//      float DensityFunc(float3 pos, const global float* paramData, const global float* transformData)


// This kernel fills a grid of density values based on the density function
kernel void GenerateDefaultField(
    //const int4 offset,
    global float* densities,
    global const float* shapeData,
    global const float4* transforms)
{
    const int x = get_global_id(0);
    const int y = get_global_id(1);
    const int z = get_global_id(2);
    
    if (x > VERTEX_SIZE || y > VERTEX_SIZE || z > VERTEX_SIZE)
        return;
    
    float3 world_pos = {x, y, z };
    // density function comes from elsewhere
    const float density = DensityFunc(world_pos, shapeData, transforms);
    
    const int4 local_pos = { x, y, z, 0 };
    const int index = vertex_index(local_pos);
    
    densities[index] = density;
}

ExampleUsage.cpp

    // Builds the vertex and index buffers, this could be ported to CL relatively easily
    //      however that would amount to a bunch of kernels that do little (2 scans, and 2 collapses based on them)
    // Based on https://github.com/nsf/mc, refer to it for offset_3d and offset_3d_slab functions

    void SurfaceNets::GenerateMesh(const float* densities, const Vec4* points, const unsigned char* writeMasks, const int size, VertexBuffer& vertexBuffer, IndexBuffer& indexBuffer)
    {
        std::vector<int> inds(65 * 65 * 2); // TODO: what's the actually maximum possible number of indices vs the actual plausible number?

        const Vec3 densitySizeVec(size + 1);
        const Vec3 sizeVec(size);
        for (int z = 0; z < size; ++z)
            for (int y = 0; y < size; ++y)
                for (int x = 0; x < size; ++x)
                {
                    Vec3 position(x, y, z);
                    unsigned cellPos = offset_3d(position, sizeVec);
                    unsigned char layout = writeMasks[cellPos];
                    if (layout == 0 || layout == 255)
                        continue;

                    const float density[8] = {
                        densities[offset_3d(position, densitySizeVec)],
                        densities[offset_3d(position + Vec3(1, 0, 0), densitySizeVec)],
                        densities[offset_3d(position + Vec3(0, 1, 0), densitySizeVec)],
                        densities[offset_3d(position + Vec3(1, 1, 0), densitySizeVec)],
                        densities[offset_3d(position + Vec3(0, 0, 1), densitySizeVec)],
                        densities[offset_3d(position + Vec3(1, 0, 1), densitySizeVec)],
                        densities[offset_3d(position + Vec3(0, 1, 1), densitySizeVec)],
                        densities[offset_3d(position + Vec3(1, 1, 1), densitySizeVec)]
                    };

                    Vec3 v = points[cellPos].xyz() / points[cellPos].w;
                    inds[offset_3d_slab(position, Vec3(65))] = vertexBuffer.size();
                    vertexBuffer.push_back({ v, Vec3(0, 0, 0) });

                    auto quad = [&](bool flip, int ia, int ib, int ic, int id)
                    {
                        if (flip)
                            std::swap(ib, id);

                        if (ia >= vertexBuffer.size() || ib >= vertexBuffer.size() || ic >= vertexBuffer.size() || id >= vertexBuffer.size())
                            return;

                        MeshVertex& va = vertexBuffer[ia];
                        MeshVertex& vb = vertexBuffer[ib];
                        MeshVertex& vc = vertexBuffer[ic];
                        MeshVertex& vd = vertexBuffer[id];

                        const Vec3 ab = va.Position - vb.Position;
                        const Vec3 cb = vc.Position - vb.Position;
                        const Vec3 n1 = cb.Cross(ab);
                        va.Normal += n1;
                        vb.Normal += n1;
                        vc.Normal += n1;

                        const Vec3 ac = va.Position - vc.Position;
                        const Vec3 dc = vd.Position - vc.Position;
                        const Vec3 n2 = dc.Cross(ac);
                        va.Normal += n2;
                        vc.Normal += n2;
                        vd.Normal += n2;

                        indexBuffer.push_back(ia);
                        indexBuffer.push_back(ib);
                        indexBuffer.push_back(ic);

                        indexBuffer.push_back(ia);
                        indexBuffer.push_back(ic);
                        indexBuffer.push_back(id);
                    };

                    const bool flip = density[0] < 0.0f;
                    if (position.y > 0 && position.z > 0 && (density[0] < 0.0f) != (density[1] < 0.0f)) {
                        quad(flip,
                            inds[offset_3d_slab(Vec3(position.x, position.y, position.z), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x, position.y, position.z - 1), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x, position.y - 1, position.z - 1), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x, position.y - 1, position.z), Vec3(65))]
                            );
                    }
                    if (position.x > 0 && position.z > 0 && (density[0] < 0.0f) != (density[2] < 0.0f)) {
                        quad(flip,
                            inds[offset_3d_slab(Vec3(position.x, position.y, position.z), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x - 1, position.y, position.z), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x - 1, position.y, position.z - 1), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x, position.y, position.z - 1), Vec3(65))]
                            );
                    }
                    if (position.x > 0 && position.y > 0 && (density[0] < 0.0f) != (density[4] < 0.0f)) {
                        quad(flip,
                            inds[offset_3d_slab(Vec3(position.x, position.y, position.z), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x, position.y - 1, position.z), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x - 1, position.y - 1, position.z), Vec3(65))],
                            inds[offset_3d_slab(Vec3(position.x - 1, position.y, position.z), Vec3(65))]
                            );
                    }
                }
    }

NaiveSurfaceNets.cl

// Calculates the vertex position for naive surface nets

// Number of vertices in the grid
#define EDGE_SIZE 65
// Number of cuberille cells in the grid
#define CELL_SIZE 64

int vertex_index(const int4 pos)
{ 
    int x = clamp(pos.x, 0, EDGE_SIZE - 1);
    int y = clamp(pos.y, 0, EDGE_SIZE - 1);
    int z = clamp(pos.z, 0, EDGE_SIZE - 1);
    return (z * EDGE_SIZE * EDGE_SIZE + y * EDGE_SIZE + x);
}

int cell_index(const int4 pos)
{
    int x = clamp(pos.x, 0, CELL_SIZE - 1);
    int y = clamp(pos.y, 0, CELL_SIZE - 1);
    int z = clamp(pos.z, 0, CELL_SIZE - 1);
    return (z * CELL_SIZE * CELL_SIZE + y * CELL_SIZE + x);
}

// Process an edge
float4 do_edge(float va, float vb, int axis, float3 pt)
{
    if (va < 0.0f == vb < 0.0f)
        return (float4){0,0,0,0};
    
    float value = va / (va - vb);
    if (axis == 0)
        return (float4){pt.x + value, pt.y, pt.z, 1.0f };
    else if (axis == 1)
        return (float4){pt.x, pt.y + value, pt.z, 1.0f };
    else
        return (float4){pt.x, pt.y, pt.z + value, 1.0f };
}

/// Takes a (DIM+1)^3 density field and computes the vertex positions in the cells using naive surface nets
kernel void CalculateVertexPositions(
    global const float* densities,      /* Density data received from the "BasicDensityCalculation.cl" kernel */
    global float4* positions,           /* contains the calculated mass-point vertex position */
    global uchar* writtenMask           /* Contains the 8bit corners mask */
)
{
    const int x = get_global_id(0);
    const int y = get_global_id(1);
    const int z = get_global_id(2);
    
    // Never trust the dispatcher, though this kernel shouldn't be able to crash and burn
    if (x > CELL_SIZE || y > CELL_SIZE || z > CELL_SIZE)
        return;
        
    int4 position = {x,y,z,0};
    
    // Grab all density at once (coherency)
    const float density[8] = {
                        densities[vertex_index(position)],
                        densities[vertex_index(position + (int4){1, 0, 0, 0})],
                        densities[vertex_index(position + (int4){0, 1, 0, 0})],
                        densities[vertex_index(position + (int4){1, 1, 0, 0})],
                        densities[vertex_index(position + (int4){0, 0, 1, 0})],
                        densities[vertex_index(position + (int4){1, 0, 1, 0})],
                        densities[vertex_index(position + (int4){0, 1, 1, 0})],
                        densities[vertex_index(position + (int4){1, 1, 1, 0})]
                    };
              
    // identify our local configuration, on the CPU we would break/return if the value was 0 or 255
    const int layout = ((density[0] < 0.0f) << 0) |
                        ((density[1] < 0.0f) << 1) |
                        ((density[2] < 0.0f) << 2) |
                        ((density[3] < 0.0f) << 3) |
                        ((density[4] < 0.0f) << 4) |
                        ((density[5] < 0.0f) << 5) |
                        ((density[6] < 0.0f) << 6) |
                        ((density[7] < 0.0f) << 7);
                    
    // Find the vertex position
    float4 vertPos = { 0.0f, 0.0f, 0.0f, 0.0f };
    vertPos += do_edge(density[0], density[1], 0, (float3){ x, y, z});
    vertPos += do_edge(density[2], density[3], 0, (float3){ x, y + 1, z});
    vertPos += do_edge(density[4], density[5], 0, (float3){ x, y, z + 1});
    vertPos += do_edge(density[6], density[7], 0, (float3){ x, y + 1, z + 1});
    vertPos += do_edge(density[0], density[2], 1, (float3){ x, y, z});
    vertPos += do_edge(density[1], density[3], 1, (float3){ x + 1, y, z});
    vertPos += do_edge(density[4], density[6], 1, (float3){ x, y, z + 1});
    vertPos += do_edge(density[5], density[7], 1, (float3){ x + 1, y, z + 1});
    vertPos += do_edge(density[0], density[4], 2, (float3){ x, y, z});
    vertPos += do_edge(density[1], density[5], 2, (float3){ x + 1, y, z});
    vertPos += do_edge(density[2], density[6], 2, (float3){ x, y + 1, z});
    vertPos += do_edge(density[3], density[7], 2, (float3){ x + 1, y + 1, z });
 
    // Write data, note that vertPos.w could very well be zero, but correct use of the "writtenMask" will never result in division by zero
    const int writeIndex = cell_index(position);
    positions[writeIndex] = vertPos;
    writtenMask[writeIndex] = layout;
    
    // Could perform 3 if tests for the need to output indices and output a uint of 1 to later run a scan on for indices, always outputs 2 triangles so there's really no need for any serious counting
    // Vertex count could be scanned off of the "writtenMask" any value that isn't 0 or 255 is a 1 in a prefix sum
}