nanovdb: add CUDA ConnectedComponents labeling for ValueOnIndex grids

Add nanovdb::tools::cuda::ConnectedComponents, which labels the active voxels
of a NanoVDB ValueOnIndex grid by 6-connectivity: two active voxels share a
label iff they are connected through a path of adjacent active voxels. The
public getVoxelLabelsAndCount() returns a per-active-voxel dense component id
in [0,N) plus the component count N.

The implementation is hierarchical:
  - per-leaf Shiloach-Vishkin union-find in shared memory (one block per leaf),
    yielding each leaf-local component's voxel Mask<3> and six face bitmasks;
  - cross-leaf edge detection by intersecting touching face masks over +X/+Y/+Z
    neighbours;
  - a global lock-free union-find over the resulting component graph.

Also add the ex_connected_components_cuda example (rasterize an .obj to a
narrow band via MeshToGrid, drop the sqrt(3)/2-voxel surface shell via
PruneGrid, then label; a CPU union-find oracle verifies the result) and a
ConnectedComponentsMultiSphere unit test.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Signed-off-by: Efty Sifakis <esifakis@nvidia.com>

Author: sifakis

nanovdb/nanovdb/CMakeLists.txt

@@ -237,6 +237,7 @@ set(NANOVDB_INCLUDE_TOOLS_CUDA_FILES
   tools/cuda/DilateGrid.cuh
   tools/cuda/CoarsenGrid.cuh
   tools/cuda/VoxelBlockManager.cuh
+  tools/cuda/ConnectedComponents.cuh
 )
 
 # NanoVDB util header files

nanovdb/nanovdb/examples/CMakeLists.txt

@@ -115,6 +115,7 @@ nanovdb_example(NAME "ex_merge_nanovdb_cuda" OPENVDB)
 nanovdb_example(NAME "ex_refine_nanovdb_cuda" OPENVDB)
 nanovdb_example(NAME "ex_coarsen_nanovdb_cuda" OPENVDB)
 nanovdb_example(NAME "ex_mesh_to_grid_cuda" OPENVDB)
+nanovdb_example(NAME "ex_connected_components_cuda")
 
 if(CUDAToolkit_FOUND)
   nanovdb_example(NAME "ex_make_mgpu_nanovdb") # requires cuRAND

nanovdb/nanovdb/examples/ex_connected_components_cuda/connected_components_cuda.cpp

@@ -0,0 +1,120 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: Apache-2.0
+
+/// @file  connected_components_cuda.cpp
+///
+/// @brief Host driver for the connected-components example (NanoVDB / CUDA only, no OpenVDB).
+///
+///        Reads a triangle mesh from a Wavefront .obj file, builds the index<->world transform,
+///        and hands the mesh to the CUDA side, which rasterizes it into a ValueOnIndex narrow-band
+///        grid, discards the surface/barrier shell, and runs connected-components labeling
+///        (nanovdb::tools::cuda::ConnectedComponents) on the result. The device side also runs a
+///        CPU union-find oracle that independently verifies the GPU labeling.
+///
+///        This example exercises ONLY the connected-components stage. For the full mesh->SDF
+///        pipeline that builds on top of it, see ex_mesh_to_sdf_cuda.
+
+#include <nanovdb/NanoVDB.h>          // host-usable: Vec3f, Vec3i, Vec3d, Map
+#include <nanovdb/GridHandle.h>       // GridHandle (header-only, host-usable)
+#include <nanovdb/cuda/DeviceBuffer.h>// nanovdb::cuda::DeviceBuffer
+
+#include <cstdint>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+// ---- Host/device seam (implemented in connected_components_cuda_kernels.cu) -----------------------
+//
+// Rasterize the mesh -> ValueOnIndex narrow band, discard the barrier shell (unsigned distance
+// within sqrt(3)/2 voxels of the surface), and run connected-components labeling. Prints topology
+// diagnostics, the component count, and a CPU-oracle PASS/FAIL, and returns the number of
+// connected components.
+uint64_t connectedComponentsFromMesh(const std::vector<nanovdb::Vec3f>& points,
+                                     const std::vector<nanovdb::Vec3i>& triangles,
+                                     const nanovdb::Map&                map,
+                                     float                              bandWidth);
+
+/// @brief Minimal Wavefront .obj reader (vertices + faces) using NanoVDB types.
+///
+///        Polygons with more than 3 vertices are fan-triangulated. Vertex references of the form
+///        `v`, `v/vt`, `v//vn`, `v/vt/vn` are accepted, as are negative (relative) indices. Lines
+///        that are not `v` or `f` are ignored.
+static void readOBJ(const std::string&            filename,
+                    std::vector<nanovdb::Vec3f>&  points,
+                    std::vector<nanovdb::Vec3i>&  triangles)
+{
+    std::ifstream file(filename);
+    if (!file.is_open())
+        throw std::runtime_error("Failed to open OBJ file: " + filename);
+
+    std::string line;
+    int lineNumber = 0;
+    while (std::getline(file, line)) {
+        ++lineNumber;
+        std::istringstream iss(line);
+        std::string type;
+        iss >> type;
+
+        if (type == "v") {
+            float x, y, z;
+            iss >> x >> y >> z;
+            points.emplace_back(x, y, z);
+        } else if (type == "f") {
+            std::vector<int> face;
+            std::string vert;
+            while (iss >> vert) {
+                const size_t slash = vert.find('/');
+                const std::string idxStr = vert.substr(0, slash);
+                if (idxStr.empty()) continue;
+                int raw = std::stoi(idxStr);
+                // OBJ indices are 1-based; negatives are relative to points read so far.
+                int idx = (raw < 0) ? int(points.size()) + raw : raw - 1;
+                if (idx < 0 || idx >= int(points.size()))
+                    throw std::runtime_error("OBJ parse error on line " +
+                                             std::to_string(lineNumber) +
+                                             ": face index out of bounds");
+                face.push_back(idx);
+            }
+            for (size_t i = 2; i < face.size(); ++i)
+                triangles.emplace_back(face[0], face[i - 1], face[i]);
+        }
+    }
+}
+
+int main(int argc, char* argv[])
+{
+    try {
+        if (argc < 2)
+            throw std::runtime_error("usage: " + std::string(argv[0]) +
+                                     " <input.obj> [voxelSize] [bandWidth]");
+
+        const std::string inputFile = argv[1];
+        const float voxelSize = (argc > 2) ? std::stof(argv[2]) : 0.01f;
+        const float bandWidth = (argc > 3) ? std::stof(argv[3]) : 3.0f;
+
+        std::vector<nanovdb::Vec3f> points;
+        std::vector<nanovdb::Vec3i> triangles;
+        std::cout << "Reading " << inputFile << "...\n";
+        readOBJ(inputFile, points, triangles);
+        std::cout << "Loaded " << points.size() << " vertices, "
+                  << triangles.size() << " triangles.\n";
+        if (points.empty() || triangles.empty())
+            throw std::runtime_error("mesh has no triangles");
+
+        // Index<->world transform: uniform voxel size, no translation.
+        nanovdb::Map map;
+        map.set(double(voxelSize), nanovdb::Vec3d(0.0), 1.0);
+
+        const uint64_t numComponents =
+            connectedComponentsFromMesh(points, triangles, map, bandWidth);
+        std::cout << "Connected components: " << numComponents << "\n";
+
+        return 0;
+    }
+    catch (const std::exception& e) {
+        std::cerr << "An exception occurred: \"" << e.what() << "\"\n";
+        return 1;
+    }
+}

nanovdb/nanovdb/examples/ex_connected_components_cuda/connected_components_cuda_kernels.cu

@@ -0,0 +1,228 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: Apache-2.0
+
+/// @file  connected_components_cuda_kernels.cu
+///
+/// @brief CUDA / NanoVDB side of the connected-components example. Rasterizes a triangle mesh into a
+///        ValueOnIndex narrow-band grid (nanovdb::tools::cuda::MeshToGrid), discards the
+///        surface/barrier shell (unsigned distance within sqrt(3)/2 voxels of the surface) with
+///        nanovdb::tools::cuda::PruneGrid, and runs connected-components labeling with
+///        nanovdb::tools::cuda::ConnectedComponents. A CPU union-find oracle independently verifies
+///        the GPU labeling (component count + per-voxel partition).
+
+#include <nanovdb/NanoVDB.h>
+#include <nanovdb/GridHandle.h>
+#include <nanovdb/cuda/DeviceBuffer.h>
+
+#include <nanovdb/tools/cuda/MeshToGrid.cuh>            // rasterize mesh -> ValueOnIndex + UDF
+#include <nanovdb/tools/cuda/PruneGrid.cuh>             // topological prune (barrier removal)
+#include <nanovdb/tools/cuda/ConnectedComponents.cuh>  // the connected-components labeling
+#include <nanovdb/util/cuda/Util.h>                     // operatorKernel, cudaCheck
+#include <nanovdb/util/cuda/DeviceGridTraits.cuh>       // DeviceGridTraits
+
+#include <thrust/universal_vector.h>
+
+#include <cstdint>
+#include <functional>
+#include <iostream>
+#include <numeric>
+#include <unordered_map>
+#include <vector>
+
+namespace {
+
+using BuildT      = nanovdb::ValueOnIndex;
+using GridHandleT = nanovdb::GridHandle<nanovdb::cuda::DeviceBuffer>;
+using Traits      = nanovdb::util::cuda::DeviceGridTraits<BuildT>;
+
+constexpr int LEAF_SIZE = 512;  // 8^3
+
+// Per-leaf retain-mask functor: a voxel is kept iff its unsigned distance to the surface exceeds the
+// barrier threshold sqrt(3)/2 voxels (i.e. UDF^2 >= 0.75 * voxelSize^2 in world units). Removing the
+// barrier shell splits each closed surface's narrow band into disjoint inner/outer shells, which is
+// what connected components then labels. One CUDA block per leaf, one thread per voxel offset.
+struct UDFBarrierPruneMaskFunctor
+{
+    static constexpr int MaxThreadsPerBlock         = LEAF_SIZE;
+    static constexpr int MinBlocksPerMultiprocessor = 1;
+
+    __device__ void operator()(const nanovdb::NanoGrid<BuildT>* d_grid,
+                               const float*                     d_udf,          // UDF sidecar, WORLD units
+                               float                            barrierSqWorld, // (sqrt(3)/2 * voxelSize)^2
+                               nanovdb::Mask<3>*                d_dstLeafMasks)
+    {
+        const int leafID   = blockIdx.x;
+        const int threadID = threadIdx.x;
+
+        const auto& leaf       = d_grid->tree().getFirstNode<0>()[leafID];
+        auto&       resultMask = d_dstLeafMasks[leafID];
+
+        // Clear the leaf's mask words in parallel, then set the retain bits.
+        if (threadID < nanovdb::Mask<3>::WORD_COUNT)
+            resultMask.words()[threadID] = 0UL;
+        __syncthreads();
+
+        if (auto n = leaf.data()->getValue(threadID)) {  // n != 0 => active voxel
+            const float udf = d_udf[n];
+            if (udf * udf >= barrierSqWorld)              // retain non-barrier voxels
+                resultMask.setOnAtomic(threadID);
+        }
+    }
+};
+
+// Pack a voxel coordinate into a sortable/ hashable int64 key (offset so negatives stay positive).
+inline int64_t encodeCoord(const nanovdb::Coord& c)
+{
+    return  (int64_t(c[0]) + (1 << 20))
+         | ((int64_t(c[1]) + (1 << 20)) << 21)
+         | ((int64_t(c[2]) + (1 << 20)) << 42);
+}
+
+// CPU union-find oracle: independently label the derived grid's active voxels by 6-connectivity and
+// verify the GPU result (a) has the same component count and (b) induces the same partition (two
+// voxels share a GPU label iff the oracle puts them in the same component). Returns true on PASS.
+bool validateAgainstOracle(const GridHandleT& derivedHandle, uint32_t leafCount, uint64_t active,
+                           const uint32_t* d_labels, uint64_t gpuCount)
+{
+    // Download the derived grid blob + the per-voxel labels to the host.
+    std::vector<char> blob(derivedHandle.bufferSize());
+    cudaCheck(cudaMemcpy(blob.data(), derivedHandle.deviceData(), blob.size(), cudaMemcpyDeviceToHost));
+    const auto* h_grid = reinterpret_cast<const nanovdb::NanoGrid<BuildT>*>(blob.data());
+
+    if (active == 0) {
+        std::cout << "CPU-oracle self-check: PASS (empty grid, 0 components)\n";
+        return gpuCount == 0;
+    }
+
+    std::vector<uint32_t> labels(active + 1);
+    cudaCheck(cudaMemcpy(labels.data(), d_labels, (active + 1) * sizeof(uint32_t), cudaMemcpyDeviceToHost));
+
+    // Gather active voxels: dense id -> {coord, gpu label}, plus a coord->id lookup for neighbours.
+    std::unordered_map<int64_t, uint32_t> coordToId;
+    std::vector<nanovdb::Coord>           idToCoord;
+    std::vector<uint32_t>                 idToGpuLabel;
+    coordToId.reserve(active * 2);
+    idToCoord.reserve(active);
+    idToGpuLabel.reserve(active);
+
+    const auto* leaves = h_grid->tree().getFirstLeaf();
+    for (uint32_t li = 0; li < leafCount; ++li) {
+        const auto& leaf = leaves[li];
+        for (uint32_t n = 0; n < uint32_t(LEAF_SIZE); ++n) {
+            if (!leaf.isActive(n)) continue;
+            const nanovdb::Coord c = leaf.origin() + nanovdb::NanoLeaf<BuildT>::OffsetToLocalCoord(n);
+            const uint64_t       slot = leaf.getValue(n);
+            const uint32_t       id   = uint32_t(idToCoord.size());
+            coordToId.emplace(encodeCoord(c), id);
+            idToCoord.push_back(c);
+            idToGpuLabel.push_back(labels[slot]);
+        }
+    }
+
+    // Union-find over 6-connectivity. Visiting only +X/+Y/+Z reaches every undirected edge once.
+    std::vector<uint32_t> parent(idToCoord.size());
+    std::iota(parent.begin(), parent.end(), 0u);
+    std::function<uint32_t(uint32_t)> find = [&](uint32_t x) {
+        while (parent[x] != x) { parent[x] = parent[parent[x]]; x = parent[x]; }
+        return x;
+    };
+    auto unite = [&](uint32_t a, uint32_t b) {
+        a = find(a); b = find(b);
+        if (a != b) parent[a > b ? a : b] = (a < b ? a : b);
+    };
+    const nanovdb::Coord dirs[3] = { {1,0,0}, {0,1,0}, {0,0,1} };
+    for (uint32_t id = 0; id < idToCoord.size(); ++id)
+        for (const auto& d : dirs) {
+            auto it = coordToId.find(encodeCoord(idToCoord[id] + d));
+            if (it != coordToId.end()) unite(id, it->second);
+        }
+
+    // Count oracle components and check the GPU labels induce the same partition.
+    std::unordered_map<uint32_t, uint32_t> rootToComp;   // oracle root -> dense component id
+    for (uint32_t id = 0; id < parent.size(); ++id) {
+        const uint32_t r = find(id);
+        if (!rootToComp.count(r)) rootToComp.emplace(r, uint32_t(rootToComp.size()));
+    }
+    const uint64_t oracleCount = rootToComp.size();
+
+    std::unordered_map<uint32_t, uint32_t> gpuToOracle;  // gpu label -> oracle component
+    uint64_t partitionViolations = 0;
+    for (uint32_t id = 0; id < parent.size(); ++id) {
+        const uint32_t gl = idToGpuLabel[id];
+        const uint32_t oc = rootToComp[find(id)];
+        auto it = gpuToOracle.find(gl);
+        if (it == gpuToOracle.end()) gpuToOracle.emplace(gl, oc);
+        else if (it->second != oc)   ++partitionViolations;
+    }
+    const uint64_t gpuDistinctLabels = gpuToOracle.size();
+
+    const bool pass = (gpuCount == oracleCount) &&
+                      (gpuDistinctLabels == oracleCount) &&
+                      (partitionViolations == 0);
+
+    std::cout << "CPU-oracle self-check: " << (pass ? "PASS" : "FAIL")
+              << "  (gpu=" << gpuCount << ", oracle=" << oracleCount
+              << ", distinct gpu labels=" << gpuDistinctLabels
+              << ", partition violations=" << partitionViolations << ")\n";
+    return pass;
+}
+
+} // anonymous namespace
+
+uint64_t connectedComponentsFromMesh(const std::vector<nanovdb::Vec3f>& points,
+                                     const std::vector<nanovdb::Vec3i>& triangles,
+                                     const nanovdb::Map&                map,
+                                     float                              bandWidth)
+{
+    const cudaStream_t stream = 0;
+
+    // ---- Step 1: rasterize the mesh -> ValueOnIndex narrow-band grid + UDF sidecar. ----
+    thrust::universal_vector<nanovdb::Vec3f> dPoints(points.begin(), points.end());
+    thrust::universal_vector<nanovdb::Vec3i> dTriangles(triangles.begin(), triangles.end());
+
+    nanovdb::tools::cuda::MeshToGrid<BuildT> converter(
+        dPoints.data().get(),    uint32_t(dPoints.size()),
+        dTriangles.data().get(), uint32_t(dTriangles.size()), map);
+    converter.setVerbose(1);
+    converter.setNarrowBandWidth(bandWidth);
+    auto [origHandle, udfSidecar] = converter.getHandleAndUDF();
+    const auto* d_orig = origHandle.template deviceGrid<BuildT>();
+
+    // World-space voxel size from the map (uniform scale here, but read it generically).
+    const nanovdb::Vec3d w0 = map.applyMap(nanovdb::Vec3d(0.0, 0.0, 0.0));
+    const nanovdb::Vec3d wx = map.applyMap(nanovdb::Vec3d(1.0, 0.0, 0.0));
+    const float voxelSize   = float(wx[0] - w0[0]);
+
+    // ---- Step 2: discard the surface/barrier shell -> derived topology. ----
+    const float    barrierSqWorld = 0.75f * voxelSize * voxelSize;  // (sqrt(3)/2 * voxelSize)^2
+    const uint32_t srcLeafCount   = Traits::getTreeData(d_orig).mNodeCount[0];
+
+    auto  retainMask   = nanovdb::cuda::DeviceBuffer::create(
+        std::size_t(srcLeafCount) * sizeof(nanovdb::Mask<3>), nullptr, false);
+    auto* d_retainMask = static_cast<nanovdb::Mask<3>*>(retainMask.deviceData());
+
+    nanovdb::util::cuda::operatorKernel<UDFBarrierPruneMaskFunctor>
+        <<<srcLeafCount, UDFBarrierPruneMaskFunctor::MaxThreadsPerBlock, 0, stream>>>(
+            d_orig, static_cast<const float*>(udfSidecar.deviceData()), barrierSqWorld, d_retainMask);
+    cudaCheckError();
+
+    nanovdb::tools::cuda::PruneGrid<BuildT> pruner(d_orig, d_retainMask, stream);
+    auto        derivedHandle = pruner.getHandle();
+    const auto* d_derived     = derivedHandle.template deviceGrid<BuildT>();
+
+    // ---- Step 3: connected-components labeling on the derived grid. ----
+    nanovdb::tools::cuda::ConnectedComponents<BuildT> cc(d_derived, stream);
+    auto [d_labels, numComponents] = cc.getVoxelLabelsAndCount();
+    cudaCheck(cudaStreamSynchronize(stream));
+
+    // Diagnostics.
+    const uint64_t derivedActive = Traits::getActiveVoxelCount(d_derived);
+    const uint32_t derivedLeaves = Traits::getTreeData(d_derived).mNodeCount[0];
+    std::cout << "Derived (barrier-removed) grid: " << derivedActive << " active voxels, "
+              << derivedLeaves << " leaves.\n";
+
+    // CPU-oracle self-check.
+    validateAgainstOracle(derivedHandle, derivedLeaves, derivedActive, d_labels, numComponents);
+
+    return numComponents;
+}