generate_sdf_mesh()

Kakshya::MeshData MayaFlux::Kinesis::generate_sdf_mesh	(	const Kinesis::SpatialField &	field,
		const glm::vec3 &	bounds_min,
		const glm::vec3 &	bounds_max,
		uint32_t	res_x,
		uint32_t	res_y,
		uint32_t	res_z,
		float	iso_level = 0.0F
	)

Extract an isosurface from a scalar field using marching cubes.

Evaluates field at every corner of a regular grid spanning [bounds_min, bounds_max] with res_x * res_y * res_z cells. Triangles are generated wherever the field crosses iso_level. Vertex normals are estimated via central finite differences on the field.

Any Kinesis::SpatialField is valid: distance(), combine(), chain(), or any user-supplied glm::vec3 -> float lambda wrapped in a SpatialField.

Parameters

field	SpatialField: glm::vec3 -> float.
bounds_min	World-space minimum corner of the evaluation volume.
bounds_max	World-space maximum corner of the evaluation volume.
res_x	Cell count along X. Clamped to minimum 1.
res_y	Cell count along Y. Clamped to minimum 1.
res_z	Cell count along Z. Clamped to minimum 1.
iso_level	Isosurface threshold value. Default 0.0.

Returns

MeshData ready for TRIANGLE_LIST draw, or empty MeshData if the field produces no crossings at iso_level.

Definition at line 68 of file MarchingCube.cpp.

{
    res_x = std::max(res_x, 1U);
    res_y = std::max(res_y, 1U);
    res_z = std::max(res_z, 1U);
 
    const glm::vec3 extent = bounds_max - bounds_min;
    const glm::vec3 step {
        extent.x / static_cast<float>(res_x),
        extent.y / static_cast<float>(res_y),
        extent.z / static_cast<float>(res_z),
    };
 
    const uint32_t nx = res_x + 1;
    const uint32_t ny = res_y + 1;
    const uint32_t nz = res_z + 1;
 
    const uint32_t px = nx + 2;
    const uint32_t py = ny + 2;
    const uint32_t pz = nz + 2;
    const size_t num_voxels = static_cast<size_t>(px) * py * pz;
 
    const size_t sx = 1;
    const size_t sy = px;
    const size_t sz = static_cast<size_t>(py) * px;
 
    std::vector<float> vals(num_voxels);
    for (uint32_t iz = 0; iz < pz; ++iz) {
        float z = bounds_min.z + static_cast<float>(static_cast<int32_t>(iz) - 1) * step.z;
        for (uint32_t iy = 0; iy < py; ++iy) {
            float y = bounds_min.y + static_cast<float>(static_cast<int32_t>(iy) - 1) * step.y;
            size_t base_idx = (static_cast<size_t>(iz) * py + iy) * px;
            for (uint32_t ix = 0; ix < px; ++ix) {
                float x = bounds_min.x + static_cast<float>(static_cast<int32_t>(ix) - 1) * step.x;
                vals[base_idx + ix] = field(glm::vec3(x, y, z));
            }
        }
    }
 
    std::vector<Kakshya::MeshVertex> verts;
    std::vector<uint32_t> indices;
    verts.reserve(static_cast<size_t>(res_x * res_y * res_z) * 3);
    indices.reserve(static_cast<size_t>(res_x * res_y * res_z) * 3);
 
    auto grad_at = [&](size_t ptr) -> glm::vec3 {
        glm::vec3 n {
            (vals[ptr + sx] - vals[ptr - sx]) / step.x,
            (vals[ptr + sy] - vals[ptr - sy]) / step.y,
            (vals[ptr + sz] - vals[ptr - sz]) / step.z
        };
        float len = glm::length(n);
        return len > 1e-7F ? (n / len) : glm::vec3(0.0F, 1.0F, 0.0F);
    };
 
    for (uint32_t iz = 0; iz < res_z; ++iz) {
        float z0 = bounds_min.z + static_cast<float>(iz) * step.z;
        float z1 = z0 + step.z;
 
        for (uint32_t iy = 0; iy < res_y; ++iy) {
            float y0 = bounds_min.y + static_cast<float>(iy) * step.y;
            float y1 = y0 + step.y;
 
            size_t voxel_ptr = (static_cast<size_t>(iz + 1) * py + (iy + 1)) * px + 1;
 
            for (uint32_t ix = 0; ix < res_x; ++ix, ++voxel_ptr) {
                const float v[8] = {
                    vals[voxel_ptr], vals[voxel_ptr + sx], vals[voxel_ptr + sy + sx], vals[voxel_ptr + sy],
                    vals[voxel_ptr + sz], vals[voxel_ptr + sz + sx], vals[voxel_ptr + sz + sy + sx], vals[voxel_ptr + sz + sy]
                };
 
                const uint32_t cube_idx = compute_cube_index(v, iso_level);
                const uint16_t edges = k_edge_table[cube_idx];
                if (edges == 0)
                    continue;
 
                float x0 = bounds_min.x + static_cast<float>(ix) * step.x;
                float x1 = x0 + step.x;
                const glm::vec3 c[8] = {
                    { x0, y0, z0 }, { x1, y0, z0 }, { x1, y1, z0 }, { x0, y1, z0 },
                    { x0, y0, z1 }, { x1, y0, z1 }, { x1, y1, z1 }, { x0, y1, z1 }
                };
 
                const glm::vec3 nrm[8] = {
                    grad_at(voxel_ptr), grad_at(voxel_ptr + sx), grad_at(voxel_ptr + sy + sx), grad_at(voxel_ptr + sy),
                    grad_at(voxel_ptr + sz), grad_at(voxel_ptr + sz + sx), grad_at(voxel_ptr + sz + sy + sx), grad_at(voxel_ptr + sz + sy)
                };
 
                InterpolatedEdge ep[12];
 
                auto interp_edge = [&](int i1, int i2) -> InterpolatedEdge {
                    float va = v[i1], vb = v[i2];
                    float d = vb - va;
                    float t = (std::abs(d) < 1e-7F) ? 0.5F : ((iso_level - va) / d);
 
                    return {
                        .position = c[i1] + t * (c[i2] - c[i1]),
                        .normal = glm::normalize(nrm[i1] + t * (nrm[i2] - nrm[i1]))
                    };
                };
 
                if (edges & 0x001)
                    ep[0] = interp_edge(0, 1);
                if (edges & 0x002)
                    ep[1] = interp_edge(1, 2);
                if (edges & 0x004)
                    ep[2] = interp_edge(2, 3);
                if (edges & 0x008)
                    ep[3] = interp_edge(3, 0);
                if (edges & 0x010)
                    ep[4] = interp_edge(4, 5);
                if (edges & 0x020)
                    ep[5] = interp_edge(5, 6);
                if (edges & 0x040)
                    ep[6] = interp_edge(6, 7);
                if (edges & 0x080)
                    ep[7] = interp_edge(7, 4);
                if (edges & 0x100)
                    ep[8] = interp_edge(0, 4);
                if (edges & 0x200)
                    ep[9] = interp_edge(1, 5);
                if (edges & 0x400)
                    ep[10] = interp_edge(2, 6);
                if (edges & 0x800)
                    ep[11] = interp_edge(3, 7);
 
                const auto& tri = k_tri_table[cube_idx];
 
                for (int i = 0; tri[i] != -1; i += 3) {
                    const auto base = static_cast<uint32_t>(verts.size());
 
                    for (int j = 0; j < 3; ++j) {
                        const auto& edge_data = ep[tri[i + j]];
                        glm::vec2 uv((edge_data.position.x - bounds_min.x) / extent.x,
                            (edge_data.position.y - bounds_min.y) / extent.y);
 
                        verts.push_back({ .position = edge_data.position,
                            .uv = uv,
                            .normal = edge_data.normal });
                    }
                    indices.insert(indices.end(), { base, base + 1, base + 2 });
                }
            }
        }
    }
 
    if (verts.empty())
        return Kakshya::MeshData::empty();
 
    auto data = Kakshya::MeshData::empty();
    Kakshya::MeshInsertion ins(data.vertex_variant, data.index_variant);
    ins.insert_flat(
        std::span<const uint8_t>(reinterpret_cast<const uint8_t*>(verts.data()), verts.size() * sizeof(Kakshya::MeshVertex)),
        std::span<const uint32_t>(indices),
        Kakshya::VertexLayout::for_meshes(sizeof(Kakshya::MeshVertex)));
 
    data.layout = Kakshya::VertexLayout::for_meshes(sizeof(Kakshya::MeshVertex));
    data.layout.vertex_count = static_cast<uint32_t>(verts.size());
    return data;
}

References compute_cube_index(), MayaFlux::Kakshya::MeshData::empty(), MayaFlux::Kakshya::VertexLayout::for_meshes(), MayaFlux::Kakshya::MeshInsertion::insert_flat(), k_edge_table, k_tri_table, and ptr.

Referenced by MayaFlux::Nodes::GpuSync::SDFNode::rebuild().

Here is the call graph for this function:

Here is the caller graph for this function: