Add support and test for doing arithmetic with Vec3H, i.e. Vec3<math::half>.

apradhana · apradhana · commit 99bb2f6ed9e8 · 2024-11-05T16:30:34.000-08:00
Signed-off-by: apradhana &lt;andre.pradhana@gmail.com&gt;
diff --git a/openvdb/openvdb/math/Math.h b/openvdb/openvdb/math/Math.h
@@ -75,6 +75,16 @@ template<> inline std::string zeroVal<std::string>() { return ""; }
 /// Return the @c bool value that corresponds to zero.
 template<> inline constexpr bool zeroVal<bool>() { return false; }
 
+
+/// @note Extends the implementation of std::is_arithmetic to support math::half
+template<typename T>
+struct is_arithmetic : std::is_arithmetic<T> {};
+template<>
+struct is_arithmetic<math::half> : std::true_type {};
+// Helper variable template (equivalent to std::is_arithmetic_v)
+template<typename T>
+inline constexpr bool is_arithmetic_v = is_arithmetic<T>::value;
+
 namespace math {
 
 /// @todo These won't be needed if we eliminate StringGrids.
@@ -987,12 +997,12 @@ enum RotationOrder {
     ZXZ_ROTATION
 };
 
-
-template <typename S, typename T>
+template <typename S, typename T, typename = std::enable_if_t<openvdb::is_arithmetic_v<S>&& openvdb::is_arithmetic_v<T>>>
 struct promote {
     using type = typename boost::numeric::conversion_traits<S, T>::supertype;
 };
 
+
 /// @brief Return the index [0,1,2] of the smallest value in a 3D vector.
 /// @note This methods assumes operator[] exists.
 /// @details The return value corresponds to the largest index of the of
diff --git a/openvdb/openvdb/math/Vec3.h b/openvdb/openvdb/math/Vec3.h
@@ -663,11 +663,13 @@ using Vec3i = Vec3<int32_t>;
 using Vec3ui = Vec3<uint32_t>;
 using Vec3s = Vec3<float>;
 using Vec3d = Vec3<double>;
+using Vec3h = Vec3<math::half>;
 
 OPENVDB_IS_POD(Vec3i)
 OPENVDB_IS_POD(Vec3ui)
 OPENVDB_IS_POD(Vec3s)
 OPENVDB_IS_POD(Vec3d)
+OPENVDB_IS_POD(Vec3h)
 
 } // namespace math
 } // namespace OPENVDB_VERSION_NAME
diff --git a/openvdb/openvdb/unittest/TestMath.cc b/openvdb/openvdb/unittest/TestMath.cc
@@ -13,6 +13,185 @@ class TestMath: public ::testing::Test
 {
 };
 
+// Testing the operators within the member functions of Vec3
+template<typename ValueT>
+void testMemberOperatorsImpl()
+{
+    using namespace openvdb;
+    using Vec3T = math::Vec3<ValueT>;
+
+    {
+        Vec3T vecA(3.14,2.18,ValueT(-299792458.f));
+
+        // Alternative indexed the elements
+        EXPECT_EQ(vecA(0), ValueT(3.14));
+        EXPECT_EQ(vecA(1), ValueT(2.18));
+        EXPECT_EQ(vecA(2), ValueT(-299792458.f));
+
+        // Assignment operator
+        Vec3T vecB = vecA;
+        EXPECT_EQ(vecB(0), vecA(0));
+        EXPECT_EQ(vecB(1), vecA(1));
+        EXPECT_EQ(vecB(2), vecA(2));
+
+        // Negation operator
+        Vec3T vecC = -vecA;
+        EXPECT_EQ(vecC(0), -vecA(0));
+        EXPECT_EQ(vecC(1), -vecA(1));
+        EXPECT_EQ(vecC(2), -vecA(2));
+
+        // Multiply each element of the vector by a scalar
+        Vec3T vecD = vecA;
+        const ValueT gr = ValueT(1.6180339887);
+        vecD *= gr;
+        EXPECT_EQ(vecD(0), ValueT(gr * vecA(0)));
+        EXPECT_EQ(vecD(1), ValueT(gr * vecA(1)));
+        EXPECT_EQ(vecD(2), ValueT(gr * vecA(2)));
+
+        // Multiply each element of the vector by the corresponding element
+        Vec3T vecE = vecA;
+        Vec3T vecF(-2.5, 1.2, 3.14159);
+        vecE *= vecF;
+        EXPECT_EQ(vecE(0), ValueT(vecA(0) * vecF(0)));
+        EXPECT_EQ(vecE(1), ValueT(vecA(1) * vecF(1)));
+        EXPECT_EQ(vecE(2), ValueT(vecA(2) * vecF(2)));
+
+        // Divide each element of the vector by a scalar
+        Vec3T vecG = vecA;
+        vecG /= gr;
+        EXPECT_EQ(vecG(0), ValueT(vecA(0) / gr));
+        EXPECT_EQ(vecG(1), ValueT(vecA(1) / gr));
+        EXPECT_EQ(vecG(2), ValueT(vecA(2) / gr));
+
+        // Divide each element of the vector by the corresponding element of the given vector
+        Vec3T vecH = vecA;
+        vecH /= vecF;
+        EXPECT_EQ(vecH(0), ValueT(vecA(0) / vecF(0)));
+        EXPECT_EQ(vecH(1), ValueT(vecA(1) / vecF(1)));
+        EXPECT_EQ(vecH(2), ValueT(vecA(2) / vecF(2)));
+
+        // Add a scalar to each element of the vector
+        Vec3T vecI = vecA;
+        vecI += gr;
+        EXPECT_EQ(vecI(0), ValueT(vecA(0) + gr));
+        EXPECT_EQ(vecI(1), ValueT(vecA(1) + gr));
+        EXPECT_EQ(vecI(2), ValueT(vecA(2) + gr));
+
+        // Add each element of the given vector to the corresponding element of this vector
+        Vec3T vecJ = vecA;
+        vecJ += vecF;
+        EXPECT_EQ(vecJ(0), ValueT(vecA(0) + vecF(0)));
+        EXPECT_EQ(vecJ(1), ValueT(vecA(1) + vecF(1)));
+        EXPECT_EQ(vecJ(2), ValueT(vecA(2) + vecF(2)));
+
+        // Subtract a scalar from each element of this vector
+        Vec3T vecK = vecA;
+        vecK -= gr;
+        EXPECT_EQ(vecK(0), ValueT(vecA(0) - gr));
+        EXPECT_EQ(vecK(1), ValueT(vecA(1) - gr));
+        EXPECT_EQ(vecK(2), ValueT(vecA(2) - gr));
+
+        // Subtract each element of the given vector from the corresponding element of this vector
+        Vec3T vecL = vecA;
+        vecL -= vecF;
+        EXPECT_EQ(vecL(0), ValueT(vecA(0) - vecF(0)));
+        EXPECT_EQ(vecL(1), ValueT(vecA(1) - vecF(1)));
+        EXPECT_EQ(vecL(2), ValueT(vecA(2) - vecF(2)));
+    }
+}
+
+TEST_F(TestMath, testMemberOperators)
+{
+    using namespace openvdb;
+
+    testMemberOperatorsImpl<math::half>();
+    testMemberOperatorsImpl<float>();
+    testMemberOperatorsImpl<double>();
+}
+
+template<typename ValueT>
+void testFreeFunctionOperatorsImpl()
+{
+    using namespace openvdb;
+    using Vec3T = math::Vec3<ValueT>;
+
+    {
+        // Vec3T vecA(3.14,2.18,ValueT(-299792458.f));
+        Vec3T vecA(1,2,3);
+        Vec3T vecB(3,4,5);
+        const ValueT gr = ValueT(1.6180339887);
+
+        /// Equality operator, does exact floating point comparisons ==
+        bool eqRes = vecA == vecB;
+        EXPECT_FALSE(eqRes);
+
+        /// Inequality operator, does exact floating point comparisons !=
+        bool ineqRes = vecA != vecB;
+        EXPECT_TRUE(ineqRes);
+
+        /// Multiply each element of the given vector by @a scalar and return the result. scalar * vec
+        Vec3T scalarMultiplyA = gr * vecA;
+        EXPECT_EQ(scalarMultiplyA(0), ValueT(vecA(0) * gr));
+        EXPECT_EQ(scalarMultiplyA(1), ValueT(vecA(1) * gr));
+        EXPECT_EQ(scalarMultiplyA(2), ValueT(vecA(2) * gr));
+
+        /// Multiply each element of the given vector by @a scalar and return the result. vec * scalar
+        Vec3T scalarMultiplyB = vecA * gr;
+        EXPECT_EQ(scalarMultiplyB(0), ValueT(vecA(0) * gr));
+        EXPECT_EQ(scalarMultiplyB(1), ValueT(vecA(1) * gr));
+        EXPECT_EQ(scalarMultiplyB(2), ValueT(vecA(2) * gr));
+
+        /// Multiply corresponding elements of @a v0 and @a v1 and return the result. vec0 * vec1
+        Vec3T multiplyRes = vecA * vecB;
+        EXPECT_EQ(multiplyRes, Vec3T(3, 8, 15));
+
+        /// Divide @a scalar by each element of the given vector and return the result. a / vec
+        Vec3T scalarDivA = gr / vecA;
+        EXPECT_EQ(scalarDivA(0), ValueT(gr / vecA(0)));
+        EXPECT_EQ(scalarDivA(1), ValueT(gr / vecA(1)));
+        EXPECT_EQ(scalarDivA(2), ValueT(gr / vecA(2)));
+
+        /// Divide each element of the given vector by @a scalar and return the result. vec / scalar
+        Vec3T scalarDivB = vecA / gr;
+        EXPECT_EQ(scalarDivB(0), ValueT(vecA(0) / gr));
+        EXPECT_EQ(scalarDivB(1), ValueT(vecA(1) / gr));
+        EXPECT_EQ(scalarDivB(2), ValueT(vecA(2) / gr));
+
+        /// Divide corresponding elements of @a v0 and @a v1 and return the result. vec0 / vec1
+        Vec3T divRes = vecA / vecB;
+        EXPECT_EQ(divRes(0), ValueT(ValueT(vecA(0)) / ValueT(vecB(0))));
+        EXPECT_EQ(divRes(1), ValueT(ValueT(vecA(1)) / ValueT(vecB(1))));
+        EXPECT_EQ(divRes(2), ValueT(ValueT(vecA(2)) / ValueT(vecB(2))));
+
+        /// Add corresponding elements of @a v0 and @a v1 and return the result. vec0 + vec1
+        Vec3T addRes = vecA + vecB;
+        EXPECT_EQ(addRes, Vec3T(4, 6, 8));
+
+        /// Add @a scalar to each element of the given vector and return the result. a + vec
+        Vec3T addScalarRes = vecA + gr;
+        EXPECT_EQ(addScalarRes(0), ValueT(vecA(0) + gr));
+        EXPECT_EQ(addScalarRes(1), ValueT(vecA(1) + gr));
+        EXPECT_EQ(addScalarRes(2), ValueT(vecA(2) + gr));
+
+        /// Subtract corresponding elements of @a v0 and @a v1 and return the result. vec0 - vec1
+        Vec3T subtractRes = vecA - vecB;
+        EXPECT_EQ(subtractRes, Vec3T(-2, -2, -2));
+
+        /// Subtract @a scalar from each element of the given vector and return the result. vec0 - a
+        Vec3T subScalarRes = vecA - gr;
+        EXPECT_EQ(subScalarRes(0), ValueT(vecA(0) - gr));
+        EXPECT_EQ(subScalarRes(1), ValueT(vecA(1) - gr));
+        EXPECT_EQ(subScalarRes(2), ValueT(vecA(2) - gr));
+    }
+}
+
+TEST_F(TestMath, testFreeFunctionOperators)
+{
+    using namespace openvdb;
+    testFreeFunctionOperatorsImpl<math::half>();
+
+}
+
 
 // This suite of tests obviously needs to be expanded!
 TEST_F(TestMath, testAll)
@@ -23,6 +202,7 @@ TEST_F(TestMath, testAll)
         EXPECT_EQ(math::Sign( 3   ), 1);
         EXPECT_EQ(math::Sign(-1.0 ),-1);
         EXPECT_EQ(math::Sign( 0.0f), 0);
+        EXPECT_EQ(math::Sign(math::half(0.)), 0);
     }
     {// SignChange
         EXPECT_TRUE( math::SignChange( -1, 1));
diff --git a/openvdb/openvdb/unittest/TestTools.cc b/openvdb/openvdb/unittest/TestTools.cc
@@ -123,7 +123,7 @@ testLevelSetSphereImpl()
     using Vec3T = typename openvdb::math::Vec3<ValueT>;
 
     const ValueT radius = 4.3f;
-    const Vec3T center(ValueT(15.8), ValueT(13.2), ValueT(16.7));
+    Vec3T center(ValueT(15.8), ValueT(13.2), ValueT(16.7));
     const ValueT voxelSize = 1.5f, width = 3.25f;
     const int dim = 32;
     ValueT tolerance = std::is_floating_point<typename GridT::ValueType>::value ? 0.0001 : 0.004;
@@ -140,8 +140,14 @@ testLevelSetSphereImpl()
     for (int i=0; i<dim; ++i) {
         for (int j=0; j<dim; ++j) {
             for (int k=0; k<dim; ++k) {
+                // Vec3T p(voxelSize*float(i), voxelSize*float(j), voxelSize*float(k));
+                // Vec3T foobar = p;
+                // foobar -= center;
+                // const float dist = foobar.length() - radius;
                 const openvdb::Vec3f p(voxelSize*float(i), voxelSize*float(j), voxelSize*float(k));
-                const float dist = (p-center).length() - radius;
+                const float dist = (p-openvdb::Vec3f(center)).length() - radius;
+                //char* foo = openvdb::math::promote<openvdb::Vec3f::ValueType, typename Vec3T::ValueType>::type(1);
+                //std::cout << foo << std::endl;
                 const float val1 = grid1->tree().getValue(openvdb::Coord(i,j,k));
                 const float val2 = grid2->tree().getValue(openvdb::Coord(i,j,k));
                 if (dist > outside) {
diff --git a/openvdb/openvdb/unittest/TestTree.cc b/openvdb/openvdb/unittest/TestTree.cc
@@ -2015,7 +2015,9 @@ testSignedFloodFillImpl()
         for (xyz[1]=0; xyz[1]<dim[1]; ++xyz[1]) {
             for (xyz[2]=0; xyz[2]<dim[2]; ++xyz[2]) {
                 const openvdb::Vec3R p =  grid->transform().indexToWorld(xyz);
-                const ValueT dist = float((p-center).length() - radius);
+                Vec3T foobar = p;
+                foobar -= center;
+                const ValueT dist = float(foobar.length() - radius);
                 if (fabs(dist) > outside) continue;
                 acc.setValue(xyz, dist);
             }
@@ -2030,7 +2032,9 @@ testSignedFloodFillImpl()
         for (xyz[1]=0; xyz[1]<dim[1]; ++xyz[1]) {
             for (xyz[2]=0; xyz[2]<dim[2]; ++xyz[2]) {
                 const openvdb::Vec3R p =  grid->transform().indexToWorld(xyz);
-                const ValueT dist = ValueT((p-center).length() - radius);
+                Vec3T foobar = p;
+                foobar -= center;
+                const ValueT dist = ValueT(foobar.length() - radius);
                 const ValueT val  =  acc.getValue(xyz);
                 if (dist < inside) {
                     ASSERT_DOUBLES_EXACTLY_EQUAL( val, outside);
@@ -2052,7 +2056,9 @@ testSignedFloodFillImpl()
         for (xyz[1]=0; xyz[1]<dim[1]; ++xyz[1]) {
             for (xyz[2]=0; xyz[2]<dim[2]; ++xyz[2]) {
                 const openvdb::Vec3R p =  grid->transform().indexToWorld(xyz);
-                const ValueT dist = ValueT((p-center).length() - radius);
+                Vec3T foobar = p;
+                foobar -= center;
+                const ValueT dist = ValueT(foobar.length() - radius);
                 const ValueT val  =  acc.getValue(xyz);
                 if (dist < inside) {
                     ASSERT_DOUBLES_EXACTLY_EQUAL( val, inside);
diff --git a/openvdb/openvdb/unittest/TestVolumeToSpheres.cc b/openvdb/openvdb/unittest/TestVolumeToSpheres.cc
@@ -174,7 +174,7 @@ testClosestSurfacePointImpl()
     using Vec3T = typename openvdb::math::Vec3<ValueT>;
 
     const ValueT voxelSize = ValueT(1.0);
-    const Vec3T center{ValueT(0.0)}; // ensure multiple internal nodes
+    const Vec3R center{ValueT(0.0)}; // ensure multiple internal nodes
 
     for (const float radius: { 8.0f, 50.0f }) {
         // Construct a spherical level set.
diff --git a/pendingchanges/half_grid_support.txt b/pendingchanges/half_grid_support.txt
@@ -21,4 +21,4 @@ Tests:
 [ ] Test explicit template instantiation
 
 # TODO, but will be addressed later:
-[ ] Add support for csgUnion for HalfGrid.
+[ ] Add support for csgUnion for HalfGrid.
