performance improvements for decimate

include/igl/always_try_never_care.cpp

@@ -0,0 +1,36 @@
+#include "always_try_never_care.h"
+
+IGL_INLINE void igl::always_try_never_care(
+    decimate_pre_collapse_func  & always_try,
+    decimate_post_collapse_func & never_care)
+{
+  always_try = [](
+    const Eigen::MatrixXd &                             ,/*V*/
+    const Eigen::MatrixXi &                             ,/*F*/
+    const Eigen::MatrixXi &                             ,/*E*/
+    const Eigen::VectorXi &                             ,/*EMAP*/
+    const Eigen::MatrixXi &                             ,/*EF*/
+    const Eigen::MatrixXi &                             ,/*EI*/
+    const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+    const Eigen::VectorXi &                             ,/*EQ*/
+    const Eigen::MatrixXd &                             ,/*C*/
+    const int                                            /*e*/
+    ) -> bool { return true;};
+  never_care = [](
+    const Eigen::MatrixXd &                             ,/*V*/
+    const Eigen::MatrixXi &                             ,/*F*/
+    const Eigen::MatrixXi &                             ,/*E*/
+    const Eigen::VectorXi &                             ,/*EMAP*/
+    const Eigen::MatrixXi &                             ,/*EF*/
+    const Eigen::MatrixXi &                             ,/*EI*/
+    const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+    const Eigen::VectorXi &                             ,/*EQ*/
+    const Eigen::MatrixXd &                             ,/*C*/
+    const int                                           ,/*e*/
+    const int                                           ,/*e1*/
+    const int                                           ,/*e2*/
+    const int                                           ,/*f1*/
+    const int                                           ,/*f2*/
+    const bool                                           /*collapsed*/
+    )-> void { };
+}

include/igl/always_try_never_care.h

@@ -0,0 +1,20 @@
+#ifndef IGL_ALWAYS_TRY_NEVER_CARE_H
+#define IGL_ALWAYS_TRY_NEVER_CARE_H
+#include "igl_inline.h"
+#include "decimate.h" // for decimate_*_func types
+namespace igl
+{
+  // Outputs:
+  //   always_try  function that always returns true
+  //   never_care  fuction that is always a no-op
+  IGL_INLINE void always_try_never_care(
+    decimate_pre_collapse_func  & always_try,
+    decimate_post_collapse_func & never_care);
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "always_try_never_care.cpp"
+#endif
+
+#endif 
+

include/igl/circulation.cpp

@@ -67,3 +67,77 @@ IGL_INLINE void igl::circulation(
   std::vector<int> N = circulation(e,ccw,EMAP,EF,EI);
   igl::list_to_matrix(N,vN);
 }
+
+IGL_INLINE void igl::circulation(
+  const int e,
+  const bool ccw,
+  const Eigen::MatrixXi & F,
+  const Eigen::VectorXi & EMAP,
+  const Eigen::MatrixXi & EF,
+  const Eigen::MatrixXi & EI,
+  /*std::vector<int> & Ne,*/
+  std::vector<int> & Nv,
+  std::vector<int> & Nf)
+{
+  //
+  // for e --> (bf) and ccw=true
+  //
+  //     c---d
+  //    / \ / \
+  //   a---b-e-f
+  //    \ / \ /
+  //     g---h
+  //
+  //  // (might start with {bhf} depending on edge)
+  //  Ne = […] -> [fd db dc cb ca ab ag gb gh hb hf fb]
+  //              {upto cylic order}
+  //  Nf = […] -> [{bfd}, {bdc}, {bca}, {bag}, {bgh}, {bhf}]
+  //  Nv = [d c a g h f]
+  //
+  // prepare output
+  //Ne.clear();Ne.reserve(2*10);
+  Nv.clear();Nv.reserve(10);
+  Nf.clear();Nf.reserve(10);
+  const int m = EMAP.size()/3;
+  assert(m*3 == EMAP.size());
+  const auto & step = [&](
+    const int e, 
+    const int ff,
+    int & ne, 
+    //int & re,
+    int & rv,
+    int & nf)
+  {
+    assert((EF(e,1) == ff || EF(e,0) == ff) && "e should touch ff");
+    //const int fside = EF(e,1)==ff?1:0;
+    const int nside = EF(e,0)==ff?1:0;
+    const int nv = EI(e,nside);
+    // get next face
+    nf = EF(e,nside);
+    // get next edge 
+    const int dir = ccw?-1:1;
+    rv = F(nf,nv);
+    ne = EMAP(nf+m*((nv+dir+3)%3));
+    //re = EMAP(nf+m*((nv+2*dir+3)%3));
+  };
+  // Always start with first face (ccw in step will be sure to turn right
+  // direction)
+  const int f0 = EF(e,0);
+  int fi = f0;
+  int ei = e;
+  while(true)
+  {
+    int re,rv;
+    step(ei,fi,ei/*,re*/,rv,fi);
+    Nf.push_back(fi);
+    //Ne.push_back(re);
+    //Ne.push_back(ei);
+    Nv.push_back(rv);
+    // back to start?
+    if(fi == f0)
+    {
+      assert(ei == e);
+      break;
+    }
+  }
+}

include/igl/circulation.h

@@ -43,6 +43,20 @@ namespace igl
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI,
     Eigen::VectorXi & vN);
+  // Outputs:
+  ////   Ne  2*#Nf list of indices into E of "next" rim-spoke-rim-spoke-...
+  //   Nv  #Nv list of "next" vertex indices
+  //   Nf  #Nf list of face indices
+  IGL_INLINE void circulation(
+    const int e,
+    const bool ccw,
+    const Eigen::MatrixXi & F,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
+    /*std::vector<int> & Ne,*/
+    std::vector<int> & Nv,
+    std::vector<int> & Nf);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/collapse_edge.cpp

@@ -8,6 +8,7 @@
 #include "collapse_edge.h"
 #include "circulation.h"
 #include "edge_collapse_is_valid.h"
+#include "always_try_never_care.h"
 #include <vector>
 
 IGL_INLINE bool igl::collapse_edge(
@@ -19,6 +20,32 @@ IGL_INLINE bool igl::collapse_edge(
   Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI,
+  int & e1,
+  int & e2,
+  int & f1,
+  int & f2)
+{
+  std::vector<int> /*Nse,*/Nsf,Nsv;
+  circulation(e, true,F,EMAP,EF,EI,/*Nse,*/Nsv,Nsf);
+  std::vector<int> /*Nde,*/Ndf,Ndv;
+  circulation(e, false,F,EMAP,EF,EI,/*Nde,*/Ndv,Ndf);
+  return collapse_edge(
+    e,p,Nsv,Nsf,Ndv,Ndf,V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+}
+
+IGL_INLINE bool igl::collapse_edge(
+  const int e,
+  const Eigen::RowVectorXd & p,
+  /*const*/ std::vector<int> & Nsv,
+  const std::vector<int> & Nsf,
+  /*const*/ std::vector<int> & Ndv,
+  const std::vector<int> & Ndf,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI,
   int & a_e1,
   int & a_e2,
   int & a_f1,
@@ -34,17 +61,19 @@ IGL_INLINE bool igl::collapse_edge(
   const int s = eflip?E(e,1):E(e,0);
   const int d = eflip?E(e,0):E(e,1);
 
-  if(!edge_collapse_is_valid(e,F,E,EMAP,EF,EI))
+  if(!edge_collapse_is_valid(Nsv,Ndv))
   {
     return false;
   }
 
+  // OVERLOAD: caller may have just computed this
+  //
   // Important to grab neighbors of d before monkeying with edges
-  const std::vector<int> nV2Fd = circulation(e,!eflip,EMAP,EF,EI);
+  const std::vector<int> & nV2Fd = (!eflip ? Nsf : Ndf);
 
   // The following implementation strongly relies on s<d
   assert(s<d && "s should be less than d");
-  // move source and destination to midpoint
+  // move source and destination to placement
   V.row(s) = p;
   V.row(d) = p;
 
@@ -112,31 +141,40 @@ IGL_INLINE bool igl::collapse_edge(
   // make sense.
   //
   // Could actually skip first and last, since those are always the two
-  // collpased faces.
-  for(auto f : nV2Fd)
+  // collpased faces. Nah, this is handled by (F(f,v) == d)
+  //
+  // Don't attempt to use Nde,Nse here because EMAP has changed
   {
-    for(int v = 0;v<3;v++)
+    int p1 = -1;
+    for(auto f : nV2Fd)
     {
-      if(F(f,v) == d)
+      for(int v = 0;v<3;v++)
       {
-        const int flip1 = (EF(EMAP(f+m*((v+1)%3)),0)==f)?1:0;
-        const int flip2 = (EF(EMAP(f+m*((v+2)%3)),0)==f)?0:1;
-        assert(
-          E(EMAP(f+m*((v+1)%3)),flip1) == d ||
-          E(EMAP(f+m*((v+1)%3)),flip1) == s);
-        E(EMAP(f+m*((v+1)%3)),flip1) = s;
-        assert(
-          E(EMAP(f+m*((v+2)%3)),flip2) == d ||
-          E(EMAP(f+m*((v+2)%3)),flip2) == s);
-        E(EMAP(f+m*((v+2)%3)),flip2) = s;
-        F(f,v) = s;
-        break;
+        if(F(f,v) == d)
+        {
+          const int e1 = EMAP(f+m*((v+1)%3));
+          const int flip1 = (EF(e1,0)==f)?1:0;
+          assert( E(e1,flip1) == d || E(e1,flip1) == s);
+          E(e1,flip1) = s;
+          const int e2 = EMAP(f+m*((v+2)%3));
+          // Skip if we just handled this edge (claim: this will be all except
+          // for the first non-trivial face)
+          if(e2 != p1)
+          {
+            const int flip2 = (EF(e2,0)==f)?0:1;
+            assert( E(e2,flip2) == d || E(e2,flip2) == s);
+            E(e2,flip2) = s;
+          }
+
+          F(f,v) = s;
+          p1 = e1;
+          break;
+        }
       }
     }
   }
   // Finally, "remove" this edge and its information
   kill_edge(e);
-
   return true;
 }
 
@@ -155,168 +193,61 @@ IGL_INLINE bool igl::collapse_edge(
 }
 
 IGL_INLINE bool igl::collapse_edge(
-  const std::function<void(
-    const int,
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    double &,
-    Eigen::RowVectorXd &)> & cost_and_placement,
+  const decimate_cost_and_placement_func & cost_and_placement,
   Eigen::MatrixXd & V,
   Eigen::MatrixXi & F,
   Eigen::MatrixXi & E,
   Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI,
-  std::set<std::pair<double,int> > & Q,
-  std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+  igl::min_heap< std::tuple<double,int,int> > & Q,
+  Eigen::VectorXi & EQ,
   Eigen::MatrixXd & C)
 {
   int e,e1,e2,f1,f2;
-  const auto always_try = [](
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi &                                         ,/*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
-    const int                                                        /*e*/
-    ) -> bool { return true;};
-  const auto never_care = [](
-    const Eigen::MatrixXd &                                         ,   /*V*/
-    const Eigen::MatrixXi &                                         ,   /*F*/
-    const Eigen::MatrixXi &                                         ,   /*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,  /*EF*/
-    const Eigen::MatrixXi &                                         ,  /*EI*/
-    const std::set<std::pair<double,int> > &                        ,   /*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-    const Eigen::MatrixXd &                                         ,   /*C*/
-    const int                                                       ,   /*e*/
-    const int                                                       ,  /*e1*/
-    const int                                                       ,  /*e2*/
-    const int                                                       ,  /*f1*/
-    const int                                                       ,  /*f2*/
-    const bool                                                  /*collapsed*/
-    )-> void { };
+  decimate_pre_collapse_func always_try;
+  decimate_post_collapse_func never_care;
+  always_try_never_care(always_try,never_care);
   return 
     collapse_edge(
       cost_and_placement,always_try,never_care,
-      V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+      V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2);
 }
 
 IGL_INLINE bool igl::collapse_edge(
-  const std::function<void(
-    const int,
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    double &,
-    Eigen::RowVectorXd &)> & cost_and_placement,
-  const std::function<bool(
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi &                                         ,/*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
-    const int                                                        /*e*/
-    )> & pre_collapse,
-  const std::function<void(
-    const Eigen::MatrixXd &                                         ,   /*V*/
-    const Eigen::MatrixXi &                                         ,   /*F*/
-    const Eigen::MatrixXi &                                         ,   /*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,  /*EF*/
-    const Eigen::MatrixXi &                                         ,  /*EI*/
-    const std::set<std::pair<double,int> > &                        ,   /*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-    const Eigen::MatrixXd &                                         ,   /*C*/
-    const int                                                       ,   /*e*/
-    const int                                                       ,  /*e1*/
-    const int                                                       ,  /*e2*/
-    const int                                                       ,  /*f1*/
-    const int                                                       ,  /*f2*/
-    const bool                                                  /*collapsed*/
-    )> & post_collapse,
+  const decimate_cost_and_placement_func & cost_and_placement,
+  const decimate_pre_collapse_func       & pre_collapse,
+  const decimate_post_collapse_func      & post_collapse,
   Eigen::MatrixXd & V,
   Eigen::MatrixXi & F,
   Eigen::MatrixXi & E,
   Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI,
-  std::set<std::pair<double,int> > & Q,
-  std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+  igl::min_heap< std::tuple<double,int,int> > & Q,
+  Eigen::VectorXi & EQ,
   Eigen::MatrixXd & C)
 {
   int e,e1,e2,f1,f2;
   return 
     collapse_edge(
       cost_and_placement,pre_collapse,post_collapse,
-      V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+      V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2);
 }
 
 
 IGL_INLINE bool igl::collapse_edge(
-  const std::function<void(
-    const int,
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    double &,
-    Eigen::RowVectorXd &)> & cost_and_placement,
-  const std::function<bool(
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi &                                         ,/*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
-    const int                                                        /*e*/
-    )> & pre_collapse,
-  const std::function<void(
-    const Eigen::MatrixXd &                                         ,   /*V*/
-    const Eigen::MatrixXi &                                         ,   /*F*/
-    const Eigen::MatrixXi &                                         ,   /*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,  /*EF*/
-    const Eigen::MatrixXi &                                         ,  /*EI*/
-    const std::set<std::pair<double,int> > &                        ,   /*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-    const Eigen::MatrixXd &                                         ,   /*C*/
-    const int                                                       ,   /*e*/
-    const int                                                       ,  /*e1*/
-    const int                                                       ,  /*e2*/
-    const int                                                       ,  /*f1*/
-    const int                                                       ,  /*f2*/
-    const bool                                                  /*collapsed*/
-    )> & post_collapse,
+  const decimate_cost_and_placement_func & cost_and_placement,
+  const decimate_pre_collapse_func       & pre_collapse,
+  const decimate_post_collapse_func      & post_collapse,
   Eigen::MatrixXd & V,
   Eigen::MatrixXi & F,
   Eigen::MatrixXi & E,
   Eigen::VectorXi & EMAP,
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI,
-  std::set<std::pair<double,int> > & Q,
-  std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+  igl::min_heap< std::tuple<double,int,int> > & Q,
+  Eigen::VectorXi & EQ,
   Eigen::MatrixXd & C,
   int & e,
   int & e1,
@@ -325,43 +256,80 @@ IGL_INLINE bool igl::collapse_edge(
   int & f2)
 {
   using namespace Eigen;
-  if(Q.empty())
+  using namespace igl;
+  std::tuple<double,int,int> p;
+  while(true)
   {
-    // no edges to collapse
-    return false;
-  }
-  std::pair<double,int> p = *(Q.begin());
-  if(p.first == std::numeric_limits<double>::infinity())
-  {
-    // min cost edge is infinite cost
-    return false;
+    // Check if Q is empty
+    if(Q.empty())
+    {
+      // no edges to collapse
+      return false;
+    }
+    // pop from Q
+    p = Q.top();
+    if(std::get<0>(p) == std::numeric_limits<double>::infinity())
+    {
+      // min cost edge is infinite cost
+      return false;
+    }
+    Q.pop();
+    e = std::get<1>(p);
+    // Check if matches timestamp
+    if(std::get<2>(p) == EQ(e))
+    {
+      break;
+    }
+    // must be stale or dead.
+    assert(std::get<2>(p)  < EQ(e) || EQ(e) == -1);
+    // try again.
   }
-  Q.erase(Q.begin());
-  e = p.second;
-  Qit[e] = Q.end();
-  std::vector<int> N  = circulation(e, true,EMAP,EF,EI);
-  std::vector<int> Nd = circulation(e,false,EMAP,EF,EI);
-  N.insert(N.begin(),Nd.begin(),Nd.end());
+
+  // Why is this computed up here?
+  // If we just need original face neighbors of edge, could we gather that more
+  // directly than gathering face neighbors of each vertex?
+  std::vector<int> /*Nse,*/Nsf,Nsv;
+  circulation(e, true,F,EMAP,EF,EI,/*Nse,*/Nsv,Nsf);
+  std::vector<int> /*Nde,*/Ndf,Ndv;
+  circulation(e, false,F,EMAP,EF,EI,/*Nde,*/Ndv,Ndf);
+
+
   bool collapsed = true;
-  if(pre_collapse(V,F,E,EMAP,EF,EI,Q,Qit,C,e))
+  if(pre_collapse(V,F,E,EMAP,EF,EI,Q,EQ,C,e))
   {
-    collapsed = collapse_edge(e,C.row(e),V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+    collapsed = collapse_edge(
+      e,C.row(e),
+      Nsv,Nsf,Ndv,Ndf,
+      V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
   }else
   {
     // Aborted by pre collapse callback
     collapsed = false;
   }
-  post_collapse(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2,collapsed);
+  post_collapse(V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2,collapsed);
   if(collapsed)
   {
-    // Erase the two, other collapsed edges
-    Q.erase(Qit[e1]);
-    Qit[e1] = Q.end();
-    Q.erase(Qit[e2]);
-    Qit[e2] = Q.end();
+    // Erase the two, other collapsed edges by marking their timestamps as -1
+    EQ(e1) = -1;
+    EQ(e2) = -1;
+    // TODO: visits edges multiple times, ~150% more updates than should
+    //
     // update local neighbors
     // loop over original face neighbors
-    for(auto n : N)
+    //
+    // Can't use previous computed Nse and Nde because those refer to EMAP
+    // before it was changed...
+    std::vector<int> Nf;
+    Nf.reserve( Nsf.size() + Ndf.size() ); // preallocate memory
+    Nf.insert( Nf.end(), Nsf.begin(), Nsf.end() );
+    Nf.insert( Nf.end(), Ndf.begin(), Ndf.end() );
+    // https://stackoverflow.com/a/1041939/148668
+    std::sort( Nf.begin(), Nf.end() );
+    Nf.erase( std::unique( Nf.begin(), Nf.end() ), Nf.end() );
+    // Collect all edges that must be updated
+    std::vector<int> Ne;
+    Ne.reserve(3*Nf.size());
+    for(auto & n : Nf)
     {
       if(F(n,0) != IGL_COLLAPSE_EDGE_NULL ||
           F(n,1) != IGL_COLLAPSE_EDGE_NULL ||
@@ -371,24 +339,33 @@ IGL_INLINE bool igl::collapse_edge(
         {
           // get edge id
           const int ei = EMAP(v*F.rows()+n);
-          // erase old entry
-          Q.erase(Qit[ei]);
-          // compute cost and potential placement
-          double cost;
-          RowVectorXd place;
-          cost_and_placement(ei,V,F,E,EMAP,EF,EI,cost,place);
-          // Replace in queue
-          Qit[ei] = Q.insert(std::pair<double,int>(cost,ei)).first;
-          C.row(ei) = place;
+          Ne.push_back(ei);
         }
       }
     }
+    // Only process edge once
+    std::sort( Ne.begin(), Ne.end() );
+    Ne.erase( std::unique( Ne.begin(), Ne.end() ), Ne.end() );
+    for(auto & ei : Ne)
+    {
+       // compute cost and potential placement
+       double cost;
+       RowVectorXd place;
+       cost_and_placement(ei,V,F,E,EMAP,EF,EI,cost,place);
+       // Increment timestamp
+       EQ(ei)++;
+       // Replace in queue
+       Q.emplace(cost,ei,EQ(ei));
+       C.row(ei) = place;
+    }
   }else
   {
     // reinsert with infinite weight (the provided cost function must **not**
     // have given this un-collapsable edge inf cost already)
-    p.first = std::numeric_limits<double>::infinity();
-    Qit[e] = Q.insert(p).first;
+    // Increment timestamp
+    EQ(e)++;
+    // Replace in queue
+    Q.emplace(std::numeric_limits<double>::infinity(),e,EQ(e));
   }
   return collapsed;
 }

include/igl/collapse_edge.h

@@ -8,6 +8,8 @@
 #ifndef IGL_COLLAPSE_EDGE_H
 #define IGL_COLLAPSE_EDGE_H
 #include "igl_inline.h"
+#include "min_heap.h"
+#include "decimate.h" // for decimate_*_func types
 #include <Eigen/Core>
 #include <vector>
 #include <set>
@@ -55,6 +57,24 @@ namespace igl
     int & e2,
     int & f1,
     int & f2);
+  // Inputs:
+  IGL_INLINE bool collapse_edge(
+    const int e,
+    const Eigen::RowVectorXd & p,
+    /*const*/ std::vector<int> & Nsv,
+    const std::vector<int> & Nsf,
+    /*const*/ std::vector<int> & Ndv,
+    const std::vector<int> & Ndf,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI,
+    int & e1,
+    int & e2,
+    int & f1,
+    int & f2);
   IGL_INLINE bool collapse_edge(
     const int e,
     const Eigen::RowVectorXd & p,
@@ -73,28 +93,19 @@ namespace igl
   //     **If the edges is collapsed** then this function will be called on all
   //     edges of all faces previously incident on the endpoints of the
   //     collapsed edge.
-  //   Q  queue containing pairs of costs and edge indices
-  //   Qit  list of iterators so that Qit[e] --> iterator of edge e in Q
+  //   Q  queue containing pairs of costs and edge indices and insertion "time"
+  //   EQ  #E list of "time" of last time pushed into Q
   //   C  #E by dim list of stored placements
   IGL_INLINE bool collapse_edge(
-    const std::function<void(
-      const int,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement,
+    const decimate_cost_and_placement_func & cost_and_placement,
     Eigen::MatrixXd & V,
     Eigen::MatrixXi & F,
     Eigen::MatrixXi & E,
     Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI,
-    std::set<std::pair<double,int> > & Q,
-    std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    igl::min_heap< std::tuple<double,int,int> > & Q,
+    Eigen::VectorXi & EQ,
     Eigen::MatrixXd & C);
   // Inputs:
   //   pre_collapse  callback called with index of edge whose collapse is about
@@ -105,103 +116,31 @@ namespace igl
   //   post_collapse  callback called with index of edge whose collapse was
   //     just attempted and a flag revealing whether this was successful.
   IGL_INLINE bool collapse_edge(
-    const std::function<void(
-      const int,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement,
-    const std::function<bool(
-      const Eigen::MatrixXd &                                         ,/*V*/
-      const Eigen::MatrixXi &                                         ,/*F*/
-      const Eigen::MatrixXi &                                         ,/*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,/*EF*/
-      const Eigen::MatrixXi &                                         ,/*EI*/
-      const std::set<std::pair<double,int> > &                        ,/*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-      const Eigen::MatrixXd &                                         ,/*C*/
-      const int                                                        /*e*/
-      )> & pre_collapse,
-    const std::function<void(
-      const Eigen::MatrixXd &                                         ,   /*V*/
-      const Eigen::MatrixXi &                                         ,   /*F*/
-      const Eigen::MatrixXi &                                         ,   /*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,  /*EF*/
-      const Eigen::MatrixXi &                                         ,  /*EI*/
-      const std::set<std::pair<double,int> > &                        ,   /*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-      const Eigen::MatrixXd &                                         ,   /*C*/
-      const int                                                       ,   /*e*/
-      const int                                                       ,  /*e1*/
-      const int                                                       ,  /*e2*/
-      const int                                                       ,  /*f1*/
-      const int                                                       ,  /*f2*/
-      const bool                                                  /*collapsed*/
-      )> & post_collapse,
+    const decimate_cost_and_placement_func & cost_and_placement,
+    const decimate_pre_collapse_func       & pre_collapse,
+    const decimate_post_collapse_func      & post_collapse,
     Eigen::MatrixXd & V,
     Eigen::MatrixXi & F,
     Eigen::MatrixXi & E,
     Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI,
-    std::set<std::pair<double,int> > & Q,
-    std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    igl::min_heap< std::tuple<double,int,int> > & Q,
+    Eigen::VectorXi & EQ,
     Eigen::MatrixXd & C);
 
   IGL_INLINE bool collapse_edge(
-    const std::function<void(
-      const int,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement,
-    const std::function<bool(
-      const Eigen::MatrixXd &                                         ,/*V*/
-      const Eigen::MatrixXi &                                         ,/*F*/
-      const Eigen::MatrixXi &                                         ,/*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,/*EF*/
-      const Eigen::MatrixXi &                                         ,/*EI*/
-      const std::set<std::pair<double,int> > &                        ,/*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-      const Eigen::MatrixXd &                                         ,/*C*/
-      const int                                                        /*e*/
-      )> & pre_collapse,
-    const std::function<void(
-      const Eigen::MatrixXd &                                         ,   /*V*/
-      const Eigen::MatrixXi &                                         ,   /*F*/
-      const Eigen::MatrixXi &                                         ,   /*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,  /*EF*/
-      const Eigen::MatrixXi &                                         ,  /*EI*/
-      const std::set<std::pair<double,int> > &                        ,   /*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-      const Eigen::MatrixXd &                                         ,   /*C*/
-      const int                                                       ,   /*e*/
-      const int                                                       ,  /*e1*/
-      const int                                                       ,  /*e2*/
-      const int                                                       ,  /*f1*/
-      const int                                                       ,  /*f2*/
-      const bool                                                  /*collapsed*/
-      )> & post_collapse,
+    const decimate_cost_and_placement_func & cost_and_placement,
+    const decimate_pre_collapse_func       & pre_collapse,
+    const decimate_post_collapse_func      & post_collapse,
     Eigen::MatrixXd & V,
     Eigen::MatrixXi & F,
     Eigen::MatrixXi & E,
     Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI,
-    std::set<std::pair<double,int> > & Q,
-    std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    igl::min_heap< std::tuple<double,int,int> > & Q,
+    Eigen::VectorXi & EQ,
     Eigen::MatrixXd & C,
     int & e,
     int & e1,

include/igl/decimate.cpp

@@ -8,11 +8,13 @@
 #include "decimate.h"
 #include "collapse_edge.h"
 #include "edge_flaps.h"
+#include "always_try_never_care.h"
 #include "is_edge_manifold.h"
 #include "remove_unreferenced.h"
 #include "slice_mask.h"
 #include "slice.h"
 #include "connect_boundary_to_infinity.h"
+#include "parallel_for.h"
 #include "max_faces_stopping_condition.h"
 #include "shortest_edge_and_midpoint.h"
 
@@ -34,18 +36,34 @@ IGL_INLINE bool igl::decimate(
   DerivedV VO;
   DerivedF FO;
   igl::connect_boundary_to_infinity(V,F,VO,FO);
+  Eigen::VectorXi EMAP;
+  Eigen::MatrixXi E,EF,EI;
+  edge_flaps(FO,E,EMAP,EF,EI);
   // decimate will not work correctly on non-edge-manifold meshes. By extension
   // this includes meshes with non-manifold vertices on the boundary since these
   // will create a non-manifold edge when connected to infinity.
-  if(!is_edge_manifold(FO))
   {
-    return false;
+    Eigen::Array<bool,Eigen::Dynamic,Eigen::Dynamic> BF;
+    Eigen::Array<bool,Eigen::Dynamic,1> BE;
+    if(!is_edge_manifold(FO,E.rows(),EMAP,BF,BE))
+    {
+      return false;
+    }
   }
+  decimate_pre_collapse_func always_try;
+  decimate_post_collapse_func never_care;
+  always_try_never_care(always_try,never_care);
   bool ret = decimate(
     VO,
     FO,
     shortest_edge_and_midpoint,
     max_faces_stopping_condition(m,orig_m,max_m),
+    always_try,
+    never_care,
+    E,
+    EMAP,
+    EF,
+    EI,
     U,
     G,
     J,
@@ -82,35 +100,9 @@ IGL_INLINE bool igl::decimate(
   Eigen::VectorXi & I
   )
 {
-  const auto always_try = [](
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi &                                         ,/*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
-    const int                                                        /*e*/
-    ) -> bool { return true;};
-  const auto never_care = [](
-    const Eigen::MatrixXd &                                         ,   /*V*/
-    const Eigen::MatrixXi &                                         ,   /*F*/
-    const Eigen::MatrixXi &                                         ,   /*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,  /*EF*/
-    const Eigen::MatrixXi &                                         ,  /*EI*/
-    const std::set<std::pair<double,int> > &                        ,   /*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-    const Eigen::MatrixXd &                                         ,   /*C*/
-    const int                                                       ,   /*e*/
-    const int                                                       ,  /*e1*/
-    const int                                                       ,  /*e2*/
-    const int                                                       ,  /*f1*/
-    const int                                                       ,  /*f2*/
-    const bool                                                  /*collapsed*/
-    )-> void { };
+  decimate_pre_collapse_func always_try;
+  decimate_post_collapse_func never_care;
+  always_try_never_care(always_try,never_care);
   return igl::decimate(
     OV,OF,cost_and_placement,stopping_condition,always_try,never_care,U,G,J,I);
 }
@@ -128,10 +120,8 @@ IGL_INLINE bool igl::decimate(
   Eigen::VectorXi & I
   )
 {
-  using namespace Eigen;
-  using namespace std;
-  VectorXi EMAP;
-  MatrixXi E,EF,EI;
+  Eigen::VectorXi EMAP;
+  Eigen::MatrixXi E,EF,EI;
   edge_flaps(OF,E,EMAP,EF,EI);
   return igl::decimate(
     OV,OF,
@@ -157,31 +147,52 @@ IGL_INLINE bool igl::decimate(
   Eigen::VectorXi & I
   )
 {
-
   // Decimate 1
   using namespace Eigen;
   using namespace std;
   // Working copies
   Eigen::MatrixXd V = OV;
   Eigen::MatrixXi F = OF;
-  Eigen::MatrixXi E = OE;
-  Eigen::VectorXi EMAP = OEMAP;
-  Eigen::MatrixXi EF = OEF;
-  Eigen::MatrixXi EI = OEI;
-  typedef std::set<std::pair<double,int> > PriorityQueue;
-  PriorityQueue Q;
-  std::vector<PriorityQueue::iterator > Qit;
-  Qit.resize(E.rows());
+  VectorXi EMAP;
+  MatrixXi E,EF,EI;
+  edge_flaps(F,E,EMAP,EF,EI);
+  {
+    Eigen::Array<bool,Eigen::Dynamic,Eigen::Dynamic> BF;
+    Eigen::Array<bool,Eigen::Dynamic,1> BE;
+    if(!is_edge_manifold(F,E.rows(),EMAP,BF,BE))
+    {
+      return false;
+    }
+  }
+
+  igl::min_heap<std::tuple<double,int,int> > Q;
+  // Could reserve with https://stackoverflow.com/a/29236236/148668
+  Eigen::VectorXi EQ = Eigen::VectorXi::Zero(E.rows());
   // If an edge were collapsed, we'd collapse it to these points:
   MatrixXd C(E.rows(),V.cols());
-  for(int e = 0;e<E.rows();e++)
+  // Pushing into a vector then using constructor was slower. Maybe using
+  // std::move + make_heap would squeeze out something?
+  
+  // Separating the cost/placement evaluation from the Q filling is a
+  // performance hit for serial but faster if we can parallelize the
+  // cost/placement.
   {
-    double cost = e;
-    RowVectorXd p(1,3);
-    cost_and_placement(e,V,F,E,EMAP,EF,EI,cost,p);
-    C.row(e) = p;
-    Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
+    Eigen::VectorXd costs(E.rows());
+    igl::parallel_for(E.rows(),[&](const int e)
+    {
+      double cost = e;
+      RowVectorXd p(1,3);
+      cost_and_placement(e,V,F,E,EMAP,EF,EI,cost,p);
+      C.row(e) = p;
+      costs(e) = cost;
+    },10000);
+    for(int e = 0;e<E.rows();e++)
+    {
+      Q.emplace(costs(e),e,0);
+    }
   }
+
+
   int prev_e = -1;
   bool clean_finish = false;
 
@@ -191,17 +202,17 @@ IGL_INLINE bool igl::decimate(
     {
       break;
     }
-    if(Q.begin()->first == std::numeric_limits<double>::infinity())
+    if(std::get<0>(Q.top()) == std::numeric_limits<double>::infinity())
     {
       // min cost edge is infinite cost
       break;
     }
     int e,e1,e2,f1,f2;
     if(collapse_edge(
-       cost_and_placement, pre_collapse, post_collapse,
-       V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2))
+      cost_and_placement, pre_collapse, post_collapse,
+      V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2))
     {
-      if(stopping_condition(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2))
+      if(stopping_condition(V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2))
       {
         clean_finish = true;
         break;
@@ -237,6 +248,6 @@ IGL_INLINE bool igl::decimate(
   F2.conservativeResize(m,F2.cols());
   J.conservativeResize(m);
   VectorXi _1;
-  remove_unreferenced(V,F2,U,G,_1,I);
+  igl::remove_unreferenced(V,F2,U,G,_1,I);
   return clean_finish;
 }

include/igl/decimate.h

@@ -8,6 +8,7 @@
 #ifndef IGL_DECIMATE_H
 #define IGL_DECIMATE_H
 #include "igl_inline.h"
+#include "min_heap.h"
 #include <Eigen/Core>
 #include <vector>
 #include <set>
@@ -16,63 +17,63 @@ namespace igl
   // Function handles used to customize the `igl::decimate` command.
   using decimate_cost_and_placement_func = 
     std::function<void(
-      const int              /*e*/,
-      const Eigen::MatrixXd &/*V*/,
-      const Eigen::MatrixXi &/*F*/,
-      const Eigen::MatrixXi &/*E*/,
-      const Eigen::VectorXi &/*EMAP*/,
-      const Eigen::MatrixXi &/*EF*/,
-      const Eigen::MatrixXi &/*EI*/,
-      double &               /*cost*/,
-      Eigen::RowVectorXd &   /*p*/
+      const int                                           ,/*e*/
+      const Eigen::MatrixXd &                             ,/*V*/
+      const Eigen::MatrixXi &                             ,/*F*/
+      const Eigen::MatrixXi &                             ,/*E*/
+      const Eigen::VectorXi &                             ,/*EMAP*/
+      const Eigen::MatrixXi &                             ,/*EF*/
+      const Eigen::MatrixXi &                             ,/*EI*/
+      double &                                            ,/*cost*/
+      Eigen::RowVectorXd &                                 /*p*/
       )>;
   using decimate_stopping_condition_func = 
     std::function<bool(
-      const Eigen::MatrixXd &                                         ,/*V*/
-      const Eigen::MatrixXi &                                         ,/*F*/
-      const Eigen::MatrixXi &                                         ,/*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,/*EF*/
-      const Eigen::MatrixXi &                                         ,/*EI*/
-      const std::set<std::pair<double,int> > &                        ,/*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-      const Eigen::MatrixXd &                                         ,/*C*/
-      const int                                                       ,/*e*/
-      const int                                                       ,/*e1*/
-      const int                                                       ,/*e2*/
-      const int                                                       ,/*f1*/
-      const int                                                        /*f2*/
+      const Eigen::MatrixXd &                             ,/*V*/
+      const Eigen::MatrixXi &                             ,/*F*/
+      const Eigen::MatrixXi &                             ,/*E*/
+      const Eigen::VectorXi &                             ,/*EMAP*/
+      const Eigen::MatrixXi &                             ,/*EF*/
+      const Eigen::MatrixXi &                             ,/*EI*/
+      const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+      const Eigen::VectorXi &                             ,/*EQ*/
+      const Eigen::MatrixXd &                             ,/*C*/
+      const int                                           ,/*e*/
+      const int                                           ,/*e1*/
+      const int                                           ,/*e2*/
+      const int                                           ,/*f1*/
+      const int                                            /*f2*/
       )>;
   using decimate_pre_collapse_func = 
     std::function<bool(
-      const Eigen::MatrixXd &                                         ,/*V*/
-      const Eigen::MatrixXi &                                         ,/*F*/
-      const Eigen::MatrixXi &                                         ,/*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,/*EF*/
-      const Eigen::MatrixXi &                                         ,/*EI*/
-      const std::set<std::pair<double,int> > &                        ,/*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-      const Eigen::MatrixXd &                                         ,/*C*/
-      const int                                                        /*e*/
+      const Eigen::MatrixXd &                             ,/*V*/
+      const Eigen::MatrixXi &                             ,/*F*/
+      const Eigen::MatrixXi &                             ,/*E*/
+      const Eigen::VectorXi &                             ,/*EMAP*/
+      const Eigen::MatrixXi &                             ,/*EF*/
+      const Eigen::MatrixXi &                             ,/*EI*/
+      const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+      const Eigen::VectorXi &                             ,/*EQ*/
+      const Eigen::MatrixXd &                             ,/*C*/
+      const int                                            /*e*/
       )>;
   using decimate_post_collapse_func = 
     std::function<void(
-      const Eigen::MatrixXd &                                         ,   /*V*/
-      const Eigen::MatrixXi &                                         ,   /*F*/
-      const Eigen::MatrixXi &                                         ,   /*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,  /*EF*/
-      const Eigen::MatrixXi &                                         ,  /*EI*/
-      const std::set<std::pair<double,int> > &                        ,   /*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-      const Eigen::MatrixXd &                                         ,   /*C*/
-      const int                                                       ,   /*e*/
-      const int                                                       ,  /*e1*/
-      const int                                                       ,  /*e2*/
-      const int                                                       ,  /*f1*/
-      const int                                                       ,  /*f2*/
-      const bool                                                  /*collapsed*/
+      const Eigen::MatrixXd &                             ,/*V*/
+      const Eigen::MatrixXi &                             ,/*F*/
+      const Eigen::MatrixXi &                             ,/*E*/
+      const Eigen::VectorXi &                             ,/*EMAP*/
+      const Eigen::MatrixXi &                             ,/*EF*/
+      const Eigen::MatrixXi &                             ,/*EI*/
+      const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+      const Eigen::VectorXi &                             ,/*EQ*/
+      const Eigen::MatrixXd &                             ,/*C*/
+      const int                                           ,/*e*/
+      const int                                           ,/*e1*/
+      const int                                           ,/*e2*/
+      const int                                           ,/*f1*/
+      const int                                           ,/*f2*/
+      const bool                                           /*collapsed*/
       )>;
   // Assumes (V,F) is a manifold mesh (possibly with boundary) Collapses edges
   // until desired number of faces is achieved. This uses default edge cost and
@@ -135,7 +136,6 @@ namespace igl
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J,
     Eigen::VectorXi & I);
-
   // Inputs:
   //   pre_collapse  callback called with index of edge whose collapse is about
   //     to be attempted (see collapse_edge)
@@ -153,7 +153,6 @@ namespace igl
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J,
     Eigen::VectorXi & I);
-
   // Inputs:
   //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
   //     unique edge in E
@@ -161,7 +160,6 @@ namespace igl
   //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
   //     e=(j->i)
   //   EI  #E by 2 list of edge flap corners (see above).
-
   IGL_INLINE bool decimate(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,

include/igl/edge_collapse_is_valid.cpp

@@ -84,3 +84,41 @@ IGL_INLINE bool igl::edge_collapse_is_valid(
   }
   return true;
 }
+
+IGL_INLINE bool igl::edge_collapse_is_valid(
+  std::vector<int> & Nsv,
+  std::vector<int> & Ndv)
+{
+  // Do we really need to check if edge is IGL_COLLAPSE_EDGE_NULL ?
+
+  if(Nsv.size()<2 || Ndv.size()<2)
+  {
+    // Bogus data
+    assert(false);
+    return false;
+  }
+  // determine if the first two vertices are the same before reordering.
+  // If they are and there are 3 each, then (I claim) this is an edge on a
+  // single tet.
+  const bool first_two_same = (Nsv[0] == Ndv[0]) && (Nsv[1] == Ndv[1]);
+  if(Nsv.size() == 3 && Ndv.size() == 3 && first_two_same)
+  {
+    // single tet
+    return false;
+  }
+  // https://stackoverflow.com/a/19483741/148668
+  std::sort(Nsv.begin(), Nsv.end());
+  std::sort(Ndv.begin(), Ndv.end());
+  std::vector<int> Nint;
+  std::set_intersection(
+    Nsv.begin(), Nsv.end(), Ndv.begin(), Ndv.end(), std::back_inserter(Nint));
+  // check if edge collapse is valid: intersection of vertex neighbors of s and
+  // d should be exactly 2+(s,d) = 4
+  // http://stackoverflow.com/a/27049418/148668
+  if(Nint.size() != 2)
+  {
+    return false;
+  }
+  
+  return true;
+}

include/igl/edge_collapse_is_valid.h

@@ -9,6 +9,7 @@
 #define IGL_EDGE_COLLAPSE_IS_VALID_H
 #include "igl_inline.h"
 #include <Eigen/Core>
+#include <vector>
 namespace igl
 {
   // Assumes (V,F) is a closed manifold mesh (except for previouslly collapsed
@@ -37,6 +38,20 @@ namespace igl
     const Eigen::VectorXi & EMAP,
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI);
+  // Inputs:
+  //   Nsv  #Nsv list of "next" vertices circulating around starting vertex of
+  //     edge
+  //   Ndv  #Ndv list of "next" vertices circulating around destination vertex of
+  //     edge
+  // Outputs:
+  //   Nsv  (side-effect: sorted by value)
+  //   Ndv  (side-effect: sorted by value)
+  // Returns true iff edge collapse is valid
+  //
+  // See also: circulation
+  IGL_INLINE bool edge_collapse_is_valid(
+    /*const*/ std::vector<int> & Nsv,
+    /*const*/ std::vector<int> & Ndv);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "edge_collapse_is_valid.cpp"

include/igl/edge_flaps.cpp

@@ -51,8 +51,7 @@ IGL_INLINE void igl::edge_flaps(
   Eigen::MatrixXi & EI)
 {
   Eigen::MatrixXi allE;
-  std::vector<std::vector<int> > uE2E;
-  igl::unique_edge_map(F,allE,uE,EMAP,uE2E);
+  igl::unique_edge_map(F,allE,uE,EMAP);
   // Const-ify to call overload
   const auto & cuE = uE;
   const auto & cEMAP = EMAP;

include/igl/edge_flaps.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2015 Alec Jacobson <[email protected]>
+// 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 

include/igl/infinite_cost_stopping_condition.cpp

@@ -8,31 +8,8 @@
 #include "infinite_cost_stopping_condition.h"
 
 IGL_INLINE void igl::infinite_cost_stopping_condition(
-  const std::function<void(
-    const int,
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    double &,
-    Eigen::RowVectorXd &)> & cost_and_placement,
-  std::function<bool(
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
-    const Eigen::MatrixXd &,
-    const int,
-    const int,
-    const int,
-    const int,
-    const int)> & stopping_condition)
+  const decimate_cost_and_placement_func & cost_and_placement,
+  decimate_stopping_condition_func & stopping_condition)
 {
   stopping_condition = 
     [&cost_and_placement]
@@ -43,9 +20,9 @@ IGL_INLINE void igl::infinite_cost_stopping_condition(
     const Eigen::VectorXi & EMAP,
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI,
-    const std::set<std::pair<double,int> > & Q,
-    const std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
-    const Eigen::MatrixXd & C,
+    const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+    const Eigen::VectorXi &                             ,/*EQ*/
+    const Eigen::MatrixXd & /*C*/,
     const int e,
     const int /*e1*/,
     const int /*e2*/,
@@ -59,49 +36,11 @@ IGL_INLINE void igl::infinite_cost_stopping_condition(
     };
 }
 
-IGL_INLINE 
-  std::function<bool(
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
-    const Eigen::MatrixXd &,
-    const int,
-    const int,
-    const int,
-    const int,
-    const int)> 
+IGL_INLINE igl::decimate_stopping_condition_func
   igl::infinite_cost_stopping_condition(
-    const std::function<void(
-      const int,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement)
+  const decimate_cost_and_placement_func & cost_and_placement)
 {
-  std::function<bool(
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
-    const Eigen::MatrixXd &,
-    const int,
-    const int,
-    const int,
-    const int,
-    const int)> stopping_condition;
+  decimate_stopping_condition_func stopping_condition;
   infinite_cost_stopping_condition(cost_and_placement,stopping_condition);
   return stopping_condition;
 }

include/igl/infinite_cost_stopping_condition.h

@@ -8,6 +8,7 @@
 #ifndef IGL_INFINITE_COST_STOPPING_CONDITION_H
 #define IGL_INFINITE_COST_STOPPING_CONDITION_H
 #include "igl_inline.h"
+#include "decimate.h" // decimate_*_func type definitions
 #include <Eigen/Core>
 #include <vector>
 #include <set>
@@ -23,58 +24,11 @@ namespace igl
   //   stopping_condition
   //
   IGL_INLINE void infinite_cost_stopping_condition(
-    const std::function<void(
-      const int,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement,
-    std::function<bool(
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const std::set<std::pair<double,int> > &,
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,
-      const Eigen::MatrixXd &,
-      const int,
-      const int,
-      const int,
-      const int,
-      const int)> & stopping_condition);
-  IGL_INLINE 
-    std::function<bool(
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const std::set<std::pair<double,int> > &,
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,
-      const Eigen::MatrixXd &,
-      const int,
-      const int,
-      const int,
-      const int,
-      const int)> 
+    const decimate_cost_and_placement_func & cost_and_placement,
+    decimate_stopping_condition_func & stopping_condition);
+  IGL_INLINE decimate_stopping_condition_func
     infinite_cost_stopping_condition(
-      const std::function<void(
-        const int,
-        const Eigen::MatrixXd &,
-        const Eigen::MatrixXi &,
-        const Eigen::MatrixXi &,
-        const Eigen::VectorXi &,
-        const Eigen::MatrixXi &,
-        const Eigen::MatrixXi &,
-        double &,
-        Eigen::RowVectorXd &)> & cost_and_placement);
+    const decimate_cost_and_placement_func & cost_and_placement);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/is_edge_manifold.cpp

@@ -8,42 +8,32 @@
 #include "is_edge_manifold.h"
 #include "oriented_facets.h"
 #include "unique_simplices.h"
+#include "unique_edge_map.h"
 
 #include <algorithm>
 #include <vector>
 
 template <
   typename DerivedF,
-  typename DerivedBF,
-  typename DerivedE,
   typename DerivedEMAP,
+  typename DerivedBF,
   typename DerivedBE>
 IGL_INLINE bool igl::is_edge_manifold(
   const Eigen::MatrixBase<DerivedF>& F,
+  const typename DerivedF::Index ne,
+  const Eigen::MatrixBase<DerivedEMAP>& EMAP,
   Eigen::PlainObjectBase<DerivedBF>& BF,
-  Eigen::PlainObjectBase<DerivedE>& E,
-  Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
   Eigen::PlainObjectBase<DerivedBE>& BE)
 {
-  using namespace Eigen;
   typedef typename DerivedF::Index Index;
-  typedef Matrix<typename DerivedF::Scalar,Dynamic,1> VectorXF;
-  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixXF2;
-  MatrixXF2 allE;
-  oriented_facets(F,allE);
-  // Find unique undirected edges and mapping
-  VectorXF _;
-  unique_simplices(allE,E,_,EMAP);
-  // could use "face_occurrences.h", but that implementation uses
-  // vector<vector>>
-  std::vector<typename DerivedF::Index> count(E.rows(),0);
+  std::vector<Index> count(ne,0);
   for(Index e = 0;e<EMAP.rows();e++)
   {
     count[EMAP[e]]++;
   }
   const Index m = F.rows();
   BF.resize(m,3);
-  BE.resize(E.rows(),1);
+  BE.resize(ne,1);
   bool all = true;
   for(Index e = 0;e<EMAP.rows();e++)
   {
@@ -54,6 +44,26 @@ IGL_INLINE bool igl::is_edge_manifold(
   return all;
 }
 
+template <
+  typename DerivedF,
+  typename DerivedBF,
+  typename DerivedE,
+  typename DerivedEMAP,
+  typename DerivedBE>
+IGL_INLINE bool igl::is_edge_manifold(
+  const Eigen::MatrixBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedBF>& BF,
+  Eigen::PlainObjectBase<DerivedE>& E,
+  Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  Eigen::PlainObjectBase<DerivedBE>& BE)
+{
+  using namespace Eigen;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixXF2;
+  MatrixXF2 allE;
+  unique_edge_map(F,allE,E,EMAP);
+  return is_edge_manifold(F,E.rows(),EMAP,BF,BE);
+}
+
 template <typename DerivedF>
 IGL_INLINE bool igl::is_edge_manifold(
   const Eigen::MatrixBase<DerivedF>& F)

include/igl/is_edge_manifold.h

@@ -18,12 +18,20 @@ namespace igl
   //
   // Inputs:
   //   F  #F by 3 list of triangle indices
-  // Outputs:
-  //   BF  #F by 3 list of flags revealing if edge opposite corresponding vertex
-  //   is non-manifold.
   // Returns true iff all edges are manifold
   //
   // See also: is_vertex_manifold
+  template <typename DerivedF>
+  IGL_INLINE bool is_edge_manifold(
+    const Eigen::MatrixBase<DerivedF>& F);
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   BF  #F by 3 list of flags revealing if edge opposite corresponding vertex
+  //   is non-manifold.
+  //   E  #E by 2 list of unique edges
+  //   EMAP  3*#F list of indices of opposite edges in "E"
+  //   BE  #E list of flages whether edge is non-manifold
   template <
     typename DerivedF, 
     typename DerivedBF,
@@ -36,9 +44,25 @@ namespace igl
     Eigen::PlainObjectBase<DerivedE>& E,
     Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
     Eigen::PlainObjectBase<DerivedBE>& BE);
-  template <typename DerivedF>
+  // Inputs:
+  //   F  #F by 3 list of triangle indices
+  //   ne  number of edges (#E)
+  //   EMAP  3*#F list of indices of opposite edges in "E"
+  // Outputs:
+  //   BF  #F by 3 list of flags revealing if edge opposite corresponding vertex
+  //     is non-manifold.
+  //   BE  ne list of flages whether edge is non-manifold
+  template <
+    typename DerivedF,
+    typename DerivedEMAP,
+    typename DerivedBF,
+    typename DerivedBE>
   IGL_INLINE bool is_edge_manifold(
-    const Eigen::MatrixBase<DerivedF>& F);
+    const Eigen::MatrixBase<DerivedF>& F,
+    const typename DerivedF::Index ne,
+    const Eigen::MatrixBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedBF>& BF,
+    Eigen::PlainObjectBase<DerivedBE>& BE);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/max_faces_stopping_condition.cpp

@@ -11,21 +11,7 @@ IGL_INLINE void igl::max_faces_stopping_condition(
   int & m,
   const int orig_m,
   const int max_m,
-  std::function<bool(
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
-    const Eigen::MatrixXd &,
-    const int,
-    const int,
-    const int,
-    const int,
-    const int)> & stopping_condition)
+  decimate_stopping_condition_func & stopping_condition)
 {
   stopping_condition = 
     [orig_m,max_m,&m](
@@ -35,8 +21,8 @@ IGL_INLINE void igl::max_faces_stopping_condition(
     const Eigen::VectorXi &,
     const Eigen::MatrixXi &,
     const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
+    const igl::min_heap< std::tuple<double,int,int> > & ,
+    const Eigen::VectorXi &                             ,
     const Eigen::MatrixXd &,
     const int,
     const int,
@@ -51,43 +37,13 @@ IGL_INLINE void igl::max_faces_stopping_condition(
     };
 }
 
-IGL_INLINE 
-  std::function<bool(
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
-    const Eigen::MatrixXd &,
-    const int,
-    const int,
-    const int,
-    const int,
-    const int)> 
-  igl::max_faces_stopping_condition(
-    int & m,
-    const int orig_m,
-    const int max_m)
+IGL_INLINE igl::decimate_stopping_condition_func
+igl::max_faces_stopping_condition(
+  int & m,
+  const int orig_m,
+  const int max_m)
 {
-  std::function<bool(
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const std::set<std::pair<double,int> > &,
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,
-    const Eigen::MatrixXd &,
-    const int,
-    const int,
-    const int,
-    const int,
-    const int)> stopping_condition;
-  max_faces_stopping_condition(
-      m,orig_m,max_m,stopping_condition);
+  decimate_stopping_condition_func stopping_condition;
+  max_faces_stopping_condition(m,orig_m,max_m,stopping_condition);
   return stopping_condition;
 }

include/igl/max_faces_stopping_condition.h

@@ -8,6 +8,7 @@
 #ifndef IGL_MAX_FACES_STOPPING_CONDITION_H
 #define IGL_MAX_FACES_STOPPING_CONDITION_H
 #include "igl_inline.h"
+#include "decimate.h" // for decimate_*_func types
 #include <Eigen/Core>
 #include <vector>
 #include <set>
@@ -31,37 +32,8 @@ namespace igl
     int & m,
     const int orig_m,
     const int max_m,
-    std::function<bool(
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const std::set<std::pair<double,int> > &,
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,
-      const Eigen::MatrixXd &,
-      const int,
-      const int,
-      const int,
-      const int,
-      const int)> & stopping_condition);
-  IGL_INLINE 
-    std::function<bool(
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const std::set<std::pair<double,int> > &,
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,
-      const Eigen::MatrixXd &,
-      const int,
-      const int,
-      const int,
-      const int,
-      const int)> 
+    decimate_stopping_condition_func & stopping_condition);
+  IGL_INLINE decimate_stopping_condition_func
     max_faces_stopping_condition(
       int & m,
       const int orign_m,

include/igl/min_heap.h

@@ -0,0 +1,20 @@
+// 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_MIN_HEAP_H
+#define IGL_MIN_HEAP_H
+#include <queue>
+#include <vector>
+#include <functional>
+namespace igl
+{
+  // Templated min heap (reverses sort order of std::priority_queue)
+  template<class T> using min_heap = 
+    std::priority_queue< T, std::vector<T >, std::greater<T > >;
+}
+#endif 
+

include/igl/qslim.cpp

@@ -58,45 +58,9 @@ IGL_INLINE bool igl::qslim(
   int v1 = -1;
   int v2 = -1;
   // Callbacks for computing and updating metric
-  std::function<void(
-    const int e,
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    double &,
-    Eigen::RowVectorXd &)> cost_and_placement;
-  std::function<bool(
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi &                                         ,/*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
-    const int                                                        /*e*/
-    )> pre_collapse;
-  std::function<void(
-    const Eigen::MatrixXd &                                         ,   /*V*/
-    const Eigen::MatrixXi &                                         ,   /*F*/
-    const Eigen::MatrixXi &                                         ,   /*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,  /*EF*/
-    const Eigen::MatrixXi &                                         ,  /*EI*/
-    const std::set<std::pair<double,int> > &                        ,   /*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-    const Eigen::MatrixXd &                                         ,   /*C*/
-    const int                                                       ,   /*e*/
-    const int                                                       ,  /*e1*/
-    const int                                                       ,  /*e2*/
-    const int                                                       ,  /*f1*/
-    const int                                                       ,  /*f2*/
-    const bool                                                  /*collapsed*/
-    )> post_collapse;
+  decimate_cost_and_placement_func cost_and_placement;
+  decimate_pre_collapse_func       pre_collapse;
+  decimate_post_collapse_func      post_collapse;
   qslim_optimal_collapse_edge_callbacks(
     E,quadrics,v1,v2, cost_and_placement, pre_collapse,post_collapse);
   // Call to greedy decimator

include/igl/qslim_optimal_collapse_edge_callbacks.cpp

@@ -15,45 +15,9 @@ IGL_INLINE void igl::qslim_optimal_collapse_edge_callbacks(
     quadrics,
   int & v1,
   int & v2,
-  std::function<void(
-    const int e,
-    const Eigen::MatrixXd &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::VectorXi &,
-    const Eigen::MatrixXi &,
-    const Eigen::MatrixXi &,
-    double &,
-    Eigen::RowVectorXd &)> & cost_and_placement,
-  std::function<bool(
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi &                                         ,/*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
-    const int                                                        /*e*/
-    )> & pre_collapse,
-  std::function<void(
-    const Eigen::MatrixXd &                                         ,   /*V*/
-    const Eigen::MatrixXi &                                         ,   /*F*/
-    const Eigen::MatrixXi &                                         ,   /*E*/
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,  /*EF*/
-    const Eigen::MatrixXi &                                         ,  /*EI*/
-    const std::set<std::pair<double,int> > &                        ,   /*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-    const Eigen::MatrixXd &                                         ,   /*C*/
-    const int                                                       ,   /*e*/
-    const int                                                       ,  /*e1*/
-    const int                                                       ,  /*e2*/
-    const int                                                       ,  /*f1*/
-    const int                                                       ,  /*f2*/
-    const bool                                                  /*collapsed*/
-    )> & post_collapse)
+  decimate_cost_and_placement_func & cost_and_placement,
+  decimate_pre_collapse_func       & pre_collapse,
+  decimate_post_collapse_func      & post_collapse)
 {
   typedef std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> Quadric;
   cost_and_placement = [&quadrics,&v1,&v2](
@@ -87,15 +51,15 @@ IGL_INLINE void igl::qslim_optimal_collapse_edge_callbacks(
   };
   // Remember endpoints
   pre_collapse = [&v1,&v2](
-    const Eigen::MatrixXd &                                         ,/*V*/
-    const Eigen::MatrixXi &                                         ,/*F*/
-    const Eigen::MatrixXi & E,
-    const Eigen::VectorXi &                                         ,/*EMAP*/
-    const Eigen::MatrixXi &                                         ,/*EF*/
-    const Eigen::MatrixXi &                                         ,/*EI*/
-    const std::set<std::pair<double,int> > &                        ,/*Q*/
-    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-    const Eigen::MatrixXd &                                         ,/*C*/
+    const Eigen::MatrixXd &                             ,/*V*/
+    const Eigen::MatrixXi &                             ,/*F*/
+    const Eigen::MatrixXi & E                           ,
+    const Eigen::VectorXi &                             ,/*EMAP*/
+    const Eigen::MatrixXi &                             ,/*EF*/
+    const Eigen::MatrixXi &                             ,/*EI*/
+    const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+    const Eigen::VectorXi &                             ,/*EQ*/
+    const Eigen::MatrixXd &                             ,/*C*/
     const int e)->bool
   {
     v1 = E(e,0);
@@ -104,21 +68,21 @@ IGL_INLINE void igl::qslim_optimal_collapse_edge_callbacks(
   };
   // update quadric
   post_collapse = [&v1,&v2,&quadrics](
-      const Eigen::MatrixXd &                                         ,   /*V*/
-      const Eigen::MatrixXi &                                         ,   /*F*/
-      const Eigen::MatrixXi &                                         ,   /*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,  /*EF*/
-      const Eigen::MatrixXi &                                         ,  /*EI*/
-      const std::set<std::pair<double,int> > &                        ,   /*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-      const Eigen::MatrixXd &                                         ,   /*C*/
-      const int                                                       ,   /*e*/
-      const int                                                       ,  /*e1*/
-      const int                                                       ,  /*e2*/
-      const int                                                       ,  /*f1*/
-      const int                                                       ,  /*f2*/
-      const bool                                                  collapsed
+      const Eigen::MatrixXd &                             ,   /*V*/
+      const Eigen::MatrixXi &                             ,   /*F*/
+      const Eigen::MatrixXi &                             ,   /*E*/
+      const Eigen::VectorXi &                             ,/*EMAP*/
+      const Eigen::MatrixXi &                             ,  /*EF*/
+      const Eigen::MatrixXi &                             ,  /*EI*/
+      const igl::min_heap< std::tuple<double,int,int> > & ,/*Q*/
+      const Eigen::VectorXi &                             ,/*EQ*/
+      const Eigen::MatrixXd &                             ,   /*C*/
+      const int                                           ,   /*e*/
+      const int                                           ,  /*e1*/
+      const int                                           ,  /*e2*/
+      const int                                           ,  /*f1*/
+      const int                                           ,  /*f2*/
+      const bool                                          collapsed
       )->void
   {
     if(collapsed)

include/igl/qslim_optimal_collapse_edge_callbacks.h

@@ -8,6 +8,7 @@
 #ifndef IGL_QSLIM_OPTIMAL_COLLAPSE_EDGE_CALLBACKS_H
 #define IGL_QSLIM_OPTIMAL_COLLAPSE_EDGE_CALLBACKS_H
 #include "igl_inline.h"
+#include "decimate.h" // decimate_*_func type definitions
 #include <Eigen/Core>
 #include <functional>
 #include <vector>
@@ -15,7 +16,6 @@
 #include <set>
 namespace igl
 {
-
   // Prepare callbacks for decimating edges using the qslim optimal placement
   // metric.
   //
@@ -35,45 +35,9 @@ namespace igl
       quadrics,
     int & v1,
     int & v2,
-    std::function<void(
-      const int e,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement,
-    std::function<bool(
-      const Eigen::MatrixXd &                                         ,/*V*/
-      const Eigen::MatrixXi &                                         ,/*F*/
-      const Eigen::MatrixXi &                                         ,/*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,/*EF*/
-      const Eigen::MatrixXi &                                         ,/*EI*/
-      const std::set<std::pair<double,int> > &                        ,/*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
-      const Eigen::MatrixXd &                                         ,/*C*/
-      const int                                                        /*e*/
-      )> & pre_collapse,
-    std::function<void(
-      const Eigen::MatrixXd &                                         ,   /*V*/
-      const Eigen::MatrixXi &                                         ,   /*F*/
-      const Eigen::MatrixXi &                                         ,   /*E*/
-      const Eigen::VectorXi &                                         ,/*EMAP*/
-      const Eigen::MatrixXi &                                         ,  /*EF*/
-      const Eigen::MatrixXi &                                         ,  /*EI*/
-      const std::set<std::pair<double,int> > &                        ,   /*Q*/
-      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
-      const Eigen::MatrixXd &                                         ,   /*C*/
-      const int                                                       ,   /*e*/
-      const int                                                       ,  /*e1*/
-      const int                                                       ,  /*e2*/
-      const int                                                       ,  /*f1*/
-      const int                                                       ,  /*f2*/
-      const bool                                                  /*collapsed*/
-      )> & post_collapse);
+    decimate_cost_and_placement_func & cost_and_placement,
+    decimate_pre_collapse_func & pre_collapse,
+    decimate_post_collapse_func & post_collapse);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "qslim_optimal_collapse_edge_callbacks.cpp"

include/igl/unique_edge_map.cpp

@@ -106,6 +106,14 @@ IGL_INLINE void igl::unique_edge_map(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long>(Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);

include/igl/unique_rows.cpp

@@ -75,6 +75,8 @@ IGL_INLINE void igl::unique_rows(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::unique_rows<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::unique_rows<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::unique_rows<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/unique_simplices.cpp

@@ -53,6 +53,8 @@ IGL_INLINE void igl::unique_simplices(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::unique_simplices<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::unique_simplices<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::unique_simplices<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

tutorial/703_Decimation/main.cpp

@@ -1,7 +1,9 @@
 #include <igl/circulation.h>
 #include <igl/collapse_edge.h>
 #include <igl/edge_flaps.h>
+#include <igl/decimate.h>
 #include <igl/shortest_edge_and_midpoint.h>
+#include <igl/parallel_for.h>
 #include <igl/read_triangle_mesh.h>
 #include <igl/opengl/glfw/Viewer.h>
 #include <Eigen/Core>
@@ -33,9 +35,8 @@ int main(int argc, char * argv[])
   // Prepare array-based edge data structures and priority queue
   VectorXi EMAP;
   MatrixXi E,EF,EI;
-  typedef std::set<std::pair<double,int> > PriorityQueue;
-  PriorityQueue Q;
-  std::vector<PriorityQueue::iterator > Qit;
+  igl::min_heap< std::tuple<double,int,int> > Q;
+  Eigen::VectorXi EQ;
   // If an edge were collapsed, we'd collapse it to these points:
   MatrixXd C;
   int num_collapsed;
@@ -46,19 +47,28 @@ int main(int argc, char * argv[])
     F = OF;
     V = OV;
     edge_flaps(F,E,EMAP,EF,EI);
-    Qit.resize(E.rows());
-
     C.resize(E.rows(),V.cols());
     VectorXd costs(E.rows());
-    Q.clear();
-    for(int e = 0;e<E.rows();e++)
+    // https://stackoverflow.com/questions/2852140/priority-queue-clear-method
+    // Q.clear();
+    Q = {};
+    EQ = Eigen::VectorXi::Zero(E.rows());
     {
-      double cost = e;
-      RowVectorXd p(1,3);
-      shortest_edge_and_midpoint(e,V,F,E,EMAP,EF,EI,cost,p);
-      C.row(e) = p;
-      Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
+      Eigen::VectorXd costs(E.rows());
+      igl::parallel_for(E.rows(),[&](const int e)
+      {
+        double cost = e;
+        RowVectorXd p(1,3);
+        shortest_edge_and_midpoint(e,V,F,E,EMAP,EF,EI,cost,p);
+        C.row(e) = p;
+        costs(e) = cost;
+      },10000);
+      for(int e = 0;e<E.rows();e++)
+      {
+        Q.emplace(costs(e),e,0);
+      }
     }
+
     num_collapsed = 0;
     viewer.data().clear();
     viewer.data().set_mesh(V,F);
@@ -75,8 +85,7 @@ int main(int argc, char * argv[])
       const int max_iter = std::ceil(0.01*Q.size());
       for(int j = 0;j<max_iter;j++)
       {
-        if(!collapse_edge(
-          shortest_edge_and_midpoint, V,F,E,EMAP,EF,EI,Q,Qit,C))
+        if(!collapse_edge(shortest_edge_and_midpoint,V,F,E,EMAP,EF,EI,Q,EQ,C))
         {
           break;
         }