diff --git a/openvdb/openvdb/Grid.h b/openvdb/openvdb/Grid.h
index 2823ff3dc6..f90f81d056 100644
--- a/openvdb/openvdb/Grid.h
+++ b/openvdb/openvdb/Grid.h
@@ -1750,7 +1750,7 @@ createLevelSet(Real voxelSize, Real halfWidth)
     using ValueType = typename GridType::ValueType;
 
     // GridType::ValueType is required to be a floating-point scalar.
-    static_assert(std::is_floating_point<ValueType>::value,
+    static_assert(openvdb::is_floating_point<ValueType>::value,
         "level-set grids must be floating-point-valued");
 
     typename GridType::Ptr grid = GridType::create(
diff --git a/openvdb/openvdb/Types.h b/openvdb/openvdb/Types.h
index cfba40476a..eb9a604f91 100644
--- a/openvdb/openvdb/Types.h
+++ b/openvdb/openvdb/Types.h
@@ -8,27 +8,7 @@
 #include "Platform.h"
 #include "TypeList.h" // backwards compat
 
-#ifdef OPENVDB_USE_IMATH_HALF
-#ifdef OPENVDB_IMATH_VERSION
-#include <Imath/half.h>
-#else
-#include <OpenEXR/half.h>
-#endif
-namespace openvdb {
-OPENVDB_USE_VERSION_NAMESPACE
-namespace OPENVDB_VERSION_NAME {
-namespace math {
-using half = half;
-}}}
-#else
-#include <openvdb/math/Half.h>
-namespace openvdb {
-OPENVDB_USE_VERSION_NAMESPACE
-namespace OPENVDB_VERSION_NAME {
-namespace math {
-using half = internal::half;
-}}}
-#endif
+#include <openvdb/math/HalfDecl.h>
 
 #include <openvdb/math/Math.h>
 #include <openvdb/math/BBox.h>
@@ -454,6 +434,14 @@ struct is_floating_point : std::is_floating_point<T> { };
 template<>
 struct is_floating_point<Half> : std::is_floating_point<float> { };
 
+
+template<class T>
+struct is_signed : std::is_signed<T> { };
+
+template<>
+struct is_signed<Half> : std::is_signed<float> { };
+
+
 /// @brief Maps one type (e.g. the value types) to other types
 template<typename T>
 struct ValueToComputeMap
diff --git a/openvdb/openvdb/math/HalfDecl.h b/openvdb/openvdb/math/HalfDecl.h
new file mode 100644
index 0000000000..f9d654557f
--- /dev/null
+++ b/openvdb/openvdb/math/HalfDecl.h
@@ -0,0 +1,30 @@
+// Copyright Contributors to the OpenVDB Project
+// SPDX-License-Identifier: MPL-2.0
+
+#ifndef OPENVDB_HALFDECL_HAS_BEEN_INCLUDED
+#define OPENVDB_HALFDECL_HAS_BEEN_INCLUDED
+
+#ifdef OPENVDB_USE_IMATH_HALF
+#ifdef OPENVDB_IMATH_VERSION
+#include <Imath/half.h>
+#else
+#include <OpenEXR/half.h>
+#endif
+namespace openvdb {
+OPENVDB_USE_VERSION_NAMESPACE
+namespace OPENVDB_VERSION_NAME {
+namespace math {
+using half = half;
+}}}
+#else
+#include <openvdb/math/Half.h>
+namespace openvdb {
+OPENVDB_USE_VERSION_NAMESPACE
+namespace OPENVDB_VERSION_NAME {
+namespace math {
+using half = internal::half;
+}}}
+#endif
+
+
+#endif // OPENVDB_HALFDECL_HAS_BEEN_INCLUDED
diff --git a/openvdb/openvdb/math/Math.h b/openvdb/openvdb/math/Math.h
index 7593118c2c..2138f0b7cb 100644
--- a/openvdb/openvdb/math/Math.h
+++ b/openvdb/openvdb/math/Math.h
@@ -10,6 +10,7 @@
 
 #include <openvdb/Platform.h>
 #include <openvdb/version.h>
+#include <openvdb/math/HalfDecl.h>
 #include <openvdb/util/Assert.h>
 #include <boost/numeric/conversion/conversion_traits.hpp>
 #include <algorithm> // for std::max()
@@ -146,9 +147,10 @@ template<> inline std::string negative(const std::string& val) { return val; }
 
 //@{
 /// Tolerance for floating-point comparison
-template<typename T> struct Tolerance { static T value() { return zeroVal<T>(); } };
-template<> struct Tolerance<float>    { static float value() { return 1e-8f; } };
-template<> struct Tolerance<double>   { static double value() { return 1e-15; } };
+template<typename T> struct Tolerance   { static T value() { return zeroVal<T>(); } };
+template<> struct Tolerance<math::half> { static math::half value() { return 0.00097656; } };
+template<> struct Tolerance<float>      { static float value() { return 1e-8f; } };
+template<> struct Tolerance<double>     { static double value() { return 1e-15; } };
 //@}
 
 //@{
@@ -362,6 +364,10 @@ isApproxZero(const Type& x, const Type& tolerance)
 }
 
 
+/// Return @c true if @a x is less than zero.
+inline bool
+isNegative(math::half& x) { return x.isNegative(); }
+
 /// Return @c true if @a x is less than zero.
 template<typename Type>
 inline bool
@@ -375,6 +381,10 @@ template<> inline bool isNegative<bool>(const bool&) { return false; }
 inline bool
 isFinite(const float x) { return std::isfinite(x); }
 
+/// Return @c true if @a x is finite.
+inline bool
+isFinite(const math::half x) { return x.isFinite(); }
+
 /// Return @c true if @a x is finite.
 template<typename Type, typename std::enable_if<std::is_arithmetic<Type>::value, int>::type = 0>
 inline bool
@@ -385,6 +395,10 @@ isFinite(const Type& x) { return std::isfinite(static_cast<double>(x)); }
 inline bool
 isInfinite(const float x) { return std::isinf(x); }
 
+/// Return @c true if @a x is an infinity value (either positive infinity or negative infinity).
+inline bool
+isInfinite(const math::half x) { return x.isInfinity(); }
+
 /// Return @c true if @a x is an infinity value (either positive infinity or negative infinity).
 template<typename Type, typename std::enable_if<std::is_arithmetic<Type>::value, int>::type = 0>
 inline bool
@@ -395,6 +409,10 @@ isInfinite(const Type& x) { return std::isinf(static_cast<double>(x)); }
 inline bool
 isNan(const float x) { return std::isnan(x); }
 
+/// Return @c true if @a x is a NaN (Not-A-Number) value.
+inline bool
+isNan(const math::half x) { return x.isNan(); }
+
 /// Return @c true if @a x is a NaN (Not-A-Number) value.
 template<typename Type, typename std::enable_if<std::is_arithmetic<Type>::value, int>::type = 0>
 inline bool
@@ -570,6 +588,15 @@ Pow(Type x, int n)
     return ans;
 }
 
+//@{
+/// Return @a b<sup>e</sup>.
+inline math::half
+Pow(math::half b, math::half e)
+{
+    OPENVDB_ASSERT( b >= 0.0f && "Pow(half,half): base is negative" );
+    return math::half(powf(float(b),float(e)));
+}
+
 //@{
 /// Return @a b<sup>e</sup>.
 inline float
@@ -590,12 +617,29 @@ Pow(double b, double e)
 
 // ==========> Max <==================
 
+namespace internal {
+
+inline const math::half&
+max_impl(const math::half& a, const math::half& b)
+{
+    return a > b ? a : b;
+}
+
+template<typename Type>
+inline const Type&
+max_impl(const Type& a, const Type& b)
+{
+    return std::max(a,b);
+}
+
+} // namespace internal
+
 /// Return the maximum of two values
 template<typename Type>
 inline const Type&
 Max(const Type& a, const Type& b)
 {
-    return std::max(a,b);
+    return internal::max_impl(a,b);
 }
 
 /// Return the maximum of three values
@@ -603,7 +647,7 @@ template<typename Type>
 inline const Type&
 Max(const Type& a, const Type& b, const Type& c)
 {
-    return std::max(std::max(a,b), c);
+    return internal::max_impl(internal::max_impl(a,b), c);
 }
 
 /// Return the maximum of four values
@@ -611,7 +655,7 @@ template<typename Type>
 inline const Type&
 Max(const Type& a, const Type& b, const Type& c, const Type& d)
 {
-    return std::max(std::max(a,b), std::max(c,d));
+    return internal::max_impl(internal::max_impl(a,b), internal::max_impl(c,d));
 }
 
 /// Return the maximum of five values
@@ -619,7 +663,7 @@ template<typename Type>
 inline const Type&
 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
 {
-    return std::max(std::max(a,b), Max(c,d,e));
+    return internal::max_impl(internal::max_impl(a,b), Max(c,d,e));
 }
 
 /// Return the maximum of six values
@@ -627,7 +671,7 @@ template<typename Type>
 inline const Type&
 Max(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
 {
-    return std::max(Max(a,b,c), Max(d,e,f));
+    return internal::max_impl(Max(a,b,c), Max(d,e,f));
 }
 
 /// Return the maximum of seven values
@@ -636,7 +680,7 @@ inline const Type&
 Max(const Type& a, const Type& b, const Type& c, const Type& d,
     const Type& e, const Type& f, const Type& g)
 {
-    return std::max(Max(a,b,c,d), Max(e,f,g));
+    return internal::max_impl(Max(a,b,c,d), Max(e,f,g));
 }
 
 /// Return the maximum of eight values
@@ -645,28 +689,48 @@ inline const Type&
 Max(const Type& a, const Type& b, const Type& c, const Type& d,
     const Type& e, const Type& f, const Type& g, const Type& h)
 {
-    return std::max(Max(a,b,c,d), Max(e,f,g,h));
+    return internal::max_impl(Max(a,b,c,d), Max(e,f,g,h));
 }
 
 
 // ==========> Min <==================
 
+namespace internal {
+
+inline const math::half&
+min_impl(const math::half& a, const math::half& b)
+{
+    return a < b ? a : b;
+}
+
+template<typename Type>
+inline const Type&
+min_impl(const Type& a, const Type& b)
+{
+    return std::min(a,b);
+}
+
+} // namespace internal
+
 /// Return the minimum of two values
 template<typename Type>
 inline const Type&
-Min(const Type& a, const Type& b) { return std::min(a, b); }
+Min(const Type& a, const Type& b) { return internal::min_impl(a, b); }
 
 /// Return the minimum of three values
 template<typename Type>
 inline const Type&
-Min(const Type& a, const Type& b, const Type& c) { return std::min(std::min(a, b), c); }
+Min(const Type& a, const Type& b, const Type& c)
+{
+    return internal::min_impl(internal::min_impl(a, b), c);
+}
 
 /// Return the minimum of four values
 template<typename Type>
 inline const Type&
 Min(const Type& a, const Type& b, const Type& c, const Type& d)
 {
-    return std::min(std::min(a, b), std::min(c, d));
+    return internal::min_impl(internal::min_impl(a, b), internal::min_impl(c, d));
 }
 
 /// Return the minimum of five values
@@ -674,7 +738,7 @@ template<typename Type>
 inline const Type&
 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e)
 {
-    return std::min(std::min(a,b), Min(c,d,e));
+    return internal::min_impl(internal::min_impl(a,b), Min(c,d,e));
 }
 
 /// Return the minimum of six values
@@ -682,7 +746,7 @@ template<typename Type>
 inline const Type&
 Min(const Type& a, const Type& b, const Type& c, const Type& d, const Type& e, const Type& f)
 {
-    return std::min(Min(a,b,c), Min(d,e,f));
+    return internal::min_impl(Min(a,b,c), Min(d,e,f));
 }
 
 /// Return the minimum of seven values
@@ -691,7 +755,7 @@ inline const Type&
 Min(const Type& a, const Type& b, const Type& c, const Type& d,
     const Type& e, const Type& f, const Type& g)
 {
-    return std::min(Min(a,b,c,d), Min(e,f,g));
+    return internal::min_impl(Min(a,b,c,d), Min(e,f,g));
 }
 
 /// Return the minimum of eight values
@@ -700,7 +764,7 @@ inline const Type&
 Min(const Type& a, const Type& b, const Type& c, const Type& d,
     const Type& e, const Type& f, const Type& g, const Type& h)
 {
-    return std::min(Min(a,b,c,d), Min(e,f,g,h));
+    return internal::min_impl(Min(a,b,c,d), Min(e,f,g,h));
 }
 
 
diff --git a/openvdb/openvdb/math/Stencils.h b/openvdb/openvdb/math/Stencils.h
index 58d96d51d4..7688843429 100644
--- a/openvdb/openvdb/math/Stencils.h
+++ b/openvdb/openvdb/math/Stencils.h
@@ -1545,7 +1545,7 @@ class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, Grid
     {
         Real alpha, normGrad;
         return this->meanCurvature(alpha, normGrad) ?
-               ValueType(alpha*mInv2Dx/math::Pow3(normGrad)) : 0;
+               ValueType(alpha*mInv2Dx/math::Pow3(normGrad)) : ValueType(0);
     }
 
     /// @brief Return the Gaussian curvature at the previously buffered location.
@@ -1556,7 +1556,7 @@ class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, Grid
     {
         Real alpha, normGrad;
         return this->gaussianCurvature(alpha, normGrad) ?
-               ValueType(alpha*mInvDx2/math::Pow4(normGrad)) : 0;
+               ValueType(alpha*mInvDx2/math::Pow4(normGrad)) : ValueType(0);
     }
 
     /// @brief Return both the mean and the Gaussian curvature at the
@@ -1571,7 +1571,7 @@ class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, Grid
           mean  = ValueType(alphaM*mInv2Dx/math::Pow3(normGrad));
           gauss = ValueType(alphaG*mInvDx2/math::Pow4(normGrad));
         } else {
-          mean = gauss = 0;
+          mean = gauss = ValueType(0);
         }
     }
 
@@ -1585,7 +1585,7 @@ class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, Grid
     {
         Real alpha, normGrad;
         return this->meanCurvature(alpha, normGrad) ?
-               ValueType(alpha*mInvDx2/(2*math::Pow2(normGrad))) : 0;
+               ValueType(alpha*mInvDx2/(2*math::Pow2(normGrad))) : ValueType(0);
     }
 
     /// Return the mean Gaussian multiplied by the norm of the
@@ -1597,7 +1597,7 @@ class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, Grid
     {
         Real alpha, normGrad;
         return this->gaussianCurvature(alpha, normGrad) ?
-               ValueType(2*alpha*mInv2Dx*mInvDx2/math::Pow3(normGrad)) : 0;
+               ValueType(2*alpha*mInv2Dx*mInvDx2/math::Pow3(normGrad)) : ValueType(0);
     }
 
     /// @brief Return both the mean and the Gaussian curvature at the
@@ -1627,7 +1627,7 @@ class CurvatureStencil: public BaseStencil<CurvatureStencil<GridT, IsSafe>, Grid
         Real alphaM, alphaG, normGrad;
         if (this->curvatures(alphaM, alphaG, normGrad)) {
           const Real mean = alphaM*mInv2Dx/math::Pow3(normGrad);
-          const Real tmp = std::sqrt(mean*mean - alphaG*mInvDx2/math::Pow4(normGrad));
+          const Real tmp = math::Sqrt(mean*mean - alphaG*mInvDx2/math::Pow4(normGrad));
           pair.first  = ValueType(mean - tmp);
           pair.second = ValueType(mean + tmp);
         }
diff --git a/openvdb/openvdb/tools/FastSweeping.h b/openvdb/openvdb/tools/FastSweeping.h
index c457700f77..ee74491733 100644
--- a/openvdb/openvdb/tools/FastSweeping.h
+++ b/openvdb/openvdb/tools/FastSweeping.h
@@ -30,6 +30,7 @@
 //#define BENCHMARK_FAST_SWEEPING
 
 #include <openvdb/openvdb.h>
+#include <openvdb/Types.h>
 #include <openvdb/Platform.h>
 #include <openvdb/math/Math.h> // for Abs() and isExactlyEqual()
 #include <openvdb/math/Stencils.h> // for GradStencil
@@ -460,7 +461,7 @@ maskSdf(const GridT &sdfGrid,
 template<typename SdfGridT, typename ExtValueT = typename SdfGridT::ValueType>
 class FastSweeping
 {
-    static_assert(std::is_floating_point<typename SdfGridT::ValueType>::value,
+    static_assert(openvdb::is_floating_point<typename SdfGridT::ValueType>::value,
                   "FastSweeping requires SdfGridT to have floating-point values");
     // Defined types related to the signed distance (or fog) grid
     using SdfValueT = typename SdfGridT::ValueType;
diff --git a/openvdb/openvdb/tools/LevelSetFilter.h b/openvdb/openvdb/tools/LevelSetFilter.h
index a31b7344be..5011da8b26 100644
--- a/openvdb/openvdb/tools/LevelSetFilter.h
+++ b/openvdb/openvdb/tools/LevelSetFilter.h
@@ -45,7 +45,7 @@ class LevelSetFilter : public LevelSetTracker<GridT, InterruptT>
     using TreeType = typename GridType::TreeType;
     using ValueType = typename TreeType::ValueType;
     using AlphaType = typename MaskType::ValueType;
-    static_assert(std::is_floating_point<AlphaType>::value,
+    static_assert(openvdb::is_floating_point<AlphaType>::value,
         "LevelSetFilter requires a mask grid with floating-point values");
 
     /// @brief Main constructor from a grid
diff --git a/openvdb/openvdb/tools/LevelSetMeasure.h b/openvdb/openvdb/tools/LevelSetMeasure.h
index a4ad847704..c658164739 100644
--- a/openvdb/openvdb/tools/LevelSetMeasure.h
+++ b/openvdb/openvdb/tools/LevelSetMeasure.h
@@ -108,7 +108,7 @@ class LevelSetMeasure
     using ValueType = typename TreeType::ValueType;
     using ManagerType = typename tree::LeafManager<const TreeType>;
 
-    static_assert(std::is_floating_point<ValueType>::value,
+    static_assert(openvdb::is_floating_point<ValueType>::value,
         "level set measure is supported only for scalar, floating-point grids");
 
     /// @brief Main constructor from a grid
@@ -408,7 +408,7 @@ MeasureCurvatures::operator()(const LeafRange& range) const
 
 template<class GridT>
 inline
-typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, Real>::type
+typename std::enable_if<openvdb::is_floating_point<typename GridT::ValueType>::value, Real>::type
 doLevelSetArea(const GridT& grid, bool useWorldUnits)
 {
     LevelSetMeasure<GridT> m(grid);
@@ -417,7 +417,7 @@ doLevelSetArea(const GridT& grid, bool useWorldUnits)
 
 template<class GridT>
 inline
-typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, Real>::type
+typename std::enable_if<!openvdb::is_floating_point<typename GridT::ValueType>::value, Real>::type
 doLevelSetArea(const GridT&, bool)
 {
     OPENVDB_THROW(TypeError,
@@ -441,7 +441,7 @@ levelSetArea(const GridT& grid, bool useWorldUnits)
 
 template<class GridT>
 inline
-typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, Real>::type
+typename std::enable_if<openvdb::is_floating_point<typename GridT::ValueType>::value, Real>::type
 doLevelSetVolume(const GridT& grid, bool useWorldUnits)
 {
     LevelSetMeasure<GridT> m(grid);
@@ -450,7 +450,7 @@ doLevelSetVolume(const GridT& grid, bool useWorldUnits)
 
 template<class GridT>
 inline
-typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, Real>::type
+typename std::enable_if<!openvdb::is_floating_point<typename GridT::ValueType>::value, Real>::type
 doLevelSetVolume(const GridT&, bool)
 {
     OPENVDB_THROW(TypeError,
@@ -474,7 +474,7 @@ levelSetVolume(const GridT& grid, bool useWorldUnits)
 
 template<class GridT>
 inline
-typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, int>::type
+typename std::enable_if<openvdb::is_floating_point<typename GridT::ValueType>::value, int>::type
 doLevelSetEulerCharacteristic(const GridT& grid)
 {
     LevelSetMeasure<GridT> m(grid);
@@ -483,7 +483,7 @@ doLevelSetEulerCharacteristic(const GridT& grid)
 
 template<class GridT>
 inline
-typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, int>::type
+typename std::enable_if<!openvdb::is_floating_point<typename GridT::ValueType>::value, int>::type
 doLevelSetEulerCharacteristic(const GridT&)
 {
     OPENVDB_THROW(TypeError,
@@ -508,7 +508,7 @@ levelSetEulerCharacteristic(const GridT& grid)
 
 template<class GridT>
 inline
-typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, int>::type
+typename std::enable_if<openvdb::is_floating_point<typename GridT::ValueType>::value, int>::type
 doLevelSetEuler(const GridT& grid)
 {
     LevelSetMeasure<GridT> m(grid);
@@ -518,7 +518,7 @@ doLevelSetEuler(const GridT& grid)
 
 template<class GridT>
 inline
-typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, int>::type
+typename std::enable_if<openvdb::is_floating_point<typename GridT::ValueType>::value, int>::type
 doLevelSetGenus(const GridT& grid)
 {
     LevelSetMeasure<GridT> m(grid);
@@ -527,7 +527,7 @@ doLevelSetGenus(const GridT& grid)
 
 template<class GridT>
 inline
-typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, int>::type
+typename std::enable_if<!openvdb::is_floating_point<typename GridT::ValueType>::value, int>::type
 doLevelSetGenus(const GridT&)
 {
     OPENVDB_THROW(TypeError,
diff --git a/openvdb/openvdb/tools/LevelSetMorph.h b/openvdb/openvdb/tools/LevelSetMorph.h
index 6812999bbf..5faee232ec 100644
--- a/openvdb/openvdb/tools/LevelSetMorph.h
+++ b/openvdb/openvdb/tools/LevelSetMorph.h
@@ -498,7 +498,7 @@ sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
                 ValueType& s = speed[voxelIter.pos()];
                 s -= target.wsSample(map.applyMap(voxelIter.getCoord().asVec3d()));
                 if (!math::isApproxZero(s)) isZero = false;
-                mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
+                mMaxAbsS = math::Max(mMaxAbsS, ValueType(math::Abs(s)));
             }
             if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
         }
@@ -515,9 +515,9 @@ sampleXformedSpeed(const LeafRange& range, Index speedBuffer)
                 const ValueType a = math::SmoothUnitStep((mask.wsSample(xyz)-min)*invNorm);
                 ValueType& s = speed[voxelIter.pos()];
                 s -= target.wsSample(xyz);
-                s *= invMask ? 1 - a : a;
+                s *= invMask ? ValueType(1 - a) : a;
                 if (!math::isApproxZero(s)) isZero = false;
-                mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
+                mMaxAbsS = math::Max(mMaxAbsS, ValueType(math::Abs(s)));
             }
             if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
         }
@@ -546,7 +546,7 @@ sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
                 ValueType& s = speed[voxelIter.pos()];
                 s -= target.getValue(voxelIter.getCoord());
                 if (!math::isApproxZero(s)) isZero = false;
-                mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
+                mMaxAbsS = math::Max(mMaxAbsS, ValueType(math::Abs(s)));
             }
             if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
         }
@@ -562,9 +562,9 @@ sampleAlignedSpeed(const LeafRange& range, Index speedBuffer)
                 const ValueType a = math::SmoothUnitStep((mask.getValue(ijk)-min)*invNorm);
                 ValueType& s = speed[voxelIter.pos()];
                 s -= target.getValue(ijk);
-                s *= invMask ? 1 - a : a;
+                s *= invMask ? ValueType(1 - a) : a;
                 if (!math::isApproxZero(s)) isZero = false;
-                mMaxAbsS = math::Max(mMaxAbsS, math::Abs(s));
+                mMaxAbsS = math::Max(mMaxAbsS, ValueType(math::Abs(s)));
             }
             if (isZero) speed[0] = std::numeric_limits<ValueType>::max();//tag first voxel
         }
@@ -632,7 +632,7 @@ euler(const LeafRange& range, ValueType dt,
             if (math::isApproxZero(speed[n])) continue;
             stencil.moveTo(voxelIter);
             const ValueType v = stencil.getValue() - dt * speed[n] * NumGrad::result(map, stencil);
-            result[n] = Nominator ? Alpha * phi[n] + Beta * v : v;
+            result[n] = Nominator ? ValueType(Alpha * phi[n] + Beta * v) : v;
         }//loop over active voxels in the leaf of the mask
     }//loop over leafs of the level set
 }
diff --git a/openvdb/openvdb/tools/LevelSetSphere.h b/openvdb/openvdb/tools/LevelSetSphere.h
index 86c7bbe950..5014c8b4f7 100644
--- a/openvdb/openvdb/tools/LevelSetSphere.h
+++ b/openvdb/openvdb/tools/LevelSetSphere.h
@@ -90,8 +90,9 @@ class LevelSetSphere
 public:
     using TreeT  = typename GridT::TreeType;
     using ValueT = typename GridT::ValueType;
-    using Vec3T  = typename math::Vec3<ValueT>;
-    static_assert(std::is_floating_point<ValueT>::value,
+    using ComputeT = typename GridT::ComputeType;
+    using Vec3T  = typename math::Vec3<ComputeT>;
+    static_assert(openvdb::is_floating_point<ValueT>::value,
         "level set grids must have scalar, floating-point value types");
 
     /// @brief Constructor
@@ -104,7 +105,7 @@ class LevelSetSphere
     /// @note If the radius of the sphere is smaller than
     /// 1.5*voxelSize, i.e. the sphere is smaller than the Nyquist
     /// frequency of the grid, it is ignored!
-    LevelSetSphere(ValueT radius, const Vec3T &center, InterruptT* interrupt = nullptr)
+    LevelSetSphere(ComputeT radius, const Vec3T &center, InterruptT* interrupt = nullptr)
         : mRadius(radius), mCenter(center), mInterrupt(interrupt)
     {
         if (mRadius<=0) OPENVDB_THROW(ValueError, "radius must be positive");
@@ -115,7 +116,7 @@ class LevelSetSphere
     /// @param voxelSize  Size of voxels in world units
     /// @param halfWidth  Half-width of narrow-band in voxel units
     /// @param threaded   If true multi-threading is enabled (true by default)
-    typename GridT::Ptr getLevelSet(ValueT voxelSize, ValueT halfWidth, bool threaded = true)
+    typename GridT::Ptr getLevelSet(ComputeT voxelSize, ComputeT halfWidth, bool threaded = true)
     {
         mGrid = createLevelSet<GridT>(voxelSize, halfWidth);
         this->rasterSphere(voxelSize, halfWidth, threaded);
@@ -124,7 +125,7 @@ class LevelSetSphere
     }
 
 private:
-    void rasterSphere(ValueT dx, ValueT w, bool threaded)
+    void rasterSphere(ComputeT dx, ComputeT w, bool threaded)
     {
         if (!(dx>0.0f)) OPENVDB_THROW(ValueError, "voxel size must be positive");
         if (!(w>1)) OPENVDB_THROW(ValueError, "half-width must be larger than one");
@@ -158,16 +159,16 @@ class LevelSetSphere
             // Compute signed distances to sphere using leapfrogging in k
             for (i = r.begin(); i != r.end(); ++i) {
                 if (util::wasInterrupted(mInterrupt)) return;
-                const auto x2 = math::Pow2(ValueT(i) - c[0]);
+                const auto x2 = math::Pow2(ComputeT(i) - c[0]);
                 for (j = jmin; j <= jmax; ++j) {
-                    const auto x2y2 = math::Pow2(ValueT(j) - c[1]) + x2;
+                    const auto x2y2 = math::Pow2(ComputeT(j) - c[1]) + x2;
                     for (k = kmin; k <= kmax; k += m) {
                         m = 1;
                         // Distance in voxel units to sphere
-                        const auto v = math::Sqrt(x2y2 + math::Pow2(ValueT(k)-c[2]))-r0;
+                        const auto v = math::Sqrt(x2y2 + math::Pow2(ComputeT(k)-c[2]))-r0;
                         const auto d = math::Abs(v);
                         if (d < w) { // inside narrow band
-                            acc.setValue(ijk, dx*v);// distance in world units
+                            acc.setValue(ijk, ValueT(dx*v));// distance in world units
                         } else { // outside narrow band
                             m += math::Floor(d-w);// leapfrog
                         }
@@ -221,12 +222,12 @@ createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize
     float halfWidth, InterruptT* interrupt, bool threaded)
 {
     // GridType::ValueType is required to be a floating-point scalar.
-    static_assert(std::is_floating_point<typename GridType::ValueType>::value,
+    static_assert(openvdb::is_floating_point<typename GridType::ValueType>::value,
         "level set grids must have scalar, floating-point value types");
 
-    using ValueT = typename GridType::ValueType;
-    LevelSetSphere<GridType, InterruptT> factory(ValueT(radius), center, interrupt);
-    return factory.getLevelSet(ValueT(voxelSize), ValueT(halfWidth), threaded);
+    using ComputeT = typename GridType::ComputeType;
+    LevelSetSphere<GridType, InterruptT> factory(ComputeT(radius), center, interrupt);
+    return factory.getLevelSet(ComputeT(voxelSize), ComputeT(halfWidth), threaded);
 }
 
 
diff --git a/openvdb/openvdb/tools/LevelSetTracker.h b/openvdb/openvdb/tools/LevelSetTracker.h
index 37bf176172..bdf881f248 100644
--- a/openvdb/openvdb/tools/LevelSetTracker.h
+++ b/openvdb/openvdb/tools/LevelSetTracker.h
@@ -67,7 +67,7 @@ class LevelSetTracker
     using LeafRange = typename LeafManagerType::LeafRange;
     using BufferType = typename LeafManagerType::BufferType;
     using MaskTreeType = typename TreeType::template ValueConverter<ValueMask>::Type;
-    static_assert(std::is_floating_point<ValueType>::value,
+    static_assert(openvdb::is_floating_point<ValueType>::value,
         "LevelSetTracker requires a level set grid with floating-point values");
 
     /// Lightweight struct that stores the state of the LevelSetTracker
@@ -638,8 +638,8 @@ eval(StencilT& stencil, const ValueType* phi, ValueType* result, Index n) const
     const ValueType phi0 = stencil.getValue();
     ValueType v = phi0 / ( math::Sqrt(math::Pow2(phi0) + normSqGradPhi) +
                            math::Tolerance<ValueType>::value() );
-    v = phi0 - mDt * v * (math::Sqrt(normSqGradPhi) * mInvDx - 1.0f);
-    result[n] = Nominator ? alpha * phi[n] + beta * v : v;
+    v = phi0 - mDt * v * (math::Sqrt(normSqGradPhi) * mInvDx - ValueType(1));
+    result[n] = Nominator ? ValueType(alpha * phi[n] + beta * v) : v;
 }
 
 template<typename GridT, typename InterruptT>
diff --git a/openvdb/openvdb/tools/LevelSetUtil.h b/openvdb/openvdb/tools/LevelSetUtil.h
index fd74c8ce12..8305b82cbf 100644
--- a/openvdb/openvdb/tools/LevelSetUtil.h
+++ b/openvdb/openvdb/tools/LevelSetUtil.h
@@ -461,7 +461,8 @@ struct SDFVoxelsToFogVolume {
 
             ValueType* values = node.buffer().data();
             for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
-                values[i] = values[i] > ValueType(0.0) ? ValueType(0.0) : values[i] * mWeight;
+                values[i] = values[i] > ValueType(0.0) ?
+                    ValueType(0.0) : ValueType(values[i] * mWeight);
                 if (values[i] > ValueType(0.0)) node.setValueOn(i);
             }
 
diff --git a/openvdb/openvdb/tools/MeshToVolume.h b/openvdb/openvdb/tools/MeshToVolume.h
index 931d01a0f9..1bfaa244ee 100644
--- a/openvdb/openvdb/tools/MeshToVolume.h
+++ b/openvdb/openvdb/tools/MeshToVolume.h
@@ -19,6 +19,7 @@
 #include <openvdb/Platform.h> // for OPENVDB_HAS_CXX11
 #include <openvdb/Types.h>
 #include <openvdb/math/FiniteDifference.h> // for GodunovsNormSqrd
+#include <openvdb/math/Math.h> // for isFinite(), isNan()
 #include <openvdb/math/Proximity.h> // for closestPointOnTriangleToPoint
 #include <openvdb/util/NullInterrupter.h>
 #include <openvdb/util/Util.h>
@@ -39,7 +40,6 @@
 #include <tbb/task_arena.h>
 
 #include <algorithm> // for std::sort()
-#include <cmath> // for std::isfinite(), std::isnan()
 #include <deque>
 #include <limits>
 #include <memory>
@@ -3348,7 +3348,7 @@ meshToVolume(
 
     // Note: inf interior width is all right, this value makes the converter fill
     // interior regions with distance values.
-    if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
+    if (!math::isFinite(exteriorWidth) || math::isNan(interiorWidth)) {
         std::stringstream msg;
         msg << "Illegal narrow band width: exterior = " << exteriorWidth
             << ", interior = " << interiorWidth;
@@ -3358,7 +3358,7 @@ meshToVolume(
 
     const ValueType voxelSize = ValueType(transform.voxelSize()[0]);
 
-    if (!std::isfinite(voxelSize) || math::isZero(voxelSize)) {
+    if (!math::isFinite(voxelSize) || math::isZero(voxelSize)) {
         std::stringstream msg;
         msg << "Illegal transform, voxel size = " << voxelSize;
         OPENVDB_LOG_DEBUG(msg.str());
diff --git a/openvdb/openvdb/tools/RayIntersector.h b/openvdb/openvdb/tools/RayIntersector.h
index d9b2d3a31e..7028b88286 100644
--- a/openvdb/openvdb/tools/RayIntersector.h
+++ b/openvdb/openvdb/tools/RayIntersector.h
@@ -91,7 +91,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
diff --git a/openvdb/openvdb/tools/RayTracer.h b/openvdb/openvdb/tools/RayTracer.h
index 659270176b..a3c93d2ab0 100644
--- a/openvdb/openvdb/tools/RayTracer.h
+++ b/openvdb/openvdb/tools/RayTracer.h
@@ -155,7 +155,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.
diff --git a/openvdb/openvdb/tools/SignedFloodFill.h b/openvdb/openvdb/tools/SignedFloodFill.h
index 228376ff4d..700ab2f0cd 100644
--- a/openvdb/openvdb/tools/SignedFloodFill.h
+++ b/openvdb/openvdb/tools/SignedFloodFill.h
@@ -85,7 +85,7 @@ class SignedFloodFillOp
     using ValueT = typename TreeOrLeafManagerT::ValueType;
     using RootT = typename TreeOrLeafManagerT::RootNodeType;
     using LeafT = typename TreeOrLeafManagerT::LeafNodeType;
-    static_assert(std::is_signed<ValueT>::value,
+    static_assert(openvdb::is_signed<ValueT>::value,
         "signed flood fill is supported only for signed value grids");
 
     SignedFloodFillOp(const TreeOrLeafManagerT& tree, Index minLevel = 0)
@@ -215,7 +215,7 @@ class SignedFloodFillOp
 
 template<typename TreeOrLeafManagerT>
 inline
-typename std::enable_if<std::is_signed<typename TreeOrLeafManagerT::ValueType>::value, void>::type
+typename std::enable_if_t<openvdb::is_signed<typename TreeOrLeafManagerT::ValueType>::value, void>
 doSignedFloodFill(TreeOrLeafManagerT& tree,
                   typename TreeOrLeafManagerT::ValueType outsideValue,
                   typename TreeOrLeafManagerT::ValueType insideValue,
@@ -231,7 +231,7 @@ doSignedFloodFill(TreeOrLeafManagerT& tree,
 // Dummy (no-op) implementation for unsigned types
 template <typename TreeOrLeafManagerT>
 inline
-typename std::enable_if<!std::is_signed<typename TreeOrLeafManagerT::ValueType>::value, void>::type
+typename std::enable_if_t<!openvdb::is_signed<typename TreeOrLeafManagerT::ValueType>::value, void>
 doSignedFloodFill(TreeOrLeafManagerT&,
                   const typename TreeOrLeafManagerT::ValueType&,
                   const typename TreeOrLeafManagerT::ValueType&,