Draft PR: Half Grid Support for OpenVDB

CMakeLists.txt

@@ -68,7 +68,7 @@ include(GNUInstallDirs)
 option(OPENVDB_BUILD_CORE "Enable the core OpenVDB library. Both static and shared versions are enabled by default" ON)
 option(OPENVDB_BUILD_BINARIES "Enable the vdb binaries. Only vdb_print is enabled by default" ON)
 option(OPENVDB_BUILD_PYTHON_MODULE "Build the pyopenvdb Python module" OFF)
-option(OPENVDB_BUILD_UNITTESTS "Build the OpenVDB unit tests" OFF)
+option(OPENVDB_BUILD_UNITTESTS "Build the OpenVDB unit tests" ON)
 option(OPENVDB_BUILD_DOCS "Build the OpenVDB documentation" OFF)
 option(OPENVDB_BUILD_HOUDINI_PLUGIN "Build the Houdini plugin" OFF)
 option(OPENVDB_BUILD_HOUDINI_ABITESTS "Build the Houdini ABI tests" OFF)

openvdb/openvdb/Metadata.h

@@ -358,6 +358,7 @@ operator<<(std::ostream& ostr, const Metadata& metadata)
 using BoolMetadata   = TypedMetadata<bool>;
 using DoubleMetadata = TypedMetadata<double>;
 using FloatMetadata  = TypedMetadata<float>;
+using HalfMetadata   = TypedMetadata<math::half>;
 using Int32Metadata  = TypedMetadata<int32_t>;
 using Int64Metadata  = TypedMetadata<int64_t>;
 using StringMetadata = TypedMetadata<std::string>;

openvdb/openvdb/Types.h

@@ -58,6 +58,7 @@ using Int64   = int64_t;
 using Int     = Int32;
 using Byte    = unsigned char;
 using Real    = double;
+using Half    = math::half;
 
 // Two-dimensional vector types
 using Vec2R = math::Vec2<Real>;
@@ -446,6 +447,29 @@ template<typename FromType, typename ToType> struct CopyConstness<const FromType
 /// @endcond
 
 
+////////////////////////////////////////
+template<class T>
+struct is_floating_point : std::is_floating_point<T> { };
+
+template<>
+struct is_floating_point<math::half> : std::is_floating_point<float> { };
+
+/// @brief Maps one type (e.g. the value types) to other types
+template<typename T>
+struct ValueToComputeMap
+{
+    using Type = T;
+    using type = T;
+};
+
+template<>
+struct ValueToComputeMap<math::half>
+{
+    using Type = float;
+    using type = float;
+};
+
+
 ////////////////////////////////////////
 
 
@@ -688,6 +712,23 @@ class Steal {};
 /// @brief Tag dispatch class that distinguishes constructors during file input
 class PartialCreate {};
 
+// For half compilation
+namespace math {
+template<>
+inline auto cwiseAdd(const math::Vec3<math::half>& v, const float s)
+{
+    math::Vec3<math::half> out;
+    const math::half* ip = v.asPointer();
+    math::half* op = out.asPointer();
+    for (unsigned i = 0; i < 3; ++i, ++op, ++ip) {
+        OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
+        *op = *ip + s;
+        OPENVDB_NO_TYPE_CONVERSION_WARNING_END
+    }
+    return out;
+}
+} // namespace math
+
 } // namespace OPENVDB_VERSION_NAME
 } // namespace openvdb
 

openvdb/openvdb/math/Vec3.h

@@ -11,7 +11,6 @@
 #include <cmath>
 #include <type_traits>
 
-
 namespace openvdb {
 OPENVDB_USE_VERSION_NAMESPACE
 namespace OPENVDB_VERSION_NAME {
@@ -668,6 +667,7 @@ OPENVDB_IS_POD(Vec3ui)
 OPENVDB_IS_POD(Vec3s)
 OPENVDB_IS_POD(Vec3d)
 
+
 } // namespace math
 } // namespace OPENVDB_VERSION_NAME
 } // namespace openvdb

openvdb/openvdb/openvdb.h

@@ -53,6 +53,7 @@ namespace io { class DelayedLoadMetadata; }
 using BoolTree     = tree::Tree4<bool,        5, 4, 3>::Type;
 using DoubleTree   = tree::Tree4<double,      5, 4, 3>::Type;
 using FloatTree    = tree::Tree4<float,       5, 4, 3>::Type;
+using HalfTree     = tree::Tree4<Half,        5, 4, 3>::Type;
 using Int32Tree    = tree::Tree4<int32_t,     5, 4, 3>::Type;
 using Int64Tree    = tree::Tree4<int64_t,     5, 4, 3>::Type;
 using MaskTree     = tree::Tree4<ValueMask,   5, 4, 3>::Type;
@@ -73,6 +74,7 @@ using VectorTree   = Vec3fTree;
 using BoolGrid     = Grid<BoolTree>;
 using DoubleGrid   = Grid<DoubleTree>;
 using FloatGrid    = Grid<FloatTree>;
+using HalfGrid     = Grid<HalfTree>;
 using Int32Grid    = Grid<Int32Tree>;
 using Int64Grid    = Grid<Int64Tree>;
 using MaskGrid     = Grid<MaskTree>;
@@ -88,7 +90,7 @@ using VectorGrid   = Vec3fGrid;
 /// @name Lists of native Grid Types
 /// @{
 /// The floating point Grid types which OpenVDB will register by default.
-using RealGridTypes   = TypeList<FloatGrid, DoubleGrid>;
+using RealGridTypes   = TypeList<HalfGrid, FloatGrid, DoubleGrid>;
 /// The integer Grid types which OpenVDB will register by default.
 using IntegerGridTypes = TypeList<Int32Grid, Int64Grid>;
 /// The scalar Grid types which OpenVDB will register by default. This is a
@@ -204,6 +206,7 @@ using MetaTypes = TypeList<
     BoolMetadata,
     DoubleMetadata,
     FloatMetadata,
+    HalfMetadata,
     Int32Metadata,
     Int64Metadata,
     StringMetadata,

openvdb/openvdb/tools/Interpolation.h

@@ -40,14 +40,16 @@
 #ifndef OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
 #define OPENVDB_TOOLS_INTERPOLATION_HAS_BEEN_INCLUDED
 
-#include <openvdb/version.h> // for OPENVDB_VERSION_NAME
+#include <openvdb/Grid.h>
 #include <openvdb/Platform.h> // for round()
+#include <openvdb/Types.h> // for ValueToComputeMap
 #include <openvdb/math/Math.h>// for SmoothUnitStep
 #include <openvdb/math/Transform.h> // for Transform
-#include <openvdb/Grid.h>
 #include <openvdb/tree/ValueAccessor.h>
+#include <openvdb/version.h> // for OPENVDB_VERSION_NAME
 #include <cmath>
 #include <type_traits>
+#include <iostream>
 
 namespace openvdb {
 OPENVDB_USE_VERSION_NAMESPACE
@@ -313,7 +315,7 @@ class GridSampler
     /// @brief Sample value in integer index space
     /// @param i Integer x-coordinate in index space
     /// @param j Integer y-coordinate in index space
-    /// @param k Integer x-coordinate in index space
+    /// @param k Integer z-coordinate in index space
     ValueType sampleVoxel(typename Coord::ValueType i,
                           typename Coord::ValueType j,
                           typename Coord::ValueType k) const
@@ -392,7 +394,7 @@ class GridSampler<tree::ValueAccessor<TreeT>, SamplerType>
     /// @brief Sample value in integer index space
     /// @param i Integer x-coordinate in index space
     /// @param j Integer y-coordinate in index space
-    /// @param k Integer x-coordinate in index space
+    /// @param k Integer z-coordinate in index space
     ValueType sampleVoxel(typename Coord::ValueType i,
                           typename Coord::ValueType j,
                           typename Coord::ValueType k) const
@@ -744,17 +746,31 @@ BoxSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
                    typename TreeT::ValueType& result)
 {
     using ValueT = typename TreeT::ValueType;
+    using ComputeT = typename ValueToComputeMap<ValueT>::Type;
 
     const Vec3i inIdx = local_util::floorVec3(inCoord);
     const Vec3R uvw = inCoord - inIdx;
 
     // Retrieve the values of the eight voxels surrounding the
     // fractional source coordinates.
-    ValueT data[2][2][2];
-
-    const bool hasActiveValues = BoxSampler::probeValues(data, inTree, Coord(inIdx));
+    ValueT probeData[2][2][2];
+    ComputeT computeData[2][2][2];
+    bool hasActiveValues;
+
+    if constexpr(std::is_same<ValueT, ComputeT>::value) {
+        hasActiveValues = BoxSampler::probeValues(computeData, inTree, Coord(inIdx));
+    } else {
+        hasActiveValues = BoxSampler::probeValues(probeData, inTree, Coord(inIdx));
+        for (int dx = 0; dx < 2; ++dx) {
+            for (int dy = 0; dy < 2; ++dy) {
+                for (int dz = 0; dz < 2; ++dz) {
+                    computeData[dx][dy][dz] = probeData[dx][dy][dz];
+                }
+            }
+        }
+    }
 
-    result = BoxSampler::trilinearInterpolation(data, uvw);
+    result = BoxSampler::trilinearInterpolation(computeData, uvw);
 
     return hasActiveValues;
 }
@@ -831,6 +847,7 @@ QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
     typename TreeT::ValueType& result)
 {
     using ValueT = typename TreeT::ValueType;
+    using ComputeT = typename ValueToComputeMap<ValueT>::Type;
 
     const Vec3i inIdx = local_util::floorVec3(inCoord), inLoIdx = inIdx - Vec3i(1, 1, 1);
     const Vec3R uvw = inCoord - inIdx;
@@ -839,15 +856,17 @@ QuadraticSampler::sample(const TreeT& inTree, const Vec3R& inCoord,
     // fractional source coordinates.
     bool active = false;
     ValueT data[3][3][3];
+    ComputeT computeData[3][3][3];
     for (int dx = 0, ix = inLoIdx.x(); dx < 3; ++dx, ++ix) {
         for (int dy = 0, iy = inLoIdx.y(); dy < 3; ++dy, ++iy) {
             for (int dz = 0, iz = inLoIdx.z(); dz < 3; ++dz, ++iz) {
                 if (inTree.probeValue(Coord(ix, iy, iz), data[dx][dy][dz])) active = true;
+                computeData[dx][dy][dz] = data[dx][dy][dz];
             }
         }
     }
 
-    result = QuadraticSampler::triquadraticInterpolation(data, uvw);
+    result = QuadraticSampler::triquadraticInterpolation(computeData, uvw);
 
     return active;
 }

openvdb/openvdb/tools/RayIntersector.h

@@ -90,7 +90,7 @@ class LevelSetRayIntersector
     using TreeT = typename GridT::TreeType;
 
     static_assert(NodeLevel >= -1 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range");
-    static_assert(std::is_floating_point<ValueT>::value,
+    static_assert(openvdb::is_floating_point<ValueT>::value,
         "level set grids must have scalar, floating-point value types");
 
     /// @brief Constructor

openvdb/openvdb/tools/RayTracer.h

@@ -154,7 +154,7 @@ class VolumeRender
     using ValueType = typename GridType::ValueType;
     using AccessorType = typename GridType::ConstAccessor;
     using SamplerType = tools::GridSampler<AccessorType, SamplerT>;
-    static_assert(std::is_floating_point<ValueType>::value,
+    static_assert(openvdb::is_floating_point<ValueType>::value,
         "VolumeRender requires a floating-point-valued grid");
 
     /// @brief Constructor taking an intersector and a base camera.

openvdb/openvdb/unittest/TestGridIO.cc

@@ -25,6 +25,7 @@ class TestGridIO: public ::testing::Test
 
 protected:
     template<typename GridType> void readAllTest();
+    void readWriteHalfGridTest();
 };
 
 
@@ -190,7 +191,82 @@ TestGridIO::readAllTest()
     ::remove("something.vdb2");
 }
 
+
+////////////////////////////////////////
+
+
+void
+TestGridIO::readWriteHalfGridTest()
+{
+    using namespace openvdb;
+
+    float TEST_FLOAT = 65504.f + 100.f;
+    openvdb::math::half TEST_HALF = 65504.f + 100.f;
+    openvdb::Coord xyz(1000, -200000000, 30000000);
+
+    std::string fileName = "testHalfAndFloat.vdb";
+    openvdb::FloatGrid::Ptr fgfd = openvdb::FloatGrid::create(2.f); // float grid stored as float on disk
+    fgfd->setTransform(openvdb::math::Transform::createLinearTransform(/*voxel size=*/1.0));
+    fgfd->setName("float_grid");
+
+    openvdb::HalfGrid::Ptr hghd = openvdb::HalfGrid::create(2.f); // half grid stored as half on disk
+    hghd->setTransform(openvdb::math::Transform::createLinearTransform(/*voxel size=*/1.0));
+    hghd->setName("half_grid");
+
+    openvdb::FloatGrid::Ptr fghd = openvdb::FloatGrid::create(2.f); // float grid stored as half on disk
+    fghd->setTransform(openvdb::math::Transform::createLinearTransform(/*voxel size=*/1.0));
+    fghd->setName("float_grid_as_half");
+    fghd->setSaveFloatAsHalf(true);
+
+    {
+        openvdb::FloatGrid::Accessor fgfdAcc = fgfd->getAccessor();
+        openvdb::HalfGrid::Accessor hghdAcc = hghd->getAccessor();
+        openvdb::FloatGrid::Accessor fghdAcc = fghd->getAccessor();
+        fgfdAcc.setValue(xyz, TEST_FLOAT);
+        hghdAcc.setValue(xyz, TEST_HALF);
+        fghdAcc.setValue(xyz, TEST_FLOAT);
+    }
+
+    openvdb::io::File writeFile(fileName);
+    openvdb::GridPtrVec grids;
+    grids.push_back(fgfd);
+    grids.push_back(hghd);
+    grids.push_back(fghd);
+    writeFile.write(grids);
+    writeFile.close();
+
+    openvdb::io::File readFile(fileName);
+    readFile.open();
+    openvdb::FloatGrid::Ptr readFgfd = openvdb::gridPtrCast<openvdb::FloatGrid>(readFile.readGrid("float_grid"));
+    openvdb::HalfGrid::Ptr readHghd = openvdb::gridPtrCast<openvdb::HalfGrid>(readFile.readGrid("half_grid"));
+    openvdb::FloatGrid::Ptr readFghdFloat = openvdb::gridPtrCast<openvdb::FloatGrid>(readFile.readGrid("float_grid_as_half"));
+    openvdb::HalfGrid::Ptr readFghdHalf = openvdb::gridPtrCast<openvdb::HalfGrid>(readFile.readGrid("float_grid_as_half"));
+
+    EXPECT_TRUE(readFgfd != NULL);
+    EXPECT_TRUE(readHghd != NULL);
+    EXPECT_TRUE(readFghdFloat != NULL);
+    EXPECT_TRUE(readFghdHalf == NULL);
+
+    {
+        openvdb::FloatGrid::Accessor fgfdAcc = readFgfd->getAccessor();
+        openvdb::HalfGrid::Accessor hghdAcc = readHghd->getAccessor();
+        openvdb::FloatGrid::Accessor fghdAcc = readFghdFloat->getAccessor();
+
+        EXPECT_EQ(fgfdAcc.getValue(xyz), TEST_FLOAT);
+        EXPECT_EQ(hghdAcc.getValue(xyz), math::half::posInf());
+        EXPECT_EQ(fghdAcc.getValue(xyz), math::half::posInf());
+    }
+
+    readFile.close();
+    ::remove(fileName.c_str());
+}
+
+
+////////////////////////////////////////
+
+
 TEST_F(TestGridIO, testReadAllBool) { readAllTest<openvdb::BoolGrid>(); }
 TEST_F(TestGridIO, testReadAllFloat) { readAllTest<openvdb::FloatGrid>(); }
 TEST_F(TestGridIO, testReadAllVec3S) { readAllTest<openvdb::Vec3SGrid>(); }
 TEST_F(TestGridIO, testReadAllFloat5432) { Float5432Grid::registerGrid(); readAllTest<Float5432Grid>(); }
+TEST_F(TestGridIO, testReadWriteHalfGrid) { readWriteHalfGridTest(); }