Fix bug in random_point_on_mesh + crash with blue_noise on Windows.

include/igl/blue_noise.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2020 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "blue_noise.h"
 #include "doublearea.h"
@@ -126,6 +126,7 @@ namespace igl
     std::unordered_map<BlueNoiseKeyType,int> & S,
     std::vector<int> & active)
   {
+    assert(M.count(nk));
     auto & Mvec = M.find(nk)->second;
     auto miter = Mvec.begin();
     while(miter != Mvec.end())
@@ -149,7 +150,14 @@ namespace igl
         // back)
         //std::swap(*miter,Mvec.back());
         *miter = Mvec.back();
+        bool was_last = (std::next(miter) == Mvec.end());
         Mvec.pop_back();
+        if (was_last) {
+          // popping from the vector can invalidate the iterator, if it was
+          // pointing to the last element that was popped. Alternatively,
+          // one could use indices directly...
+          miter = Mvec.end();
+        }
       }
     }
     return false;
@@ -259,20 +267,22 @@ IGL_INLINE void igl::blue_noise(
   // function of r and s; and removing the `if S=-1 checks`)
   const Scalar s = r/sqrt(3.0);
 
-  const double area = 
+  const double area =
     [&](){Eigen::VectorXd A;igl::doublearea(V,F,A);return A.array().sum()/2;}();
   // Circle packing in the plane has igl::PI*sqrt(3)/6 efficiency
-  const double expected_number_of_points = 
+  const double expected_number_of_points =
     area * (igl::PI * sqrt(3.0) / 6.0) / (igl::PI * min_r * min_r / 4.0);
 
   // Make a uniform random sampling with 30*expected_number_of_points.
   const int nx = 30.0*expected_number_of_points;
   MatrixX3S X,XB;
   Eigen::VectorXi XFI;
+  std::cout << V.rows() << " " << V.cols() << std::endl;
+  std::cout << F.rows() << " " << F.cols() << std::endl;
   igl::random_points_on_mesh(nx,V,F,XB,XFI,X);
 
   // Rescale so that s = 1
-  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> Xs = 
+  Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> Xs =
     ((X.rowwise()-X.colwise().minCoeff())/s).template cast<int>();
   const int w = Xs.maxCoeff()+1;
   {
@@ -294,7 +304,7 @@ IGL_INLINE void igl::blue_noise(
   S.reserve(Xs.rows());
   for(int i = 0;i<Xs.rows();i++)
   {
-    BlueNoiseKeyType k = blue_noise_key(w,Xs(i,0),Xs(i,1),Xs(i,2)); 
+    BlueNoiseKeyType k = blue_noise_key(w,Xs(i,0),Xs(i,1),Xs(i,2));
     const auto Miter = M.find(k);
     if(Miter  == M.end())
     {
@@ -309,7 +319,7 @@ IGL_INLINE void igl::blue_noise(
 
   std::vector<int> active;
   // precompute r²
-  // Q: is this necessary? 
+  // Q: is this necessary?
   const double rr = r*r;
   std::vector<int> collected;
   collected.reserve(2.0*expected_number_of_points);

include/igl/random_points_on_mesh.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 Alec Jacobson <[email protected]>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "random_points_on_mesh.h"
 #include "doublearea.h"
@@ -42,6 +42,7 @@ IGL_INLINE void igl::random_points_on_mesh(
   assert(R.minCoeff() >= 0);
   assert(R.maxCoeff() <= 1);
   histc(R,C,FI);
+  FI = FI.array().min(F.rows() - 1); // fix the bin when R(i) == 1 exactly
   const VectorXs S = (VectorXs::Random(n,1).array() + 1.)/2.;
   const VectorXs T = (VectorXs::Random(n,1).array() + 1.)/2.;
   B.resize(n,3);
@@ -51,9 +52,9 @@ IGL_INLINE void igl::random_points_on_mesh(
 }
 
 template <
-  typename DerivedV, 
-  typename DerivedF, 
-  typename DerivedB, 
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedB,
   typename DerivedFI,
   typename DerivedX>
 IGL_INLINE void igl::random_points_on_mesh(
@@ -70,7 +71,10 @@ IGL_INLINE void igl::random_points_on_mesh(
   {
     for(int b = 0;b<B.cols();b++)
     {
-      X.row(x) += B(x,b)*V.row(F(FI(x),b));
+        auto fi = FI(x);
+        auto f = F(fi, b);
+        auto bval = B(x, b);
+        X.row(x) += bval*V.row(f);
     }
   }
 }