refactored per_face_normals with option for degenerate fill in. Fixed compile issue with orient_outward_ao

Former-commit-id: b9ac0746c623a90c36abd839dc84ab668bba8e9f

examples/patches/example.cpp

@@ -56,6 +56,7 @@
 #include <stack>
 #include <iostream>
 
+int cc_selected = -1;
 
 Eigen::MatrixXd V;
 Eigen::VectorXd Vmid,Vmin,Vmax;
@@ -330,6 +331,38 @@ void display()
     {
       draw_mesh(V,F,s.N,s.C);
     }
+  
+    // visualize selected patch
+    glLineWidth(10);
+    glBegin(GL_TRIANGLES);
+    glColor3d(0, 0, 0);
+    // loop over faces
+    for(int i = 0; i<F.rows();i++)
+    {
+      if (CC(i) != cc_selected) continue;
+      // loop over corners of triangle
+      for(int j = 0;j<3;j++)
+      {
+        glVertex3d(V(F(i,j),0),V(F(i,j),1),V(F(i,j),2));
+      }
+    }
+    glEnd();
+    glLineWidth(1);
+    glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+    glBegin(GL_TRIANGLES);
+    glColor3d(1, 0, 0);
+    // loop over faces
+    for(int i = 0; i<F.rows();i++)
+    {
+      if (CC(i) != cc_selected) continue;
+      // loop over corners of triangle
+      glNormal3d(s.N(i,0),s.N(i,1),s.N(i,2));
+      for(int j = 0;j<3;j++)
+      {
+        glVertex3d(V(F(i,j),0),V(F(i,j),1),V(F(i,j),2));
+      }
+    }
+    glEnd();
   }
   glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
   if(fill_visible)
@@ -675,6 +708,11 @@ void key(unsigned char key, int mouse_x, int mouse_y)
     case 'j':
         mouse_wheel(0,-1,mouse_x,mouse_y);
         break;
+    case 'n':
+      cc_selected = (cc_selected + 1) % (CC.maxCoeff() + 2);
+      cout << "selected cc: " << cc_selected << endl;
+      glutPostRedisplay();
+      break;
     default:
       if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
       {

examples/patches/example.vcxproj

@@ -71,7 +71,7 @@
       <WarningLevel>Level3</WarningLevel>
       <PrecompiledHeader>
       </PrecompiledHeader>
-      <Optimization>MaxSpeed</Optimization>
+      <Optimization>Disabled</Optimization>
       <FunctionLevelLinking>true</FunctionLevelLinking>
       <IntrinsicFunctions>true</IntrinsicFunctions>
       <PreprocessorDefinitions>IGL_HEADER_ONLY;WIN32;NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
@@ -83,10 +83,11 @@
       <SubSystem>Console</SubSystem>
       <GenerateDebugInformation>true</GenerateDebugInformation>
       <EnableCOMDATFolding>true</EnableCOMDATFolding>
-      <OptimizeReferences>true</OptimizeReferences>
+      <OptimizeReferences>false</OptimizeReferences>
     </Link>
   </ItemDefinitionGroup>
   <ItemGroup>
+    <ClCompile Include="..\..\..\kt84\dbg.cpp" />
     <ClCompile Include="..\..\external\embree\embree\builders\heuristic_binning.cpp" />
     <ClCompile Include="..\..\external\embree\embree\builders\heuristic_spatial.cpp" />
     <ClCompile Include="..\..\external\embree\embree\builders\primrefgen.cpp" />

examples/patches/example.vcxproj.filters

@@ -143,6 +143,7 @@
     <ClCompile Include="..\..\external\embree\embree\embree.cpp">
       <Filter>embree</Filter>
     </ClCompile>
+    <ClCompile Include="..\..\..\kt84\dbg.cpp" />
   </ItemGroup>
   <ItemGroup>
     <Filter Include="igl">

include/igl/embree/orient_outward_ao.cpp

@@ -5,7 +5,6 @@
 #include "EmbreeIntersector.h"
 #include <iostream>
 #include <random>
-#include <limits>
 
 template <
   typename DerivedV, 
@@ -28,12 +27,8 @@ IGL_INLINE void igl::orient_outward_ao(
   assert(F.cols() == 3);
   assert(V.cols() == 3);
   
-  // pass both sides of faces to Embree
-  MatrixXi F2;
-  F2.resize(F.rows()*2,F.cols());
-  F2 << F, F.rowwise().reverse().eval();
   EmbreeIntersector ei;
-  ei.init(V.template cast<float>(),F2.template cast<int>());
+  ei.init(V.template cast<float>(),F);
   
   // number of faces
   const int m = F.rows();
@@ -46,7 +41,7 @@ IGL_INLINE void igl::orient_outward_ao(
   }
   
   // face normal
-  PlainObjectBase<DerivedV> N;
+  MatrixXd N;
   per_face_normals(V,F,N);
   
   // face area
@@ -71,12 +66,12 @@ IGL_INLINE void igl::orient_outward_ao(
   
   // generate all the rays
   cout << "generating rays... ";
-  uniform_real_distribution<double> rdist;
+  uniform_real_distribution<float> rdist;
   mt19937 prng;
-  prng.seed(time(0));
+  prng.seed(0);
   vector<int     > ray_face;
-  vector<Vector3d> ray_ori;
-  vector<Vector3d> ray_dir;
+  vector<Vector3f> ray_ori;
+  vector<Vector3f> ray_dir;
   ray_face.reserve(total_num_rays);
   ray_ori .reserve(total_num_rays);
   ray_dir .reserve(total_num_rays);
@@ -98,21 +93,32 @@ IGL_INLINE void igl::orient_outward_ao(
     for (int i = 0; i < num_rays_per_component[c]; ++i)
     {
       int f     = CF[ddist(prng)];    // select face with probability proportional to face area
-      double t0 = rdist(prng);        // random barycentric coordinate
-      double t1 = rdist(prng);
-      double t2 = rdist(prng);
-      double t_sum = t0 + t1 + t2;
+      float t0 = rdist(prng);        // random barycentric coordinate
+      float t1 = rdist(prng);
+      float t2 = rdist(prng);
+      float t_sum = t0 + t1 + t2;
       t0 /= t_sum;
       t1 /= t_sum;
       t2 /= t_sum;
-      Vector3d p = t0 * V.row(F(f,0))       // be careful with the index!!!
-                 + t1 * V.row(F(f,1))
-                 + t2 * V.row(F(f,2));
-      Vector3d n = N.row(f);
-      Vector3d d = random_dir();
-      if (n.dot(d) < 0)
-      {
-        d *= -1;
+      Vector3f p = t0 * V.row(F(f,0)).template cast<float>().eval()       // be careful with the index!!!
+                 + t1 * V.row(F(f,1)).template cast<float>().eval()
+                 + t2 * V.row(F(f,2)).template cast<float>().eval();
+      Vector3f n = N.row(f).cast<float>();
+      assert(n != Vector3f::Zero());
+      // random direction in hemisphere around n (avoid too grazing angle)
+      Vector3f d;
+      while (true) {
+        d = random_dir().cast<float>();
+        float ndotd = n.dot(d);
+        if (fabsf(ndotd) < 0.1f)
+        {
+          continue;
+        }
+        if (ndotd < 0)
+        {
+          d *= -1.0f;
+        }
+        break;
       }
       ray_face.push_back(f);
       ray_ori .push_back(p);
@@ -120,40 +126,53 @@ IGL_INLINE void igl::orient_outward_ao(
     }
   }
   
-  // per component accumulation of occlusion distance
-  double dist_large = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm() * 1000;
-  vector<double> C_occlude_dist_front(num_cc, 0);
-  vector<double> C_occlude_dist_back (num_cc, 0);
-
-  auto get_dist = [&] (Hit hit, const Vector3d& origin) {
-    Vector3d p0 = V.row(F2(hit.id, 0));
-    Vector3d p1 = V.row(F2(hit.id, 1));
-    Vector3d p2 = V.row(F2(hit.id, 2));
-    Vector3d p = (1 - hit.u - hit.v) * p0 + hit.u * p1 + hit.v * p2;
-    return (p - origin).norm();
-  };
+  // per component voting: first=front, second=back
+  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... ";
 #pragma omp parallel for
   for (int i = 0; i < (int)ray_face.size(); ++i)
   {
     int      f = ray_face[i];
-    Vector3d o = ray_ori [i];
-    Vector3d d = ray_dir [i];
+    Vector3f o = ray_ori [i];
+    Vector3f d = ray_dir [i];
     int c = C(f);
-    Hit hit_front;
-    Hit hit_back;
-    double dist_front = ei.intersectRay(o.template cast<float>(),  d.template cast<float>(), hit_front) ? get_dist(hit_front, o) : dist_large;
-    double dist_back  = ei.intersectRay(o.template cast<float>(), -d.template cast<float>(), hit_back ) ? get_dist(hit_back , o) : dist_large;
+    
+    // shoot ray toward front & back
+    vector<Hit> hits_front;
+    vector<Hit> hits_back;
+    int num_rays_front;
+    int num_rays_back;
+    ei.intersectRay(o,  d, hits_front, num_rays_front);
+    ei.intersectRay(o, -d, hits_back , num_rays_back );
+    if (!hits_front.empty() && hits_front[0].id == f) hits_front.erase(hits_front.begin());
+    if (!hits_back .empty() && hits_back [0].id == f) hits_back .erase(hits_back .begin());
+    
+    if (hits_front.empty())
+    {
+#pragma omp atomic
+      C_vote_infinity[c].first++;
+    } else {
 #pragma omp atomic
-    C_occlude_dist_front[c] += dist_front;
+      C_vote_distance[c].first += hits_front[0].t;
+    }
+    
+    if (hits_back.empty())
+    {
 #pragma omp atomic
-    C_occlude_dist_back [c] += dist_back;
+      C_vote_infinity[c].second++;
+    } else {
+#pragma omp atomic
+      C_vote_distance[c].second += hits_back[0].t;
+    }
   }
   
   for(int c = 0;c<num_cc;c++)
   {
-    I(c) = C_occlude_dist_front[c] < C_occlude_dist_back[c];
+    I(c) = C_vote_infinity[c].first == C_vote_infinity[c].second &&
+           C_vote_distance[c].first <  C_vote_distance[c].second ||
+           C_vote_infinity[c].first <  C_vote_infinity[c].second;
   }
   // flip according to I
   for(int f = 0;f<m;f++)

include/igl/orient_outward.cpp

@@ -36,7 +36,8 @@ IGL_INLINE void igl::orient_outward(
   PlainObjectBase<DerivedV> N,BC,BCmean;
   Matrix<typename DerivedV::Scalar,Dynamic,1> A;
   VectorXd totA(num_cc), dot(num_cc);
-  per_face_normals(V,F,N);
+  Matrix<typename DerivedV::Scalar,3,1> Z(1,1,1);
+  per_face_normals(V,F,Z.normalized(),N);
   barycenter(V,F,BC);
   doublearea(V,F,A);
   BCmean.setConstant(num_cc,3,0);

include/igl/per_face_normals.cpp

@@ -2,11 +2,12 @@
 #include <Eigen/Geometry>
 
 #define SQRT_ONE_OVER_THREE 0.57735026918962573
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedV, typename DerivedF, typename DerivedZ, typename DerivedN>
 IGL_INLINE void igl::per_face_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
-  Eigen::PlainObjectBase<DerivedV> & N)
+  const Eigen::PlainObjectBase<DerivedZ> & Z,
+  Eigen::PlainObjectBase<DerivedN> & N)
 {
   N.resize(F.rows(),3);
   // loop over faces
@@ -17,20 +18,30 @@ IGL_INLINE void igl::per_face_normals(
     Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = V.row(F(i,1)) - V.row(F(i,0));
     Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v2 = V.row(F(i,2)) - V.row(F(i,0));
     N.row(i) = v1.cross(v2);//.normalized();
-    if(N.row(i).sum() == 0)
+    typename DerivedV::Scalar r = N.row(i).norm();
+    if(r == 0)
     {
-      N(i,0) = SQRT_ONE_OVER_THREE;
-      N(i,1) = SQRT_ONE_OVER_THREE;
-      N(i,2) = SQRT_ONE_OVER_THREE;
+      N.row(i) = Z;
     }else
     {
-      N.row(i) = N.row(i).normalized().eval();
+      N.row(i) /= r;
     }
   }
 }
+
+template <typename DerivedV, typename DerivedF, typename DerivedN>
+IGL_INLINE void igl::per_face_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedN> & N)
+{
+  using namespace Eigen;
+  Matrix<typename DerivedN::Scalar,3,1> Z(0,0,0);
+  return per_face_normals(V,F,Z,N);
+}
+
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
-// generated by autoexplicit.sh
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::per_face_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/per_face_normals.h

@@ -7,14 +7,27 @@ namespace igl
   // Compute face normals via vertex position list, face list
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3 eigen Matrix of face (triangle) indices
+  //   Z  3 vector normal given to faces with degenerate normal.
   // Output:
   //   N  #F by 3 eigen Matrix of mesh face (triangle) 3D normals
-  template <typename DerivedV, typename DerivedF>
+  //
+  // Example:
+  //   // Give degenerate faces (1/3,1/3,1/3)^0.5
+  //   per_face_normals(V,F,Vector3d(1,1,1).normalized(),N);
+  template <typename DerivedV, typename DerivedF, typename DerivedZ, typename DerivedN>
   IGL_INLINE void per_face_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
-    Eigen::PlainObjectBase<DerivedV> & N);
+    const Eigen::PlainObjectBase<DerivedZ> & Z,
+    Eigen::PlainObjectBase<DerivedN> & N);
+  // Wrapper with Z = (0,0,0). Note that this means that row norms will be zero
+  // (i.e. not 1) for degenerate normals.
+  template <typename DerivedV, typename DerivedF, typename DerivedN>
+  IGL_INLINE void per_face_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedN> & N);
 }
 
 #ifdef IGL_HEADER_ONLY