Massive improvements to vdb tool

nanovdb/nanovdb/NanoVDB.h

@@ -4719,7 +4719,7 @@ using OnIndexGrid = Grid<OnIndexTree>;
 * @endcode
 */
 
-/// @brief Use this function, which depends a pointer to GridData, to call
+/// @brief Use this function, which depends on a pointer to GridData, to call
 ///        other functions that depend on a NanoGrid of a known ValueType.
 /// @details This function allows for generic programming by converting GridData
 ///          to a NanoGrid of the type encoded in GridData::mGridType.
@@ -5762,6 +5762,26 @@ class ChannelAccessor : public DefaultReadAccessor<IndexT>
 
 }; // ChannelAccessor
 
+/// @brief Generic Accessor type that maps to either a ReadAccessor or ChannelAccessor
+/// @tparam BuildT Build type, e.g. float or ValueOnIndex
+/// @tparam ValueT Value type, e.g. float or Vec3f
+template <typename BuildT, typename ValueT>
+using AccType = typename util::conditional<BuildTraits<BuildT>::is_index,
+                ChannelAccessor<ValueT, BuildT>, DefaultReadAccessor<BuildT>>::type;
+
+/// @brief Generic template functions that return an Accessor to either an index grid or a regular grid
+template <typename GridT, typename ValueT>
+inline __hostdev__ auto getAccessor(const GridT &grid, ValueT *sideCar = nullptr)
+{
+    using BuildT = typename GridT::BuildType;
+    if constexpr(BuildTraits<BuildT>::is_index) {
+        return sideCar ? ChannelAccessor<ValueT, BuildT>(grid, sideCar) : ChannelAccessor<ValueT, BuildT>(grid);
+    } else {
+        static_assert(util::is_same<ValueT, typename GridT::ValueType>::value, "wrong ValueT for regular GridT");
+        return DefaultReadAccessor<BuildT>(grid);
+    }
+}
+
 #if 0
 // This MiniGridHandle class is only included as a stand-alone example. Note that aligned_alloc is a C++17 feature!
 // Normally we recommend using GridHandle defined in util/GridHandle.h but this minimal implementation could be an

nanovdb/nanovdb/examples/CMakeLists.txt

@@ -109,6 +109,7 @@ nanovdb_example(NAME "ex_bump_pool_buffer")
 nanovdb_example(NAME "ex_collide_level_set")
 nanovdb_example(NAME "ex_raytrace_fog_volume")
 nanovdb_example(NAME "ex_raytrace_level_set")
+nanovdb_example(NAME "ex_raytrace_iso_surface")
 nanovdb_example(NAME "ex_dilate_nanovdb_cuda" OPENVDB)
 nanovdb_example(NAME "ex_merge_nanovdb_cuda" OPENVDB)
 nanovdb_example(NAME "ex_refine_nanovdb_cuda" OPENVDB)

nanovdb/nanovdb/examples/ex_raytrace_iso_surface/common.h

@@ -0,0 +1,187 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: Apache-2.0
+
+#pragma once
+
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <chrono>
+#include <fstream>
+#include <nanovdb/NanoVDB.h>
+#include "ComputePrimitives.h"
+
+struct RenderOp;
+template<typename GridT>
+__global__ void renderIsoSurfacePersistentKernel(RenderOp renderOp, float* image, const GridT* grid, int numPixels, int* nextPixel);
+
+struct RenderOp
+{
+    using Vec3T  = nanovdb::math::Vec3<float>;
+    using RayT   = nanovdb::math::Ray<float>;
+    int mWidth, mHeight;
+    float mDx, mIso, mWBBoxDimZ;
+    Vec3T mWBBoxCenter;
+
+    template<typename BufferT>
+    RenderOp(nanovdb::GridHandle<BufferT>& handle, int width, int height)
+    {
+        mWidth = width;
+        mHeight = height;
+        const auto *metaData = handle.gridMetaData();
+        mDx = float(metaData->voxelSize()[0]);
+        mIso = mDx;
+        mWBBoxDimZ = (float)metaData->worldBBox().dim()[2] * 2;
+        mWBBoxCenter = Vec3T(metaData->worldBBox().min() + metaData->worldBBox().dim() * 0.5f);
+    }
+
+    template<typename GridT>
+    inline float renderImage(bool useCuda, float* image, const GridT* grid)
+    {
+        using ClockT = std::chrono::high_resolution_clock;
+        auto t0 = ClockT::now();
+
+        computeForEach(
+            useCuda, mWidth * mHeight, 512, __FILE__, __LINE__, [this, image, grid] __hostdev__(int start, int end) {
+            (*this)(start, end, image, grid);
+        });
+        computeSync(useCuda, __FILE__, __LINE__);
+
+        auto t1 = ClockT::now();
+        auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0).count() / 1000.f;
+        return duration;
+    }
+
+    template <typename GridT>
+    inline __hostdev__ void operator()(int start, int end, float* image, const GridT* grid) const
+    {
+        static_assert(nanovdb::util::is_same<typename GridT::BuildType, float, nanovdb::ValueOnIndex, nanovdb::ValueIndex>::value, "only works for float and OnIndex grids");
+        auto acc = nanovdb::getAccessor<GridT, float>(*grid);
+        for (int i = start; i < end; ++i) {
+            this->renderPixel(i, image, grid, acc);
+        }
+    }
+
+    template <typename GridT, typename AccT>
+    inline __hostdev__ void renderPixel(int i, float* image, const GridT* grid, AccT& acc) const
+    {
+        float          t0, v;
+        nanovdb::Coord ijk;
+        RayT           iRay = this->getIndexRay(i, grid);
+        if (nanovdb::math::isoCrossing(iRay, acc, ijk, v, t0, mIso)) {// intersect...
+            this->composite(image, i, (t0 * mDx) / (mWBBoxDimZ * 2), 1.0f);
+        } else {
+            this->composite(image, i, 0.0f, 0.0f);// write background value.
+        }
+    }
+
+    template<typename GridT>
+    inline float renderImagePersistent(float* image, const GridT* grid, int* nextPixel) const
+    {
+        int device = 0;
+        NANOVDB_CUDA_CHECK_ERROR(cudaGetDevice(&device), __FILE__, __LINE__);
+
+        cudaDeviceProp properties;
+        NANOVDB_CUDA_CHECK_ERROR(cudaGetDeviceProperties(&properties, device), __FILE__, __LINE__);
+
+        constexpr int blockSize = 256;
+        // Launch a small, fixed pool of blocks that persists on the GPU and
+        // pulls pixel work from a global counter instead of launching one
+        // logical thread per pixel up front.
+        int           blockCount = properties.multiProcessorCount * 4;
+        if (blockCount < 1) blockCount = 1;
+
+        // Reset the work queue before each timed render. The kernel advances
+        // this counter by one warp of pixels at a time.
+        NANOVDB_CUDA_CHECK_ERROR(cudaMemset(nextPixel, 0, sizeof(int)), __FILE__, __LINE__);
+
+        using ClockT = std::chrono::high_resolution_clock;
+        auto t0 = ClockT::now();
+
+        const int numPixels = mWidth * mHeight;
+        renderIsoSurfacePersistentKernel<GridT><<<blockCount, blockSize>>>(*this, image, grid, numPixels, nextPixel);
+        NANOVDB_CUDA_CHECK_ERROR(cudaGetLastError(), __FILE__, __LINE__);
+        NANOVDB_CUDA_CHECK_ERROR(cudaDeviceSynchronize(), __FILE__, __LINE__);
+
+        auto t1 = ClockT::now();
+        auto duration = std::chrono::duration_cast<std::chrono::microseconds>(t1 - t0).count() / 1000.f;
+        return duration;
+    }
+
+    template <typename GridT>
+    inline __hostdev__ RayT getIndexRay(int i, const GridT *grid) const
+    {
+        // perspective camera along Z-axis...
+        const uint32_t x = i % mWidth, y = i / mWidth;
+        const float fov = 45.f;
+        const float u = (float(x) + 0.5f) / mWidth;
+        const float v = (float(y) + 0.5f) / mHeight;
+        const float aspect = mWidth / float(mHeight);
+        const float Px = (2.f * u - 1.f) * tanf(fov / 2 * 3.14159265358979323846f / 180.f) * aspect;
+        const float Py = (2.f * v - 1.f) * tanf(fov / 2 * 3.14159265358979323846f / 180.f);
+        const Vec3T origin = mWBBoxCenter + Vec3T(0, 0, mWBBoxDimZ);
+        Vec3T       dir(Px, Py, -1.f);
+        dir.normalize();
+        RayT wRay(origin, dir);
+        return wRay.worldToIndexF(*grid);// transform the ray to the grid's index-space.
+    }
+
+    inline __hostdev__ void composite(float* outImage, int offset, float value, float alpha) const
+    {
+        const uint32_t x = offset % mWidth, y = offset / mWidth;
+
+        // checkerboard background...
+        const int   mask = 1 << 7;
+        const float bg = ((x & mask) ^ (y & mask)) ? 1.0f : 0.5f;
+        outImage[offset] = alpha * value + (1.0f - alpha) * bg;
+    }
+
+    inline void saveImage(const std::string& filename, const float* image) const
+    {
+        const auto isLittleEndian = []() -> bool {
+            static int  x = 1;
+            static bool result = reinterpret_cast<uint8_t*>(&x)[0] == 1;
+            return result;
+        };
+
+        float scale = 1.0f;
+        if (isLittleEndian()) scale = -scale;
+
+        std::fstream fs(filename, std::ios::out | std::ios::binary);
+        if (!fs.is_open()) throw std::runtime_error("Unable to open file: " + filename);
+
+        fs << "Pf\n"
+           << mWidth << "\n"
+           << mHeight << "\n"
+           << scale << "\n";
+
+        for (int i = 0; i < mWidth * mHeight; ++i) {
+            float r = image[i];
+            fs.write((char*)&r, sizeof(float));
+        }
+    }
+};
+
+template<typename GridT>
+__global__ void renderIsoSurfacePersistentKernel(RenderOp renderOp, float* image, const GridT* grid, int numPixels, int* nextPixel)
+{
+    static_assert(nanovdb::util::is_same<typename GridT::BuildType, float, nanovdb::ValueOnIndex, nanovdb::ValueIndex>::value, "only works for float and OnIndex grids");
+    auto acc = nanovdb::getAccessor<GridT, float>(*grid);
+    const unsigned int lane = threadIdx.x & 31u;
+
+    // Keep the fixed set of launched threads busy until all pixels have been assigned.
+    while (true) {
+        int base = 0;
+        // Each warp asks the shared counter for the next batch of 32 pixels.
+        // Only lane 0 updates the counter; __shfl_sync copies lane 0's result
+        // to the other lanes in the warp.
+        if (lane == 0) base = atomicAdd(nextPixel, 32);
+        base = __shfl_sync(0xFFFFFFFFu, base, 0);
+
+        // Each lane renders one pixel from the batch: lane 0 renders base,
+        // lane 1 renders base + 1, and so on.
+        const int i = base + int(lane);
+        if (i >= numPixels) break;
+
+        renderOp.renderPixel(i, image, grid, acc);
+    }
+}

nanovdb/nanovdb/examples/ex_raytrace_iso_surface/main.cc

@@ -0,0 +1,51 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: Apache-2.0
+
+#include <algorithm>
+#include <iostream>
+#include <nanovdb/io/IO.h>
+#include <nanovdb/tools/CreatePrimitives.h>
+#include <nanovdb/cuda/DeviceBuffer.h>
+
+#if defined(NANOVDB_USE_CUDA)
+using BufferT = nanovdb::cuda::DeviceBuffer;
+#else
+using BufferT = nanovdb::HostBuffer;
+#endif
+
+extern void runNanoVDB(nanovdb::GridHandle<BufferT>& handle, int numIterations, int width, int height, BufferT& imageBuffer, bool usePersistentThreads);
+
+int main(int ac, char** av)
+{
+    try {
+        bool        usePersistentThreads = false;
+        const char* gridName = nullptr;
+        for (int i = 1; i < ac; ++i) {
+            if (std::strcmp(av[i], "--persistent") == 0) {
+                usePersistentThreads = true;
+            } else if (!gridName) {
+                gridName = av[i];
+            } else {
+                throw std::runtime_error("Usage: ex_raytrace_iso_surface [--persistent] [grid.nvdb]");
+            }
+        }
+        nanovdb::GridHandle<BufferT> handle;
+        if (gridName) {
+            handle = nanovdb::io::readGrid<BufferT>(gridName);
+            std::cout << "Loaded NanoVDB grid[" << handle.gridMetaData()->shortGridName() << "]...\n";
+        } else {
+            handle = nanovdb::tools::createLevelSetSphere<float, BufferT>(100.0f, nanovdb::Vec3d(-20, 0, 0), 1.0, 3.0, nanovdb::Vec3d(0), "sphere");
+        }
+
+        const int numIterations = 50;
+        const int width  = 4096;
+        const int height = 4096;
+        BufferT imageBuffer(width * height * sizeof(float));
+
+        runNanoVDB(handle, numIterations, width, height, imageBuffer, usePersistentThreads);
+    }
+    catch (const std::exception& e) {
+        std::cerr << "An exception occurred: \"" << e.what() << "\"" << std::endl;
+    }
+    return 0;
+}

nanovdb/nanovdb/examples/ex_raytrace_iso_surface/nanovdb.cu

@@ -0,0 +1,66 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: Apache-2.0
+
+#define _USE_MATH_DEFINES
+#include <cmath>
+#include <chrono>
+
+#if defined(NANOVDB_USE_CUDA)
+#include <nanovdb/cuda/DeviceBuffer.h>
+using BufferT = nanovdb::cuda::DeviceBuffer;
+#else
+using BufferT = nanovdb::HostBuffer;
+#endif
+#include <nanovdb/GridHandle.h>
+#include <nanovdb/io/IO.h>
+#include <nanovdb/math/Ray.h>
+#include <nanovdb/math/HDDA.h>
+
+#include "common.h"
+
+void runNanoVDB(nanovdb::GridHandle<BufferT>& handle, int numIterations, int width, int height, BufferT& imageBuffer, bool usePersistentThreads)
+{
+    float *h_outImage = reinterpret_cast<float*>(imageBuffer.data());
+    RenderOp renderOp(handle, width, height);
+
+    auto kernel = [&](auto *h_grid){
+        float sum = 0;
+        for (int i = 0; i < numIterations; ++i, sum += renderOp.renderImage(false/*useCuda*/, h_outImage, h_grid));
+        std::cout << "Average of " << numIterations << " renderings (NanoVDB-Host) = " << (sum/numIterations) << " ms" << std::endl;
+        renderOp.saveImage("raytrace_iso_surface-nanovdb-host.pfm", (float*)imageBuffer.data());
+
+#if defined(NANOVDB_USE_CUDA)
+        handle.deviceUpload();
+        using BuildT = typename nanovdb::util::remove_pointer_t<decltype(h_grid)>::BuildType;
+        auto* d_grid = handle.deviceGrid<BuildT>();
+        if (!d_grid) throw std::runtime_error("GridHandle does not contain a valid device grid");
+        imageBuffer.deviceUpload();
+        float* d_outImage = reinterpret_cast<float*>(imageBuffer.deviceData());
+        sum = 0;
+        if (usePersistentThreads) {
+            int* d_nextPixel = nullptr;
+            NANOVDB_CUDA_CHECK_ERROR(cudaMalloc(&d_nextPixel, sizeof(int)), __FILE__, __LINE__);
+            for (int i = 0; i < numIterations; ++i, sum += renderOp.renderImagePersistent(d_outImage, d_grid, d_nextPixel));
+            NANOVDB_CUDA_CHECK_ERROR(cudaFree(d_nextPixel), __FILE__, __LINE__);
+            std::cout << "Average of " << numIterations << " renderings (NanoVDB-Cuda-Persistent) = " << (sum/numIterations) << " ms " << std::endl;
+            imageBuffer.deviceDownload();
+            renderOp.saveImage("raytrace_iso_surface-nanovdb-cuda-persistent.pfm", (float*)imageBuffer.data());
+        } else {
+            for (int i = 0; i < numIterations; ++i, sum += renderOp.renderImage(true/*useCuda*/, d_outImage, d_grid));
+            std::cout << "Average of " << numIterations << " renderings (NanoVDB-Cuda) = " << (sum/numIterations) << " ms " << std::endl;
+            imageBuffer.deviceDownload();
+            renderOp.saveImage("raytrace_iso_surface-nanovdb-cuda.pfm", (float*)imageBuffer.data());
+        }
+#endif
+    };// kernel
+
+    if (auto *h_grid = handle.grid<float>()) {
+        kernel(h_grid);
+    } else if (auto *h_grid = handle.grid<nanovdb::ValueIndex>()) {
+        kernel(h_grid);
+    } else if (auto *h_grid = handle.grid<nanovdb::ValueOnIndex>()) {
+        kernel(h_grid);
+    } else {
+        throw std::runtime_error("GridHandle does not contain a valid device grid");
+    }
+}// runNanoVDB

nanovdb/nanovdb/examples/ex_raytrace_level_set/nanovdb.cu

@@ -20,15 +20,14 @@ using BufferT = nanovdb::HostBuffer;
 
 void runNanoVDB(nanovdb::GridHandle<BufferT>& handle, int numIterations, int width, int height, BufferT& imageBuffer)
 {
-    using GridT = nanovdb::FloatGrid;
+    using GridT  = nanovdb::FloatGrid;
     using CoordT = nanovdb::Coord;
-    using RealT = float;
-    using Vec3T = nanovdb::math::Vec3<RealT>;
-    using RayT = nanovdb::math::Ray<RealT>;
+    using RealT  = float;
+    using Vec3T  = nanovdb::math::Vec3<RealT>;
+    using RayT   = nanovdb::math::Ray<RealT>;
 
     auto *h_grid = handle.grid<float>();
-    if (!h_grid)
-        throw std::runtime_error("GridHandle does not contain a valid host grid");
+    if (!h_grid) throw std::runtime_error("GridHandle does not contain a valid host grid");
 
     float* h_outImage = reinterpret_cast<float*>(imageBuffer.data());
 
@@ -58,7 +57,7 @@ void runNanoVDB(nanovdb::GridHandle<BufferT>& handle, int numIterations, int wid
             float  t0;
             CoordT ijk;
             float  v;
-            if (nanovdb::math::ZeroCrossing(iRay, acc, ijk, v, t0)) {
+            if (nanovdb::math::zeroCrossing(iRay, acc, ijk, v, t0)) {
                 // write distance to surface. (we assume it is a uniform voxel)
                 float wT0 = t0 * float(grid->voxelSize()[0]);
                 compositeOp(image, i, width, height, wT0 / (wBBoxDimZ * 2), 1.0f);