From ffdcbb51149fad45f19bd1017ad9f64c78ff5ffb Mon Sep 17 00:00:00 2001
From: Kacper Pluta <kacper.pluta@bsdlabs.com.br>
Date: Tue, 23 Oct 2018 12:23:22 +0300
Subject: [PATCH] avoid unsafe calls to acos by fixing potential floating-point
 errors and avoiding division by zero.

---
 include/igl/copyleft/comiso/nrosy.cpp | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/include/igl/copyleft/comiso/nrosy.cpp b/include/igl/copyleft/comiso/nrosy.cpp
index 0f3ee86bc..c334e015b 100644
--- a/include/igl/copyleft/comiso/nrosy.cpp
+++ b/include/igl/copyleft/comiso/nrosy.cpp
@@ -16,6 +16,8 @@
 #include <iostream>
 #include <fstream>
 
+#include <stdexcept>
+
 #include <Eigen/Geometry>
 #include <Eigen/Sparse>
 #include <queue>
@@ -803,8 +805,16 @@ Eigen::VectorXd igl::copyleft::comiso::NRosyField::angleDefect()
     {
       Eigen::VectorXd a = V.row(F(i,(j+1)%3)) - V.row(F(i,j));
       Eigen::VectorXd b = V.row(F(i,(j+2)%3)) - V.row(F(i,j));
-      double t = a.transpose()*b;
-      t /= (a.norm() * b.norm());
+      double t = a.transpose() * b;
+      double norm_prod = a.norm() * b.norm();
+      if ( norm_prod > 0. )
+        t /= norm_prod;
+      else
+        throw std::runtime_error ( "Error in 'igl::copyleft::comiso::NRosyField::angleDefect': division by zero! " );
+      if ( t > 1. )
+        t = 1.;
+      else if ( t < -1. )
+        t = -1.;
       A(F(i,j)) += acos(t);
     }
   }