use nan replacement doublearea

Former-commit-id: e8b945a16b57eb55311293f47c98742187ce2eb2

include/igl/circumradius.cpp

@@ -20,7 +20,7 @@ IGL_INLINE void igl::circumradius(
   Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> l;
   igl::edge_lengths(V,F,l);
   DerivedR A;
-  igl::doublearea(l,A);
+  igl::doublearea(l,0.,A);
   // use formula: R=abc/(4*area) to compute the circum radius
   R = l.col(0).array() * l.col(1).array() * l.col(2).array() / (2.0*A.array());
 }

include/igl/cotmatrix_entries.cpp

@@ -44,7 +44,7 @@ IGL_INLINE void igl::cotmatrix_entries(
       
       // double area
       Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
-      doublearea(l,dblA);
+      doublearea(l,0.,dblA);
       // cotangents and diagonal entries for element matrices
       // correctly divided by 4 (alec 2010)
       C.resize(m,3);

include/igl/doublearea.cpp

@@ -67,7 +67,7 @@ IGL_INLINE void igl::doublearea(
     default:
     {
       edge_lengths(V,F,l);
-      return doublearea(l,dblA);
+      return doublearea(l,0.,dblA);
     }
   }
 }
@@ -131,6 +131,16 @@ template <typename Derivedl, typename DeriveddblA>
 IGL_INLINE void igl::doublearea(
   const Eigen::MatrixBase<Derivedl> & ul,
   Eigen::PlainObjectBase<DeriveddblA> & dblA)
+{
+  // Default is to leave NaNs and fire asserts in debug mode
+  return doublearea(ul,0.0/0.0,dblA);
+}
+
+template <typename Derivedl, typename DeriveddblA>
+IGL_INLINE void igl::doublearea(
+  const Eigen::MatrixBase<Derivedl> & ul,
+  const typename Derivedl::Scalar nan_replacement,
+  Eigen::PlainObjectBase<DeriveddblA> & dblA)
 {
   using namespace Eigen;
   using namespace std;
@@ -152,7 +162,7 @@ IGL_INLINE void igl::doublearea(
   dblA.resize(l.rows(),1);
   parallel_for(
     m,
-    [&l,&dblA](const int i)
+    [&l,&dblA,&nan_replacement](const int i)
     {
       // Kahan's Heron's formula
       typedef typename Derivedl::Scalar Scalar;
@@ -169,12 +179,16 @@ IGL_INLINE void igl::doublearea(
       // computation may be larger than (or rather smaller) than the height of
       // the triangle. In "Lecture Notes on Geometric Robustness" Shewchuck 09,
       // Section 3.1 http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf,
-      // he recommends computing area of a triangle with vertices nD (if you
-      // have access to them) by computing the determinant of the n+1 × n+1
-      // matrix containg the vertices in n+1D homogeneous coordinates as rows
-      // and dividing by d! (factorial).
-      assert( ((l(i,2) - (l(i,0)-l(i,1)))>=-igl::EPS<Scalar>()) 
+      // he recommends computing the triangle areas for 2D and 3D using 2D
+      // signed areas computed with determinants.
+      assert( 
+        (nan_replacement == nan_replacement || 
+          (l(i,2) - (l(i,0)-l(i,1)))>=0)
           && "Side lengths do not obey the triangle inequality.");
+      if(dblA(i) != dblA(i))
+      {
+        dblA(i) = nan_replacement;
+      }
       assert(dblA(i) == dblA(i) && "DOUBLEAREA() PRODUCED NaN");
     },
     1000l);

include/igl/doublearea.h

@@ -55,17 +55,32 @@ namespace igl
     const Eigen::MatrixBase<DerivedA> & A,
     const Eigen::MatrixBase<DerivedB> & B,
     const Eigen::MatrixBase<DerivedC> & C);
-  // Same as above but use instrinsic edge lengths rather than (V,F) mesh
+  // Same as above but use instrinsic edge lengths rather than (V,F) mesh. This
+  //
   // Inputs:
   //   l  #F by dim list of edge lengths using
   //     for triangles, columns correspond to edges 23,31,12
+  //   nan_replacement  what value should be used for triangles whose given
+  //     edge lengths do not obey the triangle inequality. These may be very
+  //     wrong (e.g., [100 1 1]) or may be nearly degenerate triangles whose
+  //     floating point side length computation leads to breach of the triangle
+  //     inequality. One may wish to set this parameter to 0 if side lengths l
+  //     are _known_ to come from a valid embedding (e.g., some mesh (V,F)). In
+  //     that case, the only circumstance the triangle inequality is broken is
+  //     when the triangle is nearly degenerate and floating point error
+  //     dominates: hence replacing with zero is reasonable.
   // Outputs:
   //   dblA  #F list of triangle double areas
   template <typename Derivedl, typename DeriveddblA>
+  IGL_INLINE void doublearea(
+    const Eigen::MatrixBase<Derivedl> & l,
+    const typename Derivedl::Scalar nan_replacement,
+    Eigen::PlainObjectBase<DeriveddblA> & dblA);
+  // default behavior is to assert on NaNs and leave them in place
+  template <typename Derivedl, typename DeriveddblA>
   IGL_INLINE void doublearea(
     const Eigen::MatrixBase<Derivedl> & l,
     Eigen::PlainObjectBase<DeriveddblA> & dblA);
-
   // DOUBLEAREA_QUAD computes twice the area for each input quadrilateral
   //
   // Inputs:

include/igl/face_areas.cpp

@@ -25,6 +25,15 @@ template <typename DerivedL, typename DerivedA>
 IGL_INLINE void igl::face_areas(
   const Eigen::MatrixBase<DerivedL>& L,
   Eigen::PlainObjectBase<DerivedA>& A)
+{
+  return face_areas(L,0./0.,A);
+}
+
+template <typename DerivedL, typename DerivedA>
+IGL_INLINE void igl::face_areas(
+  const Eigen::MatrixBase<DerivedL>& L,
+  const typename DerivedL::Scalar doublearea_nan_replacement,
+  Eigen::PlainObjectBase<DerivedA>& A)
 {
   using namespace Eigen;
   assert(L.cols() == 6);

include/igl/face_areas.h

@@ -41,6 +41,12 @@ namespace igl
   IGL_INLINE void face_areas(
     const Eigen::MatrixBase<DerivedL>& L,
     Eigen::PlainObjectBase<DerivedA>& A);
+  // doublearea_nan_replacement  See doublearea.h
+  template <typename DerivedL, typename DerivedA>
+  IGL_INLINE void face_areas(
+    const Eigen::MatrixBase<DerivedL>& L,
+    const typename DerivedL::Scalar doublearea_nan_replacement,
+    Eigen::PlainObjectBase<DerivedA>& A);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/inradius.cpp

@@ -28,6 +28,6 @@ IGL_INLINE void igl::inradius(
   // r/(4*area) = 1/(2*(a+b+c))
   // r = (2*area)/(a+b+c)
   DerivedR A;
-  igl::doublearea(l,A);
+  igl::doublearea(l,0.,A);
   r = A.array() /l.array().rowwise().sum();
 }