orient_outward_ao: separate distribution data for each thread, remove two-sided faces using unique_simplices

Former-commit-id: d8cad1a86dd7368a3168fd34b302c7f9d39e4735

examples/patches/example.cpp

@@ -29,6 +29,7 @@
 #include <igl/boost/bfs_orient.h>
 #include <igl/orient_outward.h>
 #include <igl/embree/orient_outward_ao.h>
+#include <igl/unique_simplices.h>
 
 #include <Eigen/Core>
 #include <Eigen/Geometry>
@@ -697,6 +698,9 @@ int main(int argc, char * argv[])
     }
     triangulate(vF,F);
   }
+  MatrixXi F_unique;
+  unique_simplices(F, F_unique);
+  F = F_unique;
 
   init_patches();
 

include/igl/embree/orient_outward_ao.cpp

@@ -47,41 +47,44 @@ IGL_INLINE void igl::orient_outward_ao(
   PlainObjectBase<DerivedV> N;
   per_face_normals(V,F,N);
   
+  // random number generator/distribution for each thread
+  int max_threads = omp_get_max_threads();
+  
+  // prng
+  vector<mt19937> engine(max_threads);
+  for (int i = 0; i < max_threads; ++i)
+      engine[i].seed(time(0) * (i + 1));
+  
   // discrete distribution for random selection of faces with probability proportional to their areas
   Matrix<typename DerivedV::Scalar,Dynamic,1> A;
   doublearea(V,F,A);
   double minarea = A.minCoeff();
   Matrix<int, Dynamic, 1> A_int = (A * 100.0 / minarea).template cast<int>();       // only integer is allowed for weight
   auto ddist_func = [&] (double i) { return A_int(static_cast<int>(i)); };
-  discrete_distribution<int> ddist(m, 0, m, ddist_func);      // simple ctor of (Iter, Iter) not provided by the stupid VC11 impl...
+  vector<discrete_distribution<int>> ddist(max_threads, discrete_distribution<int>(m, 0, m, ddist_func));      // simple ctor of (Iter, Iter) not provided by the stupid VC11 impl...
   
-  uniform_real_distribution<double> rdist;
+  // uniform real between in [0, 1]
+  vector<uniform_real_distribution<double>> rdist(max_threads);
   
   // occlusion count per component: +1 when front ray is occluded, -1 when back ray is occluded
   Matrix<int, Dynamic, 1> C_occlude_count;
   C_occlude_count.setZero(m, 1);
   
-  // prng for each omp thread
-  int max_threads = omp_get_max_threads();
-  vector<mt19937> engine(max_threads);
-  for (int i = 0; i < max_threads; ++i)
-      engine[i].seed(time(0) * (i + 1));
-  
 #pragma omp parallel for
   for (int i = 0; i < num_samples; ++i)
   {
     int thread_num = omp_get_thread_num();
-    int f     = ddist(engine[thread_num]);   // select face with probability proportional to face area
-    double t0 = rdist(engine[thread_num]);
-    double t1 = rdist(engine[thread_num]);
-    double t2 = rdist(engine[thread_num]);
+    int f     = ddist[thread_num](engine[thread_num]);   // select face with probability proportional to face area
+    double t0 = rdist[thread_num](engine[thread_num]);
+    double t1 = rdist[thread_num](engine[thread_num]);
+    double t2 = rdist[thread_num](engine[thread_num]);
     double t_sum = t0 + t1 + t2;
     t0 /= t_sum;
     t1 /= t_sum;
     t2 /= t_sum;
     RowVector3d p = t0 * V.row(F(f,0)) + t1 * V.row(F(f,1)) + t1 * V.row(F(f,2));
     RowVector3d n = N.row(f);
-    bool is_backside = rdist(engine[thread_num]) < 0.5;
+    bool is_backside = rdist[thread_num](engine[thread_num]) < 0.5;
     if (is_backside)
     {
         n *= -1;