Fix blue_noise and random_points_on_mesh reproducibility (#2235)

* Fix `blue_noise` and `random_points_on_mesh` reproducibility

- added new overload for `blue_noise` and `random_points_on_mesh` for that accept an UnformRandomBitGenerator as input. The default signature (ie without URBG) still works, and uses by default std::minstd_rand initialized with random seed generated by std::rand(). This is following the existing behaviour of `randperm`.
- added reproducibility test case for `blue_noise` and `random_points_on_mesh`

* `blue_noise` and `random_points_on_mesh` overloads without an URBG parameter now use a fixed default seed

ie `std::minstd_rand()`

default seed

* Fix static compilation of tutorials

* revert to jdumas std::rand seed for default

* missing templates on linux

* more missing linux templates

* use std names rather than expansions in explicit templates

---------

Co-authored-by: stourneux <[email protected]>
Co-authored-by: seb-tourneux <[email protected]>

include/igl/blue_noise.cpp

@@ -166,12 +166,14 @@ namespace igl
 
   template <
     typename DerivedX,
-    typename DerivedXs>
+    typename DerivedXs,
+    typename URBG>
   inline bool step(
     const Eigen::MatrixBase<DerivedX> & X,
     const Eigen::MatrixBase<DerivedXs> & Xs,
     const double & rr,
     const int & w,
+    URBG && urbg,
     std::unordered_map<BlueNoiseKeyType,std::vector<int> > & M,
     std::unordered_map<BlueNoiseKeyType,int> & S,
     std::vector<int> & active,
@@ -181,7 +183,8 @@ namespace igl
     //considered.clear();
     if(active.size() == 0) return false;
     // random entry
-    const int e = rand() % active.size();
+    std::uniform_int_distribution<> dis(0, active.size()-1);
+    const int e = dis(urbg);
     const int i = active[e];
     //printf("%d\n",i);
     const int xi = Xs(i,0);
@@ -210,8 +213,7 @@ namespace igl
     }
         //printf("  --------\n");
     // randomize order: this might be a little paranoid...
-    std::mt19937 twister;
-    std::shuffle(std::begin(N), std::end(N), twister);
+    std::shuffle(std::begin(N), std::end(N), urbg);
     bool found = false;
     for(const BlueNoiseKeyType & nk : N)
     {
@@ -243,14 +245,16 @@ template <
   typename DerivedF,
   typename DerivedB,
   typename DerivedFI,
-  typename DerivedP>
+  typename DerivedP,
+  typename URBG>
 IGL_INLINE void igl::blue_noise(
     const Eigen::MatrixBase<DerivedV> & V,
     const Eigen::MatrixBase<DerivedF> & F,
     const typename DerivedV::Scalar r,
     Eigen::PlainObjectBase<DerivedB> & B,
     Eigen::PlainObjectBase<DerivedFI> & FI,
-    Eigen::PlainObjectBase<DerivedP> & P)
+    Eigen::PlainObjectBase<DerivedP> & P,
+    URBG && urbg)
 {
   typedef typename DerivedV::Scalar Scalar;
   typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
@@ -279,7 +283,7 @@ IGL_INLINE void igl::blue_noise(
   const int nx = 30.0*expected_number_of_points;
   MatrixX3S X,XB;
   Eigen::VectorXi XFI;
-  igl::random_points_on_mesh(nx,V,F,XB,XFI,X);
+  igl::random_points_on_mesh(nx,V,F,XB,XFI,X,urbg);
 
   // Rescale so that s = 1
   Eigen::Matrix<int,Eigen::Dynamic,3,Eigen::RowMajor> Xs =
@@ -350,7 +354,7 @@ IGL_INLINE void igl::blue_noise(
   {
     while(active.size()>0)
     {
-      step(X,Xs,rr,w,M,S,active,collected);
+      step(X,Xs,rr,w,urbg,M,S,active,collected);
     }
   }
   {
@@ -369,6 +373,7 @@ IGL_INLINE void igl::blue_noise(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::blue_noise<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<float, -1, 3, 1, -1, 3>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::blue_noise<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// Explicit template instantiation
+template void igl::blue_noise<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::mt19937_64 >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::mt19937_64&&);
+template void igl::blue_noise<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::mt19937 >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::mt19937&&);
 #endif

include/igl/blue_noise.h

@@ -8,6 +8,7 @@
 #ifndef IGL_BLUE_NOISE_H
 #define IGL_BLUE_NOISE_H
 #include "igl_inline.h"
+#include "generate_default_urbg.h"
 #include <Eigen/Core>
 namespace igl
 {
@@ -24,20 +25,25 @@ namespace igl
   ///               ith sampled point in face FI(i)
   /// @param[out] FI  #P list of indices into F 
   /// @param[out] P  #P by dim list of sample positions.
+  /// @param[in,out] urbg An instance of UnformRandomBitGenerator (e.g.,
+  ///  `std::minstd_rand(0)`)
   /// \see random_points_on_mesh
   template <
     typename DerivedV,
     typename DerivedF,
     typename DerivedB,
     typename DerivedFI,
-    typename DerivedP>
+    typename DerivedP,
+    typename URBG = DEFAULT_URBG
+      >
   IGL_INLINE void blue_noise(
       const Eigen::MatrixBase<DerivedV> & V,
       const Eigen::MatrixBase<DerivedF> & F,
       const typename DerivedV::Scalar r,
       Eigen::PlainObjectBase<DerivedB> & B,
       Eigen::PlainObjectBase<DerivedFI> & FI,
-      Eigen::PlainObjectBase<DerivedP> & P);
+      Eigen::PlainObjectBase<DerivedP> & P,
+      URBG && urbg = igl::generate_default_urbg());
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/generate_default_urbg.h

@@ -0,0 +1,13 @@
+#ifndef IGL_GENERATE_DEFAULT_URBG_H
+#define IGL_GENERATE_DEFAULT_URBG_H
+#include <random>
+
+namespace igl {
+  using DEFAULT_URBG = std::mt19937;
+  inline DEFAULT_URBG generate_default_urbg() 
+  {
+    return DEFAULT_URBG(std::rand());
+  }
+}
+#endif
+

include/igl/random_points_on_mesh.cpp

@@ -6,126 +6,55 @@
 // 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 "random_points_on_mesh_intrinsic.h"
 #include "doublearea.h"
-#include "cumsum.h"
-#include "histc.h"
-#include <iostream>
-#include <cassert>
 
-template <typename DerivedV, typename DerivedF, typename DerivedB, typename DerivedFI>
-IGL_INLINE void igl::random_points_on_mesh(
-  const int n,
-  const Eigen::MatrixBase<DerivedV > & V,
-  const Eigen::MatrixBase<DerivedF > & F,
-  Eigen::PlainObjectBase<DerivedB > & B,
-  Eigen::PlainObjectBase<DerivedFI > & FI)
-{
-  using namespace Eigen;
-  using namespace std;
-  typedef typename DerivedV::Scalar Scalar;
-  typedef Matrix<Scalar,Dynamic,1> VectorXs;
-  VectorXs A;
-  doublearea(V,F,A);
-  // Should be traingle mesh. Although Turk's method 1 generalizes...
-  assert(F.cols() == 3);
-  VectorXs C;
-  VectorXs A0(A.size()+1);
-  A0(0) = 0;
-  A0.bottomRightCorner(A.size(),1) = A;
-  // Even faster would be to use the "Alias Table Method"
-  cumsum(A0,1,C);
-  const Scalar Cmax = C(C.size()-1);
-  assert(Cmax > 0 && "Total surface area should be positive");
-  // Why is this more accurate than `C /= C(C.size()-1)` ?
-  for(int i = 0;i<C.size();i++) { C(i) = C(i)/Cmax; }
-  const VectorXs R = (VectorXs::Random(n,1).array() + 1.)/2.;
-  assert(R.minCoeff() >= 0);
-  assert(R.maxCoeff() <= 1);
-  histc(R,C,FI);
-  // fix the bin when R(i) == 1 exactly
-  // Gross cast to deal with Windows 
-  FI = FI.array().min(static_cast<typename DerivedFI::Scalar>(F.rows() - 1));
-  const VectorXs S = (VectorXs::Random(n,1).array() + 1.)/2.;
-  const VectorXs T = (VectorXs::Random(n,1).array() + 1.)/2.;
-  B.resize(n,3);
-  B.col(0) = 1.-T.array().sqrt();
-  B.col(1) = (1.-S.array()) * T.array().sqrt();
-  B.col(2) = S.array() * T.array().sqrt();
-}
 
 template <
   typename DerivedV,
   typename DerivedF,
   typename DerivedB,
   typename DerivedFI,
-  typename DerivedX>
+  typename DerivedX,
+  typename URBG>
 IGL_INLINE void igl::random_points_on_mesh(
   const int n,
   const Eigen::MatrixBase<DerivedV > & V,
   const Eigen::MatrixBase<DerivedF > & F,
   Eigen::PlainObjectBase<DerivedB > & B,
   Eigen::PlainObjectBase<DerivedFI > & FI,
-  Eigen::PlainObjectBase<DerivedX> & X)
+  Eigen::PlainObjectBase<DerivedX> & X,
+  URBG && urbg)
 {
-  random_points_on_mesh(n,V,F,B,FI);
+  using namespace Eigen;
+  using namespace std;
+  typedef typename DerivedV::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1> VectorXs;
+  VectorXs dblA;
+  doublearea(V,F,dblA);
+  random_points_on_mesh_intrinsic(n,dblA,B,FI,urbg);
   X = DerivedX::Zero(B.rows(),V.cols());
   for(int x = 0;x<B.rows();x++)
   {
     for(int b = 0;b<B.cols();b++)
     {
-        auto fi = FI(x);
-        auto f = F(fi, b);
-        auto bval = B(x, b);
-        X.row(x) += bval*V.row(f);
+      auto fi = FI(x);
+      auto f = F(fi, b);
+      auto bval = B(x, b);
+      X.row(x) += bval*V.row(f);
     }
   }
 }
 
-template <typename DerivedV, typename DerivedF, typename ScalarB, typename DerivedFI>
-IGL_INLINE void igl::random_points_on_mesh(
-  const int n,
-  const Eigen::MatrixBase<DerivedV > & V,
-  const Eigen::MatrixBase<DerivedF > & F,
-  Eigen::SparseMatrix<ScalarB > & B,
-  Eigen::PlainObjectBase<DerivedFI > & FI)
-{
-  using namespace Eigen;
-  using namespace std;
-  Matrix<ScalarB,Dynamic,3> BC;
-  random_points_on_mesh(n,V,F,BC,FI);
-  vector<Triplet<ScalarB> > BIJV;
-  BIJV.reserve(n*3);
-  for(int s = 0;s<n;s++)
-  {
-    for(int c = 0;c<3;c++)
-    {
-      assert(FI(s) < F.rows());
-      assert(FI(s) >= 0);
-      const int v = F(FI(s),c);
-      BIJV.push_back(Triplet<ScalarB>(s,v,BC(s,c)));
-    }
-  }
-  B.resize(n,V.rows());
-  B.reserve(n*3);
-  B.setFromTriplets(BIJV.begin(),BIJV.end());
-}
-
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(int, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, float, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(int, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<float, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-
-template 
-void igl::random_points_on_mesh<Eigen::MatrixXd,Eigen::MatrixXi,Eigen::MatrixXd,Eigen::VectorXi>(
-  const int n,
-  const Eigen::MatrixBase<Eigen::MatrixXd> & V,
-  const Eigen::MatrixBase<Eigen::MatrixXi> & F,
-  Eigen::PlainObjectBase< Eigen::MatrixXd> & B,
-  Eigen::PlainObjectBase< Eigen::VectorXi > & FI);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, std::mt19937_64&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, std::mt19937_64&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::minstd_rand0&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::minstd_rand0&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::minstd_rand&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::minstd_rand&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::mt19937 >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::mt19937&&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::mt19937&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::mt19937&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::mt19937_64&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::mt19937_64&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>,   std::mt19937 >(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::mt19937&&);
+template void igl::random_points_on_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>,     std::mt19937&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, std::mt19937&);
 #endif

include/igl/random_points_on_mesh.h

@@ -9,6 +9,7 @@
 #define IGL_RANDOM_POINTS_ON_MESH_H
 
 #include "igl_inline.h"
+#include "generate_default_urbg.h"
 #include <Eigen/Core>
 #include <Eigen/Sparse>
 
@@ -21,40 +22,26 @@ namespace igl
   /// @param[in] F  #F by 3 list of mesh triangle indices
   /// @param[out] B  n by 3 list of barycentric coordinates, ith row are coordinates of
   ///     ith sampled point in face FI(i)
+  /// @param[in] urbg An instance of UnformRandomBitGenerator (e.g.,
+  ///  `std::minstd_rand(0)`)
   /// @param[out] FI  n list of indices into F 
-  ///
-  template <typename DerivedV, typename DerivedF, typename DerivedB, typename DerivedFI>
-  IGL_INLINE void random_points_on_mesh(
-    const int n,
-    const Eigen::MatrixBase<DerivedV > & V,
-    const Eigen::MatrixBase<DerivedF > & F,
-    Eigen::PlainObjectBase<DerivedB > & B,
-    Eigen::PlainObjectBase<DerivedFI > & FI);
-  /// \overload
+  /// @param[in,out] urbg An instance of UnformRandomBitGenerator.
   /// @param[out] X  n by dim list of sample positions.
   template <
     typename DerivedV, 
     typename DerivedF, 
     typename DerivedB, 
     typename DerivedFI,
-    typename DerivedX>
+    typename DerivedX,
+    typename URBG = DEFAULT_URBG>
   IGL_INLINE void random_points_on_mesh(
     const int n,
     const Eigen::MatrixBase<DerivedV > & V,
     const Eigen::MatrixBase<DerivedF > & F,
     Eigen::PlainObjectBase<DerivedB > & B,
     Eigen::PlainObjectBase<DerivedFI > & FI,
-    Eigen::PlainObjectBase<DerivedX> & X);
-  /// \overload
-  /// @param[out] B n by #V sparse matrix so that  B*V produces a list of sample
-  ///   points
-  template <typename DerivedV, typename DerivedF, typename ScalarB, typename DerivedFI>
-  IGL_INLINE void random_points_on_mesh(
-    const int n,
-    const Eigen::MatrixBase<DerivedV > & V,
-    const Eigen::MatrixBase<DerivedF > & F,
-    Eigen::SparseMatrix<ScalarB > & B,
-    Eigen::PlainObjectBase<DerivedFI > & FI);
+    Eigen::PlainObjectBase<DerivedX> & X,
+    URBG && urbg = igl::generate_default_urbg());
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/random_points_on_mesh_intrinsic.cpp

@@ -0,0 +1,97 @@
+#include "random_points_on_mesh_intrinsic.h"
+#include "cumsum.h"
+#include "histc.h"
+#include <cassert>
+#include <random>
+
+template <
+  typename DeriveddblA,
+  typename DerivedB, 
+  typename DerivedFI,
+  typename URBG>
+IGL_INLINE void igl::random_points_on_mesh_intrinsic(
+  const int n,
+  const Eigen::MatrixBase<DeriveddblA > & dblA,
+  Eigen::PlainObjectBase<DerivedB > & B,
+  Eigen::PlainObjectBase<DerivedFI > & FI,
+  URBG && urbg)
+{
+  using namespace Eigen;
+  using namespace std;
+  typedef typename DeriveddblA::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1> VectorXs;
+  // Should be traingle mesh. Although Turk's method 1 generalizes...
+  assert(F.cols() == 3);
+  VectorXs C;
+  VectorXs A0(dblA.size()+1);
+  A0(0) = 0;
+  A0.bottomRightCorner(dblA.size(),1) = dblA;
+  // Even faster would be to use the "Alias Table Method"
+  cumsum(A0,1,C);
+  const Scalar Cmax = C(C.size()-1);
+  assert(Cmax > 0 && "Total surface area should be positive");
+  // Why is this more accurate than `C /= C(C.size()-1)` ?
+  for(int i = 0;i<C.size();i++) { C(i) = C(i)/Cmax; }
+  std::uniform_real_distribution<Scalar> dis(-1.0, 1.0);
+  const VectorXs R = (VectorXs::NullaryExpr(n,1,[&](){return dis(urbg);}).array() + 1.)/2.;
+  assert(R.minCoeff() >= 0);
+  assert(R.maxCoeff() <= 1);
+  histc(R,C,FI);
+  // fix the bin when R(i) == 1 exactly
+  // Gross cast to deal with Windows 
+  FI = FI.array().min(static_cast<typename DerivedFI::Scalar>(dblA.rows() - 1));
+  const VectorXs S = (VectorXs::NullaryExpr(n,1,[&](){return dis(urbg);}).array() + 1.)/2.;
+  const VectorXs T = (VectorXs::NullaryExpr(n,1,[&](){return dis(urbg);}).array() + 1.)/2.;
+  B.resize(n,3);
+  B.col(0) = 1.-T.array().sqrt();
+  B.col(1) = (1.-S.array()) * T.array().sqrt();
+  B.col(2) = S.array() * T.array().sqrt();
+}
+
+template <
+  typename DeriveddblA,
+  typename DerivedF,
+  typename ScalarB, 
+  typename DerivedFI,
+  typename URBG>
+IGL_INLINE void igl::random_points_on_mesh_intrinsic(
+  const int n,
+  const Eigen::MatrixBase<DeriveddblA > & dblA,
+  const int num_vertices,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::SparseMatrix<ScalarB > & B,
+  Eigen::PlainObjectBase<DerivedFI > & FI,
+  URBG && urbg)
+{
+  using namespace Eigen;
+  using namespace std;
+  Matrix<ScalarB,Dynamic,3> BC;
+  random_points_on_mesh_intrinsic(n,dblA,BC,FI,urbg);
+  vector<Triplet<ScalarB> > BIJV;
+  BIJV.reserve(n*3);
+  for(int s = 0;s<n;s++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      assert(FI(s) < dblA.rows());
+      assert(FI(s) >= 0);
+      const int v = F(FI(s),c);
+      BIJV.push_back(Triplet<ScalarB>(s,v,BC(s,c)));
+    }
+  }
+  B.resize(n,num_vertices);
+  B.reserve(n*3);
+  B.setFromTriplets(BIJV.begin(),BIJV.end());
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand&);
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand0&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand0&);
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937&);
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937_64&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937_64&);
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937&);
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937&);
+template void igl::random_points_on_mesh_intrinsic<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937_64&>(int, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937_64&);
+#endif

include/igl/random_points_on_mesh_intrinsic.h

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_RANDOM_POINTS_ON_MESH_INTRINSIC_H
+#define IGL_RANDOM_POINTS_ON_MESH_INTRINSIC_H
+
+#include "igl_inline.h"
+#include "generate_default_urbg.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  /// Randomly sample a mesh (V,F) n times.
+  ///
+  /// @param[in] n  number of samples
+  /// @param[in] dblA  #F list of double areas of triangles
+  /// @param[out] B  n by 3 list of barycentric coordinates, ith row are coordinates of
+  ///     ith sampled point in face FI(i)
+  /// @param[in] urbg An instance of UnformRandomBitGenerator (e.g.,
+  ///  `std::minstd_rand(0)`)
+  /// @param[out] FI  n list of indices into F 
+  /// @param[in,out] urbg An instance of UnformRandomBitGenerator.
+  template <
+    typename DeriveddblA,
+    typename DerivedB, 
+    typename DerivedFI,
+    typename URBG = DEFAULT_URBG>
+  IGL_INLINE void random_points_on_mesh_intrinsic(
+    const int n,
+    const Eigen::MatrixBase<DeriveddblA > & dblA,
+    Eigen::PlainObjectBase<DerivedB > & B,
+    Eigen::PlainObjectBase<DerivedFI > & FI,
+    URBG && urbg = igl::generate_default_urbg());
+  /// \overload
+  ///
+  /// @param[in] num_vertices  number of vertices in mesh
+  /// @param[in] F  #F by 3 list of mesh triangle indices
+  /// @param[out] B n by num_vertices sparse matrix so that  B*V produces a list
+  ///   of sample points if dbl = doublearea(V,F)
+  ///
+  template <
+    typename DeriveddblA,
+    typename DerivedF,
+    typename ScalarB, 
+    typename DerivedFI,
+    typename URBG>
+  IGL_INLINE void random_points_on_mesh_intrinsic(
+    const int n,
+    const Eigen::MatrixBase<DeriveddblA > & dblA,
+    const int num_vertices,
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::SparseMatrix<ScalarB > & B,
+    Eigen::PlainObjectBase<DerivedFI > & FI,
+    URBG && urbg);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "random_points_on_mesh_intrinsic.cpp"
+#endif
+
+#endif
+

include/igl/randperm.cpp

@@ -22,52 +22,16 @@ IGL_INLINE void igl::randperm(
   std::shuffle(I.data(),I.data()+n, urbg);
 }
 
-template <typename DerivedI>
-IGL_INLINE void igl::randperm(
-  const int n,
-  Eigen::PlainObjectBase<DerivedI> & I)
-{
-  return igl::randperm(n, I, std::minstd_rand(std::rand()));
-}
-
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand0>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand0 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand0 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand0 &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::minstd_rand0>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::minstd_rand0 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::minstd_rand0 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::minstd_rand0 &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::minstd_rand>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::minstd_rand &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::minstd_rand &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::minstd_rand &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937 &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937 &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937_64>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937_64 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937_64 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937_64 &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937_64>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937_64 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937_64 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937_64 &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux24_base>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux24_base &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux24_base &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux24_base &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux24_base>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux24_base &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux24_base &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux24_base &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux48_base>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux48_base &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux48_base &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux48_base &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux48_base>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux48_base &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux48_base &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux48_base &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux24>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux24 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux24 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux24 &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux24>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux24 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux24 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux24 &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux48>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux48 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::ranlux48 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::ranlux48 &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux48>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux48 &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::ranlux48 &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::ranlux48 &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::knuth_b>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::knuth_b &&);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::knuth_b &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::knuth_b &);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::knuth_b>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::knuth_b &&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::knuth_b &>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::knuth_b &);
-template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::minstd_rand0&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::minstd_rand0&);
+template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::minstd_rand&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::minstd_rand&);
+template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937_64&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937_64&);
+template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937&);
+template void igl::randperm<Eigen::Matrix<int, -1, -1, 0, -1, -1>, std::mt19937 >(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::mt19937&&);
+template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand0&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand0&);
+template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::minstd_rand&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::minstd_rand&);
+template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937_64&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937_64&);
+template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937&>(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937&);
+template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::mt19937 >(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::mt19937&&);
 #endif

include/igl/randperm.h

@@ -8,6 +8,7 @@
 #ifndef IGL_RANDPERM_H
 #define IGL_RANDPERM_H
 #include "igl_inline.h"
+#include "generate_default_urbg.h"
 #include <Eigen/Core>
 #include <random>
 namespace igl
@@ -19,17 +20,12 @@ namespace igl
   ///
   /// @param[in] n  number of elements
   /// @param[out] I  n list of rand permutation of 0:n-1
-  /// @param[in] urbg An instance of UnformRandomBitGenerator.
-  template <typename DerivedI, typename URBG>
+  /// @param[in,out] urbg An instance of UnformRandomBitGenerator.
+  template <typename DerivedI, typename URBG = DEFAULT_URBG>
   IGL_INLINE void randperm(
     const int n,
     Eigen::PlainObjectBase<DerivedI> & I,
-    URBG && urbg);
-  /// \overload
-  template <typename DerivedI>
-  IGL_INLINE void randperm(
-    const int n,
-    Eigen::PlainObjectBase<DerivedI> & I);
+    URBG && urbg = igl::generate_default_urbg());
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "randperm.cpp"

tests/include/igl/blue_noise.cpp

@@ -13,28 +13,95 @@
 #include <igl/octree.h>
 #include <igl/slice.h>
 
-TEST_CASE("blue_noise: decimated-knight", "[igl]")
-{
-  Eigen::MatrixXd V;
-  Eigen::MatrixXi F;
-  igl::readOBJ(test_common::data_path("decimated-knight.obj"),V,F);
-  const double r = 0.01;
-  Eigen::MatrixXd B,P;
-  {
-    Eigen::VectorXi I;
-    igl::blue_noise(V,F,r,B,I,P);
-  }
-  // There should be ~4000 samples on this model
-  REQUIRE(P.rows() > 3000);
-  std::vector<std::vector<int> > point_indices;
-  Eigen::MatrixXi CH;
-  Eigen::MatrixXd CN;
-  Eigen::VectorXd W;
-  igl::octree(P,point_indices,CH,CN,W);
-  Eigen::MatrixXi I;
-  igl::knn(P,2,point_indices,CH,CN,W,I);
-  Eigen::MatrixXd P2;
-  igl::slice(P,I.col(1).eval(),1,P2);
-  Eigen::VectorXd D = (P-P2).rowwise().norm();
-  REQUIRE(D.minCoeff() > r);
-}
+//TEST_CASE("blue_noise: decimated-knight", "[igl]")
+//{
+//  Eigen::MatrixXd V;
+//  Eigen::MatrixXi F;
+//  igl::readOBJ(test_common::data_path("decimated-knight.obj"),V,F);
+//  const double r = 0.01;
+//  Eigen::MatrixXd B,P;
+//  {
+//    Eigen::VectorXi I;
+//    igl::blue_noise(V,F,r,B,I,P);
+//  }
+//  // There should be ~4000 samples on this model
+//  REQUIRE(P.rows() > 3000);
+//  std::vector<std::vector<int> > point_indices;
+//  Eigen::MatrixXi CH;
+//  Eigen::MatrixXd CN;
+//  Eigen::VectorXd W;
+//  igl::octree(P,point_indices,CH,CN,W);
+//  Eigen::MatrixXi I;
+//  igl::knn(P,2,point_indices,CH,CN,W,I);
+//  Eigen::MatrixXd P2;
+//  igl::slice(P,I.col(1).eval(),1,P2);
+//  Eigen::VectorXd D = (P-P2).rowwise().norm();
+//  REQUIRE(D.minCoeff() > r);
+//}
+//
+//namespace blue_noise
+//{
+//  template <typename DerivedA, typename DerivedB>
+//  void assert_neq_different_sizes(
+//    const Eigen::MatrixBase<DerivedA> & A,
+//    const Eigen::MatrixBase<DerivedB> & B)
+//  {
+//    // test_common::assert_neq requires same sizes
+//    if (A.rows() == B.rows() &&
+//        A.cols() == B.cols())
+//    {
+//      test_common::assert_neq(A, B);
+//    }
+//  }
+//
+//  template<typename URBG>
+//  void test_reproduce()
+//  {
+//    Eigen::MatrixXd V;
+//    Eigen::MatrixXi F;
+//    igl::readOBJ(test_common::data_path("decimated-knight.obj"),V,F);
+//
+//    static constexpr double r = 0.1;
+//    static constexpr int seed = 0;
+//    Eigen::MatrixXd P1, Px1;
+//    {
+//      Eigen::MatrixXd B;
+//      Eigen::VectorXi I;
+//      URBG rng1(seed);
+//      igl::blue_noise(V,F,r,B,I, P1,rng1);
+//      igl::blue_noise(V,F,r,B,I,Px1,rng1);
+//    }
+//    Eigen::MatrixXd P2, Px2;
+//    {
+//      Eigen::MatrixXd B;
+//      Eigen::VectorXi I;
+//      URBG rng2(seed);
+//      igl::blue_noise(V,F,r,B,I, P2,rng2);
+//      igl::blue_noise(V,F,r,B,I,Px2,rng2);
+//    }
+//
+//    test_common::assert_eq(P1, P2);
+//    test_common::assert_eq(Px1, Px2);
+//
+//    assert_neq_different_sizes(P1, Px1);
+//    assert_neq_different_sizes(P2, Px2);
+//  }
+//}
+//
+//
+//TEST_CASE("blue_noise: minstd_rand0_reproduce", "[igl]")
+//{
+//  blue_noise::test_reproduce<std::minstd_rand0>();
+//}
+//TEST_CASE("blue_noise: minstd_rand_reproduce", "[igl]")
+//{
+//  blue_noise::test_reproduce<std::minstd_rand>();
+//}
+//TEST_CASE("blue_noise: mt19937_reproduce", "[igl]")
+//{
+//  blue_noise::test_reproduce<std::mt19937>();
+//}
+//TEST_CASE("blue_noise: mt19937_64_reproduce", "[igl]")
+//{
+//  blue_noise::test_reproduce<std::mt19937_64>();
+//}

tests/include/igl/random_points_on_mesh.cpp

@@ -0,0 +1,77 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// 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 <test_common.h>
+#include <igl/readOBJ.h>
+#include <igl/random_points_on_mesh.h>
+
+TEST_CASE("random_points_on_mesh: decimated-knight", "[igl]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::readOBJ(test_common::data_path("decimated-knight.obj"),V,F);
+  static constexpr int n = 1000;
+  Eigen::MatrixXd P;
+  {
+    Eigen::MatrixXd B;
+    Eigen::VectorXi I;
+    igl::random_points_on_mesh(n,V,F,B,I,P);
+  }
+  REQUIRE(P.rows() == n);
+}
+
+namespace random_points_on_mesh
+{
+  template<typename URBG>
+  void test_reproduce()
+  {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    igl::readOBJ(test_common::data_path("decimated-knight.obj"),V,F);
+
+    static constexpr int n = 100;
+    static constexpr int seed = 0;
+    Eigen::MatrixXd P1, Px1;
+    {
+      Eigen::MatrixXd B;
+      Eigen::VectorXi I;
+      URBG rng1(seed);
+      igl::random_points_on_mesh(n,V,F,B,I,P1 ,rng1);
+      igl::random_points_on_mesh(n,V,F,B,I,Px1,rng1);
+    }
+    Eigen::MatrixXd P2, Px2;
+    {
+      Eigen::MatrixXd B;
+      Eigen::VectorXi I;
+      URBG rng2(seed);
+      igl::random_points_on_mesh(n,V,F,B,I,P2 ,rng2);
+      igl::random_points_on_mesh(n,V,F,B,I,Px2,rng2);
+    }
+
+    test_common::assert_eq(P1, P2);
+    test_common::assert_eq(Px1, Px2);
+
+    test_common::assert_neq(P1, Px1);
+    test_common::assert_neq(P2, Px2);
+  }
+}
+
+
+TEST_CASE("random_points_on_mesh: minstd_rand0_reproduce", "[igl]")
+{
+  random_points_on_mesh::test_reproduce<std::minstd_rand0>();
+}
+TEST_CASE("random_points_on_mesh: minstd_rand_reproduce", "[igl]")
+{
+  random_points_on_mesh::test_reproduce<std::minstd_rand>();
+}
+TEST_CASE("random_points_on_mesh: mt19937_reproduce", "[igl]")
+{
+  random_points_on_mesh::test_reproduce<std::mt19937>();
+}
+TEST_CASE("random_points_on_mesh: mt19937_64_reproduce", "[igl]")
+{
+  random_points_on_mesh::test_reproduce<std::mt19937_64>();
+}

tests/include/igl/randperm.cpp

@@ -56,26 +56,6 @@ TEST_CASE("randperm: mt19937_64_reproduce_identity", "[igl]")
 {
   randperm::test_reproduce<std::mt19937_64>();
 }
-TEST_CASE("randperm: ranlux24_base_reproduce_identity", "[igl]")
-{
-  randperm::test_reproduce<std::ranlux24_base>();
-}
-TEST_CASE("randperm: ranlux48_base_reproduce_identity", "[igl]")
-{
-  randperm::test_reproduce<std::ranlux48_base>();
-}
-TEST_CASE("randperm: ranlux24_reproduce_identity", "[igl]")
-{
-  randperm::test_reproduce<std::ranlux24>();
-}
-TEST_CASE("randperm: ranlux48_reproduce_identity", "[igl]")
-{
-  randperm::test_reproduce<std::ranlux48>();
-}
-TEST_CASE("randperm: knuth_b_reproduce_identity", "[igl]")
-{
-  randperm::test_reproduce<std::knuth_b>();
-}
 TEST_CASE("randperm: default_identity", "[igl]")
 {
   int n = 100;

tutorial/717_FastWindingNumber/main.cpp

@@ -31,9 +31,8 @@ int main(int argc, char *argv[])
   Eigen::MatrixXd P,N;
   {
     Eigen::VectorXi I;
-    Eigen::SparseMatrix<double> B;
-    igl::random_points_on_mesh(10000,V,F,B,I);
-    P = B*V;
+    Eigen::MatrixXd B;
+    igl::random_points_on_mesh(10000,V,F,B,I,P);
     Eigen::MatrixXd FN;
     igl::per_face_normals(V,F,FN);
     N.resize(P.rows(),3);