performance bug in doublearea

Former-commit-id: a2ddfab299a096e98d6f37b27dad4ac69e907e84

include/igl/WindingNumberAABB.h

@@ -142,7 +142,6 @@ inline void igl::WindingNumberAABB<Point>::grow()
 
   // Blerg, why is selecting rows so difficult
 
-  vector<int> id( this->getF().rows());
   double split_value;
   // Split in longest direction
   switch(split_method)
@@ -160,6 +159,7 @@ inline void igl::WindingNumberAABB<Point>::grow()
   //cout<<"c: "<<0.5*(max_corner[max_d] + min_corner[max_d])<<" "<<
   //  "m: "<<split_value<<endl;;
 
+  vector<int> id( this->getF().rows());
   for(int i = 0;i<this->getF().rows();i++)
   {
     if(BC(i,max_d) <= split_value)
@@ -216,7 +216,7 @@ inline bool igl::WindingNumberAABB<Point>::inside(const Point & p) const
   {
     //// Perfect matching is **not** robust
     //if( p(i) < min_corner(i) || p(i) >= max_corner(i))
-    // **MUST** be conservative!!
+    // **MUST** be conservative
     if( p(i) < min_corner(i) || p(i) > max_corner(i))
     {
       return false;

include/igl/doublearea.cpp

@@ -25,18 +25,30 @@ IGL_INLINE void igl::doublearea(
   Eigen::PlainObjectBase<DerivedV> l;
   // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
   // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
+
+  // Projected area helper
+  const auto & proj_doublearea = 
+    [&V,&F](const int x, const int y, const int f)->double
+  {
+    auto rx = V(F(f,0),x)-V(F(f,2),x);
+    auto sx = V(F(f,1),x)-V(F(f,2),x);
+    auto ry = V(F(f,0),y)-V(F(f,2),y);
+    auto sy = V(F(f,1),y)-V(F(f,2),y);
+    return rx*sy - ry*sx;
+  };
+
   switch(dim)
   {
     case 3:
     {
       dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
-      for(int d = 0;d<3;d++)
+      for(int f = 0;f<F.rows();f++)
       {
-        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();
+        for(int d = 0;d<3;d++)
+        {
+          double dblAd = proj_doublearea(d,(d+1)%3,f);
+          dblA[f] += dblAd*dblAd;
+        }
       }
       dblA = dblA.array().sqrt().eval();
       break;
@@ -46,9 +58,7 @@ IGL_INLINE void igl::doublearea(
       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);
+        dblA(f) = proj_doublearea(0,1,f);
       }
       break;
     }

include/igl/doublearea.h

@@ -26,6 +26,8 @@ namespace igl
   // Outputs:
   //   dblA  #F list of triangle double areas (SIGNED only for 2D input)
   //
+  // Known bug: For dim==3 complexity is O(#V + #F)!! Not just O(#F). This is a big deal
+  // if you have 1million unreferenced vertices and 1 face
   template <typename DerivedV, typename DerivedF, typename DeriveddblA>
   IGL_INLINE void doublearea( 
     const Eigen::PlainObjectBase<DerivedV> & V,