doublearea now more robust

Former-commit-id: 2950ff7b4f27f0237dc0bcd29656e87c05446f69

Makefile

@@ -20,6 +20,8 @@ debug: OPTFLAGS= -g -Wall
 debug: DEBUG=debug
 CFLAGS += $(OPTFLAGS)
 #CFLAGS += -DIGL_NO_OPENGL -DIGL_NO_ANTTWEAKBAR
+# We use well-supported features of c++11
+CFLAGS += -std=c++11
 
 EXTRA_DIRS=
 ifeq ($(IGL_WITH_TETGEN),1)

examples/quicklook-mesh/src/render_to_buffer.cpp

@@ -34,7 +34,7 @@ static int width,height;
 static Eigen::MatrixXd V,N;
 static Eigen::MatrixXi F;
 static Eigen::Vector3d Vmean, Vmax,Vmin;
-static bool invert = false;
+//static bool invert = false;
 static float background_color[4] = {0,0,0,1};
 
 // Small viewports struct for keeping track of size and camera info
@@ -268,27 +268,36 @@ void display()
     push_scene(viewports[vp]);
     push_object(viewports[vp]);
 
+    // Draw the mesh, inverted if need be.
     // Set material properties
     glDisable(GL_COLOR_MATERIAL);
-    glMaterialfv(GL_FRONT, GL_AMBIENT,  GOLD_AMBIENT);
-    glMaterialfv(GL_FRONT, GL_DIFFUSE,  GOLD_DIFFUSE  );
-    glMaterialfv(GL_FRONT, GL_SPECULAR, GOLD_SPECULAR);
-    glMaterialf (GL_FRONT, GL_SHININESS, 128);
-    glMaterialfv(GL_BACK, GL_AMBIENT,  SILVER_AMBIENT);
-    glMaterialfv(GL_BACK, GL_DIFFUSE,  FAST_GREEN_DIFFUSE  );
-    glMaterialfv(GL_BACK, GL_SPECULAR, SILVER_SPECULAR);
-    glMaterialf (GL_BACK, GL_SHININESS, 128);
-
-    // Draw the mesh, inverted if need be.
-    if(invert)
-    {
-      glFrontFace(GL_CW);
-    }
+    //if(invert)
+    //{
+    //  glFrontFace(GL_CW);
+    //  glMaterialfv(GL_FRONT, GL_AMBIENT,  CYAN_AMBIENT);
+    //  glMaterialfv(GL_FRONT, GL_DIFFUSE,  CYAN_DIFFUSE  );
+    //  glMaterialfv(GL_FRONT, GL_SPECULAR, CYAN_SPECULAR);
+    //  glMaterialf (GL_FRONT, GL_SHININESS, 128);
+    //  glMaterialfv(GL_BACK, GL_AMBIENT,  SILVER_AMBIENT);
+    //  glMaterialfv(GL_BACK, GL_DIFFUSE,  FAST_RED_DIFFUSE);
+    //  glMaterialfv(GL_BACK, GL_SPECULAR, SILVER_SPECULAR);
+    //  glMaterialf (GL_BACK, GL_SHININESS, 128);
+    //}else
+    //{
+      glMaterialfv(GL_FRONT, GL_AMBIENT,  GOLD_AMBIENT);
+      glMaterialfv(GL_FRONT, GL_DIFFUSE,  GOLD_DIFFUSE  );
+      glMaterialfv(GL_FRONT, GL_SPECULAR, GOLD_SPECULAR);
+      glMaterialf (GL_FRONT, GL_SHININESS, 128);
+      glMaterialfv(GL_BACK, GL_AMBIENT,  SILVER_AMBIENT);
+      glMaterialfv(GL_BACK, GL_DIFFUSE,  FAST_GREEN_DIFFUSE  );
+      glMaterialfv(GL_BACK, GL_SPECULAR, SILVER_SPECULAR);
+      glMaterialf (GL_BACK, GL_SHININESS, 128);
+    //}
     draw_mesh(V,F,N);
-    if(invert)
-    {
-      glFrontFace(GL_CCW);
-    }
+    //if(invert)
+    //{
+    //  glFrontFace(GL_CCW);
+    //}
     pop_object();
 
     // Draw a nice floor unless we're looking from beneath it
@@ -411,20 +420,20 @@ bool render_to_buffer(
   Vmin = V.colwise().minCoeff();
   Vmean = 0.5*(Vmax + Vmin);
 
-  // Figure out if normals should be flipped (hopefully this is never a
-  // bottleneck)
-  MatrixXd BC;
-  VectorXd dblA;
-  barycenter(V,F,BC);
-  BC.col(0).array() -= Vmean(0,0);
-  BC.col(1).array() -= Vmean(1,0);
-  BC.col(2).array() -= Vmean(2,0);
-  doublearea(V,F,dblA);
-  VectorXd BCDN = (BC.array() * N.array()).rowwise().sum();
-  const double tot_dp = dblA.transpose() * BCDN;
-  invert = tot_dp < 0;
-  cout<<"Normals: "<<(get_seconds()-ts)<<"s"<<endl;
-  ts = get_seconds();
+  //// Figure out if normals should be flipped (hopefully this is never a
+  //// bottleneck)
+  //MatrixXd BC;
+  //VectorXd dblA;
+  //barycenter(V,F,BC);
+  //BC.col(0).array() -= Vmean(0,0);
+  //BC.col(1).array() -= Vmean(1,0);
+  //BC.col(2).array() -= Vmean(2,0);
+  //doublearea(V,F,dblA);
+  //VectorXd BCDN = (BC.array() * N.array()).rowwise().sum();
+  //const double tot_dp = dblA.transpose() * BCDN;
+  //invert = tot_dp < 0;
+  //cout<<"Normals: "<<(get_seconds()-ts)<<"s"<<endl;
+  //ts = get_seconds();
 
   // Initialize MESA
   OSMesaContext ctx;

include/igl/doublearea.cpp

@@ -1,6 +1,8 @@
 #include "doublearea.h"
 #include "edge_lengths.h"
+#include "sort.h"
 #include <cassert>
+#include <iostream>
 
 template <typename DerivedV, typename DerivedF, typename DeriveddblA>
 IGL_INLINE void igl::doublearea( 
@@ -8,33 +10,84 @@ IGL_INLINE void igl::doublearea(
   const Eigen::PlainObjectBase<DerivedF> & F, 
   Eigen::PlainObjectBase<DeriveddblA> & dblA)
 {
+  const int dim = V.cols();
   // Only support triangles
   assert(F.cols() == 3);
+  const int m = F.rows();
   // Compute edge lengths
   Eigen::PlainObjectBase<DerivedV> l;
-  edge_lengths(V,F,l);
-  return doublearea(l,dblA);
+  // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
+  // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
+  switch(dim)
+  {
+    case 3:
+    {
+      dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
+      for(int d = 0;d<3;d++)
+      {
+        Eigen::Matrix<typename DerivedV::Scalar, DerivedV::RowsAtCompileTime, 2> Vd(V.rows(),2);
+        Vd << V.col(d), V.col((d+1)%3);
+        Eigen::PlainObjectBase<DeriveddblA> dblAd;
+        doublearea(Vd,F,dblAd);
+        dblA.array() += dblAd.array().square();
+      }
+      dblA = dblA.array().sqrt().eval();
+      break;
+    }
+    case 2:
+    {
+      dblA.resize(m,1);
+      for(int f = 0;f<m;f++)
+      {
+        auto r = V.row(F(f,0))-V.row(F(f,2));
+        auto s = V.row(F(f,1))-V.row(F(f,2));
+        dblA(f) = r(0)*s(1) - r(1)*s(0);
+      }
+      break;
+    }
+    default:
+    {
+      edge_lengths(V,F,l);
+      return doublearea(l,dblA);
+    }
+  }
 }
 
 template <typename Derivedl, typename DeriveddblA>
 IGL_INLINE void igl::doublearea( 
-  const Eigen::PlainObjectBase<Derivedl> & l, 
+  const Eigen::PlainObjectBase<Derivedl> & ul,
   Eigen::PlainObjectBase<DeriveddblA> & dblA)
 {
   using namespace Eigen;
+  using namespace std;
   // Only support triangles
-  assert(l.cols() == 3);
+  assert(ul.cols() == 3);
   // Number of triangles
-  const int m = l.rows();
+  const int m = ul.rows();
+  Eigen::PlainObjectBase<Derivedl> l;
+  MatrixXi _;
+  sort(ul,2,false,l,_);
   // semiperimeters
   Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
   assert(s.rows() == m);
   // resize output
   dblA.resize(l.rows(),1);
-  // Heron's formula for area
   for(int i = 0;i<m;i++)
   {
-    dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
+    //// Heron's formula for area
+    //const typename Derivedl::Scalar arg = 
+    //  s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2));
+    //assert(arg>=0);
+    //dblA(i) = 2.0*sqrt(arg);
+    // Kahan's Heron's formula
+    const typename Derivedl::Scalar arg = 
+      (l(i,0)+(l(i,1)+l(i,2)))*
+      (l(i,2)-(l(i,0)-l(i,1)))*
+      (l(i,2)+(l(i,0)-l(i,1)))*
+      (l(i,0)+(l(i,1)-l(i,2)));
+    dblA(i) = 2.0*0.25*sqrt(arg);
+    assert( l(i,2) - (l(i,0)-l(i,1)) && "FAILED KAHAN'S ASSERTION");
+    assert(dblA(i) == dblA(i) && "DOUBLEAREA() PRODUCED NaN");
   }
 }
 

include/igl/embree/orient_outward_ao.cpp

@@ -71,7 +71,7 @@ IGL_INLINE void igl::orient_outward_ao(
   }
   
   // generate all the rays
-  cout << "generating rays... ";
+  //cout << "generating rays... ";
   uniform_real_distribution<float> rdist;
   mt19937 prng;
   prng.seed(time(nullptr));
@@ -139,7 +139,7 @@ IGL_INLINE void igl::orient_outward_ao(
   vector<pair<float, float>> C_vote_distance(num_cc, make_pair(0, 0));     // sum of distance between ray origin and intersection
   vector<pair<int  , int  >> C_vote_infinity(num_cc, make_pair(0, 0));     // number of rays reaching infinity
   
-  cout << "shooting rays... ";
+  //cout << "shooting rays... ";
 #pragma omp parallel for
   for (int i = 0; i < (int)ray_face.size(); ++i)
   {
@@ -191,7 +191,7 @@ IGL_INLINE void igl::orient_outward_ao(
       FF.row(f) = FF.row(f).reverse().eval();
     }
   }
-  cout << "done!\n";
+  //cout << "done!\n";
 }
 
 // Call with default parameters

include/igl/readOBJ.cpp

@@ -167,9 +167,13 @@ IGL_INLINE bool igl::readOBJ(
           fclose(obj_file);
           return false;
         }
-      }else if(strlen(type) >= 1 && type[0] == '#')
+      }else if(strlen(type) >= 1 && (type[0] == '#' || 
+            type[0] == 'g'  ||
+            type[0] == 's'  ||
+            strcmp("usemtl",type)==0 ||
+            strcmp("mtllib",type)==0))
       {
-        //ignore comments
+        //ignore comments or other shit
       }else
       {
         //ignore any other lines