Fix most floating point exceptions (#2266)

* better documentation, test against eigen slicing

* remove a bunch of slices

* clean up templates; static asserts on vector types

* final purges of slice

* fix templates in debug and ears test

* fix slice mask bug

* fix slice bug

* fix most fpe exceptions

* fix limit case

* more robust push, pop

include/igl/active_set.cpp

@@ -111,8 +111,6 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 
   // Keep track of previous Z for comparison
   DerivedZ old_Z;
-  old_Z = DerivedZ::Constant(
-      n,1,numeric_limits<typename DerivedZ::Scalar>::max());
 
   int iter = 0;
   while(true)
@@ -167,14 +165,17 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       cout<<"  new_as_lx: "<<new_as_lx<<endl;
       cout<<"  new_as_ux: "<<new_as_ux<<endl;
 #endif
-      const double diff = (Z-old_Z).squaredNorm();
+      if(iter > 0)
+      {
+        const double diff = (Z-old_Z).squaredNorm();
 #ifdef ACTIVE_SET_CPP_DEBUG
-      cout<<"diff: "<<diff<<endl;
+        cout<<"diff: "<<diff<<endl;
 #endif
-      if(diff < params.solution_diff_threshold)
-      {
-        ret = SOLVER_STATUS_CONVERGED;
-        break;
+        if(diff < params.solution_diff_threshold)
+        {
+          ret = SOLVER_STATUS_CONVERGED;
+          break;
+        }
       }
       old_Z = Z;
     }

include/igl/fast_winding_number.cpp

@@ -58,15 +58,22 @@ IGL_INLINE void igl::fast_winding_number(
       Eigen::Matrix<real_ec,1,3> zeroth_expansion;
       zeroth_expansion << 0,0,0;
       real_p areatotal = 0.0;
-      for(int j = 0; j < point_indices[index].size(); j++){
+      const int num_points = point_indices[index].size();
+      for(int j = 0; j < num_points; j++){
           int curr_point_index = point_indices[index][j];
         
           areatotal += A(curr_point_index);
           masscenter += A(curr_point_index)*P.row(curr_point_index);
           zeroth_expansion += A(curr_point_index)*N.row(curr_point_index);
       }
-    
-      masscenter = masscenter/areatotal;
+      // Avoid divide by zero
+      if(num_points > 0)
+      {
+        masscenter = masscenter/areatotal;
+      }else
+      {
+        masscenter.setConstant(std::numeric_limits<real_cm>::quiet_NaN());
+      }
       CM.row(index) = masscenter;
       EC.block(index,0,1,3) = zeroth_expansion;
     

include/igl/qslim_optimal_collapse_edge_callbacks.cpp

@@ -39,8 +39,17 @@ IGL_INLINE void igl::qslim_optimal_collapse_edge_callbacks(
     const auto & A = std::get<0>(quadric_p);
     const auto & b = std::get<1>(quadric_p);
     const auto & c = std::get<2>(quadric_p);
-    p = -b*A.inverse();
-    cost = p.dot(p*A) + 2*p.dot(b) + c;
+    if(b.array().isInf().any())
+    {
+      cost = std::numeric_limits<double>::infinity();
+      p.resizeLike(b);
+      p.setConstant(std::numeric_limits<double>::quiet_NaN());
+    }else
+    {
+      p = -b*A.inverse();
+      cost = p.dot(p*A) + 2*p.dot(b) + c;
+    }
+
     // Force infs and nans to infinity
     if(std::isinf(cost) || cost!=cost)
     {

include/igl/slim.cpp

@@ -410,7 +410,11 @@ IGL_INLINE void igl::slim_update_weights_and_closest_rotations_with_jacobians(co
         {
           double s1_g = 2 * (s1 - pow(s1, -3));
           double s2_g = 2 * (s2 - pow(s2, -3));
-          m_sing_new << sqrt(s1_g / (2 * (s1 - 1))), sqrt(s2_g / (2 * (s2 - 1)));
+          // Limit is 4 if s==1 according to Equation (32) in Rabinovich et al.
+          // [2017]
+          m_sing_new << 
+            (s1==1?4:sqrt(s1_g / (2 * (s1 - 1)))), 
+            (s2==1?4:sqrt(s2_g / (2 * (s2 - 1))));
           break;
         }
         case igl::MappingEnergyType::LOG_ARAP:

include/igl/solid_angle.cpp

@@ -56,6 +56,7 @@ IGL_INLINE typename DerivedA::Scalar igl::solid_angle(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>::Scalar igl::solid_angle<Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>, Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>, Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>, Eigen::Matrix<double, 1, 2, 1, 1, 2>>(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>> const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, 2, 0, -1, 2> const, 1, 2, false>> const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&);
 // generated by autoexplicit.sh
 template Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false>::Scalar igl::solid_angle<Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false>, Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false>, Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, 3, 0, -1, 3> const, 1, 3, false> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&);
 // generated by autoexplicit.sh

tests/include/igl/embree/EmbreeIntersector.cpp

@@ -3,6 +3,7 @@
 
 TEST_CASE("EmbreeIntersector: cube", "[igl/embree]")
 {
+  IGL_PUSH_FPE;
   //The allowed error for this test
   const double epsilon = 1e-6;
 
@@ -95,5 +96,6 @@ TEST_CASE("EmbreeIntersector: cube", "[igl/embree]")
     REQUIRE(hit.u == Approx(0.5).margin(epsilon));
     REQUIRE(hit.v == Approx(0.25).margin(epsilon));
   }
+  IGL_POP_FPE;
 }
 

tests/include/igl/ray_mesh_intersect.cpp

@@ -2,8 +2,10 @@
 #include <igl/ray_mesh_intersect.h>
 #include <igl/AABB.h>
 
+
 TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
 {
+  IGL_PUSH_FPE;
   Eigen::MatrixXd V(3,3);
   V.row(0) << 0.0, 0.0, 0.0;
   V.row(1) << 1.0, 0.0, 0.0;
@@ -23,11 +25,14 @@ TEST_CASE("ray_mesh_intersect: one_triangle", "[igl]")
   REQUIRE(igl::ray_mesh_intersect(source, direction, V, F, hits) == true);
   REQUIRE(hits.size() == 1);
   REQUIRE(hits.front().t == Approx(1.0));
+
+  IGL_POP_FPE;
 }
 
 
 TEST_CASE("ray_mesh_intersect: corner-case", "[igl]")
 {
+  IGL_PUSH_FPE;
   //       .      //
   //      /|\     //
   //     / | \    //
@@ -72,4 +77,5 @@ TEST_CASE("ray_mesh_intersect: corner-case", "[igl]")
     REQUIRE (is_hit == is_hit_bvh);
     REQUIRE (hits.size() == hits_bvh.size());
   }
+  IGL_POP_FPE;
 }

tests/main.cpp

@@ -6,19 +6,18 @@
 #include <catch2/catch.hpp>
 
 
-// TODO: Fix floating point exceptions raised in debug mode before re-enabling this.
-// #ifndef NDEBUG
-// #ifdef __linux__
-// #include <fenv.h>
-// #endif
-// #endif
+#ifndef NDEBUG
+#ifdef __linux__
+#include <fenv.h>
+#endif
+#endif
 
-// #ifndef NDEBUG
-// #ifdef __linux__
-// void beforeMain (void) __attribute__((constructor));
-// void beforeMain (void)
-// {
-//     feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
-// }
-// #endif
-// #endif
+#ifndef NDEBUG
+#ifdef __linux__
+void beforeMain (void) __attribute__((constructor));
+void beforeMain (void)
+{
+    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+}
+#endif
+#endif

tests/test_common.h

@@ -24,6 +24,20 @@
 #define IGL_DEBUG_OFF "[!hide]"
 #endif
 
+#if !defined(NDEBUG) && defined(__linux__)
+#include <fenv.h>
+#define IGL_PUSH_FPE \
+  fexcept_t current_exceptions; \
+  fegetexceptflag(&current_exceptions, FE_ALL_EXCEPT); \
+  fedisableexcept(FE_ALL_EXCEPT);
+#define IGL_POP_FPE \
+  /* Surprise! Use fesetexceptflag here not feenableexcept */ \
+  fesetexceptflag(&current_exceptions,FE_ALL_EXCEPT);
+#else
+// No-op
+#define IGL_PUSH_FPE 
+#define IGL_POP_FPE 
+#endif
 
 #include <igl/STR.h>
 template<>