blue noise tutorial entry

.gitignore

@@ -97,3 +97,4 @@ tutorial/cmake-build-debug
 tutorial/data
 tutorial/readme.html
 untitled
+scripts

include/igl/ambient_occlusion.cpp

@@ -128,6 +128,8 @@ IGL_INLINE void igl::ambient_occlusion(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::ambient_occlusion<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<double, -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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::ambient_occlusion<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

include/igl/barycenter.cpp

@@ -33,19 +33,12 @@ IGL_INLINE void igl::barycenter(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(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, 0, -1, 3> >&);
-// generated by autoexplicit.sh
 template void igl::barycenter<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -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::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(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, 4, 0, -1, 4> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);

include/igl/barycentric_interpolation.cpp

@@ -0,0 +1,42 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "barycentric_interpolation.h"
+#include "parallel_for.h"
+
+template <
+  typename DerivedD,
+  typename DerivedF,
+  typename DerivedB,
+  typename DerivedI,
+  typename DerivedX>
+IGL_INLINE void igl::barycentric_interpolation(
+  const Eigen::MatrixBase<DerivedD> & D,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedB> & B,
+  const Eigen::MatrixBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  assert(B.rows() == I.size());
+  assert(F.cols() == B.cols());
+  X.setZero(B.rows(),D.cols());
+  // should use parallel_for
+  //for(int i = 0;i<X.rows();i++)
+  parallel_for(X.rows(),[&X,&B,&D,&F,&I](const int i)
+  {
+    for(int j = 0;j<F.cols();j++)
+    {
+      X.row(i) += B(i,j) * D.row(F(I(i),j));
+    }
+  },1000);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::barycentric_interpolation<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::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> >&);
+#endif

include/igl/barycentric_interpolation.h

@@ -0,0 +1,42 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_BARYCENTRIC_INTERPOLATION_H
+#define IGL_BARYCENTRIC_INTERPOLATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Interpolate data on a triangle mesh using barycentric coordinates 
+  //
+  // Inputs:
+  //   D  #D by dim list of per-vertex data
+  //   F  #F by 3 list of triangle indices
+  //   B  #X by 3 list of barycentric corodinates
+  //   I  #X list of triangle indices
+  // Outputs:
+  //   X  #X by dim list of interpolated data
+  template <
+    typename DerivedD,
+    typename DerivedF,
+    typename DerivedB,
+    typename DerivedI,
+    typename DerivedX>
+  IGL_INLINE void barycentric_interpolation(
+    const Eigen::MatrixBase<DerivedD> & D,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const Eigen::MatrixBase<DerivedI> & I,
+    Eigen::PlainObjectBase<DerivedX> & X);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "barycentric_interpolation.cpp"
+#endif
+
+#endif
+

include/igl/blue_noise.cpp

@@ -0,0 +1,375 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "blue_noise.h"
+#include "doublearea.h"
+#include "random_points_on_mesh.h"
+#include "slice.h"
+#include "sortrows.h"
+#include "PI.h"
+#include "get_seconds.h"
+#include <unordered_map>
+#include <algorithm>
+#include <vector>
+
+namespace igl
+{
+  // It is very important that we use 64bit keys to avoid out of bounds (easy to
+  // get to happen with dense samplings (e.g., r = 0.0005*bbd)
+  typedef int64_t BlueNoiseKeyType;
+}
+
+// Helper functions
+namespace igl
+{
+  // Should probably find and replace with less generic name
+  //
+  // Map 3D subscripts (x,y,z) to unique index (return value)
+  //
+  // Inputs:
+  //   w  side length of w×w×w integer cube lattice
+  //   x  subscript along x direction
+  //   y  subscript along y direction
+  //   z  subscript along z direction
+  // Returns index value
+  //
+  inline BlueNoiseKeyType key(
+    const BlueNoiseKeyType w, // pass by copy --> int64_t so that multiplication is OK
+    const BlueNoiseKeyType x, // pass by copy --> int64_t so that multiplication is OK
+    const BlueNoiseKeyType y, // pass by copy --> int64_t so that multiplication is OK
+    const BlueNoiseKeyType z) // pass by copy --> int64_t so that multiplication is OK
+  {
+    return x+w*(y+w*z);
+  }
+  // Determine if a query candidate at position X.row(i) is too close to already
+  // selected sites (stored in S).
+  //
+  // Inputs:
+  //   X  #X by 3 list of raw candidate positions
+  //   Xs  #Xs by 3 list of corresponding integer cell subscripts
+  //   i  index of candidate in question
+  //   S   map from cell index to index into X of selected candidate (or -1 if
+  //     cell is currently empty)
+  //   rr  Poisson disk radius squared
+  //   w  side length of w×w×w integer cube lattice (into which Xs subscripts)
+  template <
+    typename DerivedX,
+    typename DerivedXs>
+  inline bool far_enough(
+    const Eigen::MatrixBase<DerivedX> & X,
+    const Eigen::MatrixBase<DerivedXs> & Xs,
+    const std::unordered_map<BlueNoiseKeyType,int> & S,
+    const double & rr,
+    const int & w,
+    const int i)
+  {
+    const int xi = Xs(i,0);
+    const int yi = Xs(i,1);
+    const int zi = Xs(i,2);
+    BlueNoiseKeyType k = key(w,xi,yi,zi);
+    int g = 2; // ceil(r/s)
+    for(int x = std::max(xi-g,0);x<=std::min(xi+g,w-1);x++)
+    for(int y = std::max(yi-g,0);y<=std::min(yi+g,w-1);y++)
+    for(int z = std::max(zi-g,0);z<=std::min(zi+g,w-1);z++)
+    {
+      if(x!=xi || y!=yi || z!=zi)
+      {
+        const BlueNoiseKeyType nk = key(w,x,y,z);
+        // have already selected from this cell
+        const auto Siter = S.find(nk);
+        if(Siter !=S.end() && Siter->second >= 0)
+        {
+          const int ni = Siter->second;
+          // too close
+          if( (X.row(i)-X.row(ni)).squaredNorm() < rr)
+          {
+            return false;
+          }
+        }
+      }
+    }
+    return true;
+  }
+  // Try to activate a candidate in a given cell
+  //
+  // Inputs:
+  //   X  #X by 3 list of raw candidate positions
+  //   Xs  #Xs by 3 list of corresponding integer cell subscripts
+  //   rr  Poisson disk radius squared
+  //   w  side length of w×w×w integer cube lattice (into which Xs subscripts)
+  //   nk  index of cell in which we'd like to activate a candidate
+  //   M   map from cell index to list of candidates
+  //   S   map from cell index to index into X of selected candidate (or -1 if
+  //     cell is currently empty)
+  //   active   list of indices into X of active candidates
+  // Outputs:
+  //   M  visited candidates deemed too close to already selected points are
+  //      removed
+  //   S  updated to reflect activated point (if successful)
+  //   active   updated to reflect activated point (if successful)
+  // Returns true iff activation was successful
+  template <
+    typename DerivedX,
+    typename DerivedXs>
+  inline bool activate(
+    const Eigen::MatrixBase<DerivedX> & X,
+    const Eigen::MatrixBase<DerivedXs> & Xs,
+    const double & rr,
+    const int & i,
+    const int & w,
+    const BlueNoiseKeyType & nk,
+    std::unordered_map<BlueNoiseKeyType,std::vector<int> > & M,
+    std::unordered_map<BlueNoiseKeyType,int> & S,
+    std::vector<int> & active)
+  {
+    auto & Mvec = M.find(nk)->second;
+    auto miter = Mvec.begin();
+    while(miter != Mvec.end())
+    {
+      const int mi = *miter;
+      // mi is our candidate sample. Is it far enough from all existing
+      // samples?
+      if(i>=0 && (X.row(i)-X.row(mi)).squaredNorm() > 4.*rr)
+      {
+        // too far skip (reject)
+       miter++;
+      } else if(far_enough(X,Xs,S,rr,w,mi))
+      {
+        active.push_back(mi);
+        S.find(nk)->second = mi;
+        //printf("  found %d\n",mi);
+        return true;
+      }else
+      {
+        // remove forever (instead of incrementing we swap and eat from the
+        // back)
+        //std::swap(*miter,Mvec.back());
+        *miter = Mvec.back();
+        Mvec.pop_back();
+      }
+    }
+    return false;
+  }
+
+  template <
+    typename DerivedX,
+    typename DerivedXs>
+  inline bool step(
+    const Eigen::MatrixBase<DerivedX> & X,
+    const Eigen::MatrixBase<DerivedXs> & Xs,
+    const double & rr,
+    const int & w,
+    std::unordered_map<BlueNoiseKeyType,std::vector<int> > & M,
+    std::unordered_map<BlueNoiseKeyType,int> & S,
+    std::vector<int> & active,
+    std::vector<int> & collected
+    )
+  {
+    //considered.clear();
+    if(active.size() == 0) return false;
+    // random entry
+    const int e = rand() % active.size();
+    const int i = active[e];
+    //printf("%d\n",i);
+    const int xi = Xs(i,0);
+    const int yi = Xs(i,1);
+    const int zi = Xs(i,2);
+    //printf("%d %d %d - %g %g %g\n",xi,yi,zi,X(i,0),X(i,1),X(i,2));
+    // cell indices of neighbors
+    int g = 4;
+    std::vector<BlueNoiseKeyType> N;N.reserve((1+g*1)^3-1);
+    for(int x = std::max(xi-g,0);x<=std::min(xi+g,w-1);x++)
+    for(int y = std::max(yi-g,0);y<=std::min(yi+g,w-1);y++)
+    for(int z = std::max(zi-g,0);z<=std::min(zi+g,w-1);z++)
+    {
+      if(x!=xi || y!=yi || z!=zi)
+      {
+        //printf("  %d %d %d\n",x,y,z);
+        const BlueNoiseKeyType nk = key(w,x,y,z);
+        // haven't yet selected from this cell?
+        const auto Siter = S.find(nk);
+        if(Siter !=S.end() && Siter->second < 0)
+        {
+          assert(M.find(nk) != M.end());
+          N.emplace_back(nk);
+        }
+      }
+    }
+        //printf("  --------\n");
+    // randomize order: this might be a little paranoid...
+    std::random_shuffle(std::begin(N), std::end(N));
+    bool found = false;
+    for(const BlueNoiseKeyType & nk : N)
+    {
+      assert(M.find(nk) != M.end());
+      if(activate(X,Xs,rr,i,w,nk,M,S,active))
+      {
+        found = true;
+        break;
+      }
+    }
+    if(!found)
+    {
+      // remove i from active list
+      // https://stackoverflow.com/a/60765833/148668
+      collected.push_back(i);
+      //printf("  before: "); for(const int j : active) { printf("%d ",j); } printf("\n");
+      std::swap(active[e], active.back());
+      //printf("  after : "); for(const int j : active) { printf("%d ",j); } printf("\n");
+      active.pop_back();
+      //printf("  removed %d\n",i);
+    }
+    //printf("  active: "); for(const int j : active) { printf("%d ",j); } printf("\n");
+    return true;
+  }
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedB,
+  typename DerivedFI,
+  typename DerivedP>
+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)
+{
+  const auto & tictoc = []() -> double
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  tictoc();
+  typedef typename DerivedV::Scalar Scalar;
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  // float+RowMajor is faster...
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,3,Eigen::RowMajor> MatrixX3S;
+  assert(V.cols() == 3 && "Only 3D embeddings allowed");
+  // minimum radius
+  const Scalar min_r = r;
+  // cell size based on 3D distance
+  // It works reasonably well (but is probably biased to use s=2*r/√3 here and
+  // g=1 in the outer loop below.
+  //
+  // One thing to try would be to store a list in S (rather than a single point)
+  // or equivalently a mask over M and just use M as a generic spatial hash
+  // (with arbitrary size) and then tune its size (being careful to make g a
+  // function of r and s; and removing the `if S=-1 checks`)
+  const Scalar s = r/sqrt(3.0);
+
+  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 = 
+    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;
+  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 = 
+    ((X.rowwise()-X.colwise().minCoeff())/s).template cast<int>();
+  const int w = Xs.maxCoeff()+1;
+  {
+    Eigen::VectorXi I;
+    igl::sortrows(decltype(Xs)(Xs),true,Xs,I);
+    igl::slice(decltype(X)(X),I,1,X);
+    // These two could be spun off in their own thread.
+    igl::slice(decltype(XB)(XB),I,1,XB);
+    igl::slice(decltype(XFI)(XFI),I,1,XFI);
+  }
+  // Initialization
+  std::unordered_map<BlueNoiseKeyType,std::vector<int> > M;
+  std::unordered_map<BlueNoiseKeyType, int > S;
+  // attempted to seed
+  std::unordered_map<BlueNoiseKeyType, int > A;
+  // Q: Too many?
+  // A: Seems to help though.
+  M.reserve(Xs.rows());
+  S.reserve(Xs.rows());
+  for(int i = 0;i<Xs.rows();i++)
+  {
+    BlueNoiseKeyType k = key(w,Xs(i,0),Xs(i,1),Xs(i,2)); 
+    const auto Miter = M.find(k);
+    if(Miter  == M.end())
+    {
+      M.insert({k,{i}});
+    }else
+    {
+      Miter->second.push_back(i);
+    }
+    S.emplace(k,-1);
+    A.emplace(k,false);
+  }
+
+  std::vector<int> active;
+  // precompute r²
+  // Q: is this necessary? 
+  const double rr = r*r;
+  std::vector<int> collected;
+  collected.reserve(2.0*expected_number_of_points);
+  //printf("  init: %g secs\n",tictoc());
+
+  auto Mouter = M.begin();
+  // Just take the first point as the initial seed
+  const auto initialize = [&]()->bool
+  {
+    while(true)
+    {
+      if(Mouter == M.end())
+      {
+        return false;
+      }
+      const BlueNoiseKeyType k = Mouter->first;
+      // Haven't placed in this cell yet
+      if(S[k]<0)
+      {
+        if(activate(X,Xs,rr,-1,w,k,M,S,active)) return true;
+      }
+      Mouter++;
+    }
+    assert(false && "should not be reachable.");
+  };
+
+  // important if mesh contains many connected components
+  while(initialize())
+  {
+    while(active.size()>0)
+    {
+      step(X,Xs,rr,w,M,S,active,collected);
+    }
+  }
+  //printf(" while: %g secs\n",tictoc());
+  {
+    const int n = collected.size();
+    P.resize(n,3);
+    B.resize(n,3);
+    FI.resize(n);
+    for(int i = 0;i<n;i++)
+    {
+      const int c = collected[i];
+      P.row(i) = X.row(c).template cast<typename DerivedP::Scalar>();
+      B.row(i) = XB.row(c).template cast<typename DerivedB::Scalar>();
+      FI(i) = XFI(c);
+    }
+  }
+  //printf("collec: %g secs\n",tictoc());
+}
+
+#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> >&);
+#endif

include/igl/blue_noise.h

@@ -0,0 +1,48 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_BLUE_NOISE_H
+#define IGL_BLUE_NOISE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // "Fast Poisson Disk Sampling in Arbitrary Dimensions" [Bridson 2007]
+  //
+  // For very dense samplings this is faster than (up to 2x) cyCodeBase's
+  // implementation of "Sample Elimination for Generating Poisson Disk Sample
+  // Sets" [Yuksel 2015]. YMMV
+  //
+  // Inputs:
+  //    V  #V by dim list of mesh vertex positions
+  //    F  #F by 3 list of mesh triangle indices into rows of V
+  //    r  Poisson disk radius (evaluated according to Euclidean distance on V)
+  //   B  n by 3 list of barycentric coordinates, ith row are coordinates of
+  //     ith sampled point in face FI(i)
+  //   FI  n list of indices into F 
+  //   P  n by dim list of sample positions.
+  // See also: random_points_on_mesh
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedB,
+    typename DerivedFI,
+    typename DerivedP>
+  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);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "blue_noise.cpp"
+#endif
+
+#endif

include/igl/opengl/ViewerData.cpp

@@ -369,6 +369,24 @@ IGL_INLINE void igl::opengl::ViewerData::set_edges(
   dirty |= MeshGL::DIRTY_OVERLAY_LINES;
 }
 
+IGL_INLINE void igl::opengl::ViewerData::set_edges_from_vector_field(
+  const Eigen::MatrixXd& P, 
+  const Eigen::MatrixXd& V, 
+  const Eigen::MatrixXd& C)
+{
+  assert(P.rows() == V.rows());
+  Eigen::MatrixXi E(P.rows(),2);
+  const Eigen::MatrixXd PV = 
+    (Eigen::MatrixXd(P.rows()+V.rows(),3)<<P,P+V).finished();
+  for(int i = 0;i<P.rows();i++)
+  {
+    E(i,0) = i;
+    E(i,1) = i+P.rows();
+  }
+  const Eigen::MatrixXd CC = C.replicate<2,1>();
+  set_edges(PV,E, C.rows() == 1?C:C.replicate<2,1>());
+}
+
 IGL_INLINE void igl::opengl::ViewerData::add_edges(const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C)
 {
   Eigen::MatrixXd P1_temp,P2_temp;

include/igl/opengl/ViewerData.h

@@ -156,6 +156,11 @@ public:
   // Alec: This is very confusing. Why does add_edges have a different API from
   // set_edges?
   IGL_INLINE void add_edges (const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C);
+  // Sets edges given a list of points and eminating vectors
+  IGL_INLINE void set_edges_from_vector_field(
+    const Eigen::MatrixXd& P, 
+    const Eigen::MatrixXd& V, 
+    const Eigen::MatrixXd& C);
 
   // Clear the edge data
   IGL_INLINE void clear_edges();

include/igl/random_points_on_mesh.cpp

@@ -106,6 +106,8 @@ IGL_INLINE void igl::random_points_on_mesh(
 
 #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, 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> >&);

include/igl/slice.cpp

@@ -189,6 +189,10 @@ IGL_INLINE DerivedX igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::slice<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<float, -1, 3, 1, -1, 3> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<float, -1, 3, 1, -1, 3>&);
+// generated by autoexplicit.sh
+template void igl::slice<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::Matrix<double, -1, 3, 1, -1, 3> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 3, 1, -1, 3>&);
+// generated by autoexplicit.sh
 template void igl::slice<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<double, -1, 1, 0, -1, 1> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
 // generated by autoexplicit.sh
 template void igl::slice<Eigen::Array<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Array<bool, -1, 1, 0, -1, 1> >(Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Array<bool, -1, 1, 0, -1, 1>&);

tutorial/720_BlueNoise/CMakeLists.txt

@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME}_bin main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw tutorials)

tutorial/720_BlueNoise/main.cpp

@@ -0,0 +1,105 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/blue_noise.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/barycentric_interpolation.h>
+#include <igl/blue_noise.h>
+#include <igl/doublearea.h>
+#include <igl/random_points_on_mesh.h>
+
+
+int main(int argc, char *argv[])
+{
+
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(
+    argc>1?argv[1]: TUTORIAL_SHARED_PATH "/elephant.obj",V,F);
+  Eigen::MatrixXd N;
+  igl::per_vertex_normals(V,F,N);
+
+  Eigen::MatrixXd P_blue;
+  Eigen::MatrixXd N_blue;
+  Eigen::VectorXd A_blue;
+  Eigen::MatrixXd C_blue;
+  {
+    const int n_desired = argc>2?atoi(argv[2]):50000;
+    // Heuristic to  determine radius from desired number 
+    const double r = [&V,&F](const int n)
+    {
+      Eigen::VectorXd A;
+      igl::doublearea(V,F,A);
+      return sqrt(((A.sum()*0.5/(n*0.6162910373))/M_PI));
+    }(n_desired);
+    Eigen::MatrixXd B;
+    Eigen::VectorXi I;
+    igl::blue_noise(V,F,r,B,I,P_blue);
+    igl::barycentric_interpolation(N,F,B,I,N_blue);
+    N_blue.rowwise().normalize();
+  }
+  Eigen::MatrixXd P_white;
+  Eigen::MatrixXd N_white;
+  Eigen::VectorXd A_white;
+  Eigen::MatrixXd C_white;
+  {
+    Eigen::MatrixXd B;
+    Eigen::VectorXi I;
+    igl::random_points_on_mesh(P_blue.rows(),V,F,B,I,P_white);
+    igl::barycentric_interpolation(N,F,B,I,N_white);
+    N_white.rowwise().normalize();
+  }
+
+  // Color
+  Eigen::RowVector3d C(1,1,1);
+
+  // Plot the mesh
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V,F);
+  viewer.data().show_lines = false;
+  viewer.data().show_faces = false;
+  viewer.data().point_size = 4;
+
+  bool use_blue = true;
+  const auto update = [&]()
+  {
+    const Eigen::RowVector3d orange(1.0,0.7,0.2);
+    if(use_blue)
+    {
+      viewer.data().set_points(P_blue,Eigen::RowVector3d(0.3,0.4,1.0));
+      viewer.data().set_edges_from_vector_field(P_blue,N_blue,orange);
+    }else
+    {
+      viewer.data().set_points(P_white,Eigen::RowVector3d(1,1,1));
+      viewer.data().set_edges_from_vector_field(P_white,N_white,orange);
+    }
+  };
+  update();
+
+  viewer.callback_key_pressed = 
+    [&](igl::opengl::glfw::Viewer &,unsigned char key,int)->bool
+  {
+    switch(key)
+    {
+      case ' ':
+      {
+        use_blue = !use_blue;
+        update();
+        return true;
+      }
+    }
+    return false;
+  };
+  std::cout<<R"(
+[space]  toggle between blue and white noise samplings
+)";
+
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -164,4 +164,5 @@ if(TUTORIALS_CHAPTER7)
   endif()
   add_subdirectory("718_IterativeClosestPoint")
   add_subdirectory("719_ExplodedView")
+  add_subdirectory("720_BlueNoise")
 endif()