diff --git a/openvdb/openvdb/math/Stencils.h b/openvdb/openvdb/math/Stencils.h
index 0c4ed2c0c8..5d145815d9 100644
--- a/openvdb/openvdb/math/Stencils.h
+++ b/openvdb/openvdb/math/Stencils.h
@@ -1627,7 +1627,9 @@ 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 base = mean*mean - alphaG*mInvDx2/math::Pow4(normGrad);
+          const Real clampVal = 0;
+          const Real tmp = std::sqrt(std::max(base, clampVal));
           pair.first  = ValueType(mean - tmp);
           pair.second = ValueType(mean + tmp);
         }
diff --git a/openvdb/openvdb/points/IndexFilter.h b/openvdb/openvdb/points/IndexFilter.h
index 9f0e57d708..0c44126133 100644
--- a/openvdb/openvdb/points/IndexFilter.h
+++ b/openvdb/openvdb/points/IndexFilter.h
@@ -240,7 +240,8 @@ class RandomLeafFilter
             currentPoints += iter->pointCount();
         }
 
-        const float factor = targetPoints > currentPoints ? 1.0f : float(targetPoints) / float(currentPoints);
+        const float denom = currentPoints > 0 ? float(currentPoints) : 1.0f; // Do not divide by zero.
+        const float factor = targetPoints > currentPoints ? 1.0f : float(targetPoints) / denom;
 
         std::mt19937 generator(seed);
         std::uniform_int_distribution<unsigned int> dist(0, std::numeric_limits<unsigned int>::max() - 1);
diff --git a/openvdb/openvdb/tools/impl/ConvexVoxelizer.h b/openvdb/openvdb/tools/impl/ConvexVoxelizer.h
index 5ce7226a7d..e00a09746d 100644
--- a/openvdb/openvdb/tools/impl/ConvexVoxelizer.h
+++ b/openvdb/openvdb/tools/impl/ConvexVoxelizer.h
@@ -389,7 +389,10 @@ class ConvexVoxelizer
     circleBottom(const ValueT& x0, const ValueT& y0,
                  const ValueT& r, const ValueT& x)
     {
-        return y0 - math::Sqrt(math::Pow2(r) - math::Pow2(x-x0));
+        const ValueT base = math::Pow2(r) - math::Pow2(x-x0);
+        return base <= 0
+            ? y0
+            : y0 - math::Sqrt(base);
     }
 
     /// @brief Computes the top y-coordinate of a circle at a given x position.
@@ -402,7 +405,10 @@ class ConvexVoxelizer
     circleTop(const ValueT& x0, const ValueT& y0,
               const ValueT& r, const ValueT& x)
     {
-        return y0 + math::Sqrt(math::Pow2(r) - math::Pow2(x-x0));
+        const ValueT base = math::Pow2(r) - math::Pow2(x-x0);
+        return base <= 0
+            ? y0
+            : y0 + math::Sqrt(base);
     }
 
     /// @brief Computes the bottom z-coordinate of a sphere at a given (x, y) position.
@@ -417,7 +423,10 @@ class ConvexVoxelizer
     sphereBottom(const ValueT& x0, const ValueT& y0, const ValueT& z0,
                  const ValueT& r, const ValueT& x, const ValueT& y)
     {
-        return z0 - math::Sqrt(math::Pow2(r) - math::Pow2(x-x0) - math::Pow2(y-y0));
+        const ValueT base = math::Pow2(r) - math::Pow2(x-x0) - math::Pow2(y-y0);
+        return base <= 0 // Do not take square root of negative numbers.
+            ? z0
+            : z0 - math::Sqrt(base);
     }
 
     /// @brief Computes the top z-coordinate of a sphere at a given (x, y) position.
@@ -432,7 +441,10 @@ class ConvexVoxelizer
     sphereTop(const ValueT& x0, const ValueT& y0, const ValueT& z0,
               const ValueT& r, const ValueT& x, const ValueT& y)
     {
-        return z0 + math::Sqrt(math::Pow2(r) - math::Pow2(x-x0) - math::Pow2(y-y0));
+        const ValueT base = math::Pow2(r) - math::Pow2(x-x0) - math::Pow2(y-y0);
+        return base <= 0 // Do not take square root of negative numbers.
+            ? z0
+            : z0 + math::Sqrt(base);
     }
 
     // ------------ nested classes ------------
diff --git a/openvdb/openvdb/unittest/CMakeLists.txt b/openvdb/openvdb/unittest/CMakeLists.txt
index 184e96a5d2..eed19b738c 100644
--- a/openvdb/openvdb/unittest/CMakeLists.txt
+++ b/openvdb/openvdb/unittest/CMakeLists.txt
@@ -18,6 +18,8 @@ set(OPENVDB_TESTS "" CACHE STRING [=[
   \"Activate;NodeManager\" would build just the TestActivate.cc and TestNodeManager.cc
   unit tests.]=])
 
+option(OPENVDB_TESTS_FPE "Enable FP exceptions in unit tests (linux only)" FALSE)
+
 ##########################################################################
 
 message(STATUS "----------------------------------------------------")
@@ -196,6 +198,10 @@ else()
   )
 endif()
 
+if (OPENVDB_TESTS_FPE)
+  add_compile_definitions(OPENVDB_TESTS_FPE)
+endif()
+
 add_executable(vdb_test ${UNITTEST_SOURCE_FILES})
 
 # Blosc and ZLib are hidden dependencies for the core library
diff --git a/openvdb/openvdb/unittest/main.cc b/openvdb/openvdb/unittest/main.cc
index 5066512862..6faf197c9e 100644
--- a/openvdb/openvdb/unittest/main.cc
+++ b/openvdb/openvdb/unittest/main.cc
@@ -12,11 +12,25 @@
 #include <string>
 #include <vector>
 
+#if defined(__linux__) && defined(OPENVDB_TESTS_FPE)
+#include <fenv.h>
+#endif
+
 #include <gtest/gtest.h>
 
 int
 main(int argc, char *argv[])
 {
     ::testing::InitGoogleTest(&argc, argv);
+
+#if defined(__linux__) && defined(OPENVDB_TESTS_FPE)
+    int excepts = FE_DIVBYZERO | FE_INVALID;
+    if (getenv("OPENVDB_TEST_OVERFLOW")) {
+        // when set, test for FP overflow as well.
+        excepts = excepts | FE_OVERFLOW;
+    }
+    feenableexcept(excepts);
+#endif
+
     return RUN_ALL_TESTS();
 }