pre and post edge collapse callbacks, qslim

Former-commit-id: a85dc8372805f76beba123bfe636b1dd51cc455f

include/igl/collapse_edge.cpp

@@ -165,6 +165,8 @@ IGL_INLINE bool igl::collapse_edge(
     const Eigen::MatrixXi &,
     double &,
     Eigen::RowVectorXd &)> & cost_and_placement,
+  const std::function<bool(const int )> & pre_collapse,
+  const std::function<void(const int , const bool)> & post_collapse,
   Eigen::MatrixXd & V,
   Eigen::MatrixXi & F,
   Eigen::MatrixXi & E,
@@ -177,7 +179,9 @@ IGL_INLINE bool igl::collapse_edge(
 {
   int e,e1,e2,f1,f2;
   return 
-    collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+    collapse_edge(
+      cost_and_placement,pre_collapse,post_collapse,
+      V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
 }
 
 
@@ -192,6 +196,8 @@ IGL_INLINE bool igl::collapse_edge(
     const Eigen::MatrixXi &,
     double &,
     Eigen::RowVectorXd &)> & cost_and_placement,
+  const std::function<bool(const int )> & pre_collapse,
+  const std::function<void(const int , const bool)> & post_collapse,
   Eigen::MatrixXd & V,
   Eigen::MatrixXi & F,
   Eigen::MatrixXi & E,
@@ -225,8 +231,16 @@ IGL_INLINE bool igl::collapse_edge(
   std::vector<int> N  = circulation(e, true,F,E,EMAP,EF,EI);
   std::vector<int> Nd = circulation(e,false,F,E,EMAP,EF,EI);
   N.insert(N.begin(),Nd.begin(),Nd.end());
-  const bool collapsed =
-    collapse_edge(e,C.row(e),V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+  bool collapsed = true;
+  if(pre_collapse(e))
+  {
+    collapsed = collapse_edge(e,C.row(e),V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+  }else
+  {
+    // Aborted by pre collapse callback
+    collapsed = false;
+  }
+  post_collapse(e,collapsed);
   if(collapsed)
   {
     // Erase the two, other collapsed edges

include/igl/collapse_edge.h

@@ -70,6 +70,16 @@ namespace igl
   //   cost_and_placement  function computing cost of collapsing an edge and 3d
   //     position where it should be placed:
   //     cost_and_placement(V,F,E,EMAP,EF,EI,cost,placement);
+  //     **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.
+  //   pre_collapse  callback called with index of edge whose collapse is about
+  //     to be attempted. This function should return whether to **proceed**
+  //     with the collapse: returning true means "yes, try to collapse",
+  //     returning false means "No, consider this edge 'uncollapsable', behave
+  //     as if collapse_edge(e) returned false.
+  //   post_collapse  callback called with index of edge whose collapse was
+  //     just attempted and a flag revealing whether this was successful.
   //   Q  queue containing pairs of costs and edge indices
   //   Qit  list of iterators so that Qit[e] --> iterator of edge e in Q
   //   C  #E by dim list of stored placements
@@ -84,6 +94,8 @@ namespace igl
       const Eigen::MatrixXi &,
       double &,
       Eigen::RowVectorXd &)> & cost_and_placement,
+    const std::function<bool(const int )> & pre_collapse,
+    const std::function<void(const int , const bool)> & post_collapse,
     Eigen::MatrixXd & V,
     Eigen::MatrixXi & F,
     Eigen::MatrixXi & E,
@@ -104,6 +116,8 @@ namespace igl
       const Eigen::MatrixXi &,
       double &,
       Eigen::RowVectorXd &)> & cost_and_placement,
+    const std::function<bool(const int )> & pre_collapse,
+    const std::function<void(const int , const bool)> & post_collapse,
     Eigen::MatrixXd & V,
     Eigen::MatrixXi & F,
     Eigen::MatrixXi & E,

include/igl/decimate.cpp

@@ -35,11 +35,15 @@ IGL_INLINE bool igl::decimate(
   DerivedV VO;
   DerivedF FO;
   igl::connect_boundary_to_infinity(V,F,VO,FO);
+  const auto & always_try = [](const int e)->bool{return true;};
+  const auto & never_care = [](const int e, const bool collapsed){};
   bool ret = decimate(
     VO,
     FO,
     shortest_edge_and_midpoint,
     max_faces_stopping_condition(m,orig_m,max_m),
+    always_try,
+    never_care,
     U,
     G,
     J,
@@ -93,6 +97,61 @@ IGL_INLINE bool igl::decimate(
       const int,
       const int,
       const int)> & stopping_condition,
+  const std::function<bool(const int )> & pre_collapse,
+  const std::function<void(const int , const bool)> & post_collapse,
+  Eigen::MatrixXd & U,
+  Eigen::MatrixXi & G,
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & I
+  )
+{
+  using namespace Eigen;
+  using namespace std;
+  VectorXi EMAP;
+  MatrixXi E,EF,EI;
+  edge_flaps(OF,E,EMAP,EF,EI);
+  return igl::decimate(
+    OV,OF,
+    cost_and_placement,stopping_condition,pre_collapse,post_collapse,
+    E,EMAP,EF,EI,
+    U,G,J,I);
+}
+
+
+IGL_INLINE bool igl::decimate(
+  const Eigen::MatrixXd & OV,
+  const Eigen::MatrixXi & OF,
+  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 &,
+      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 std::function<bool(const int )> & pre_collapse,
+  const std::function<void(const int , const bool)> & post_collapse,
+  const Eigen::MatrixXi & OE,
+  const Eigen::VectorXi & OEMAP,
+  const Eigen::MatrixXi & OEF,
+  const Eigen::MatrixXi & OEI,
   Eigen::MatrixXd & U,
   Eigen::MatrixXi & G,
   Eigen::VectorXi & J,
@@ -104,12 +163,14 @@ IGL_INLINE bool igl::decimate(
   // Working copies
   Eigen::MatrixXd V = OV;
   Eigen::MatrixXi F = OF;
-  VectorXi EMAP;
-  MatrixXi E,EF,EI;
+  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;
-  edge_flaps(F,E,EMAP,EF,EI);
   Qit.resize(E.rows());
   // If an edge were collapsed, we'd collapse it to these points:
   MatrixXd C(E.rows(),V.cols());
@@ -136,7 +197,9 @@ IGL_INLINE bool igl::decimate(
       break;
     }
     int e,e1,e2,f1,f2;
-    if(collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C,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))
     {
       if(stopping_condition(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2))
       {

include/igl/decimate.h

@@ -65,6 +65,11 @@ namespace igl
   //     collapsing edge e removing edges (e,e1,e2) and faces (f1,f2):
   //     bool should_stop =
   //       stopping_condition(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+  //   pre_collapse  callback called with index of edge whose collapse is about
+  //     to be attempted (see collapse_edge)
+  //   post_collapse  callback called with index of edge whose collapse was
+  //     just attempted and a flag revealing whether this was successful (see
+  //     collapse_edge)
   IGL_INLINE bool decimate(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
@@ -95,6 +100,55 @@ namespace igl
       const int                                                       ,/*f1*/
       const int                                                        /*f2*/
       )> & stopping_condition,
+    const std::function<bool(const int )> & pre_collapse,
+    const std::function<void(const int , const bool)> & post_collapse,
+    Eigen::MatrixXd & U,
+    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
+  //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  //     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,
+    const 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*/
+      )> & 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*/
+      const int                                                       ,/*e1*/
+      const int                                                       ,/*e2*/
+      const int                                                       ,/*f1*/
+      const int                                                        /*f2*/
+      )> & stopping_condition,
+    const std::function<bool(const int )> & pre_collapse,
+    const std::function<void(const int , const bool)> & post_collapse,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
     Eigen::MatrixXd & U,
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J,

include/igl/per_vertex_point_to_plane_quadrics.cpp

@@ -0,0 +1,141 @@
+#include "per_vertex_point_to_plane_quadrics.h"
+#include <Eigen/QR>
+#include <cassert>
+
+
+IGL_INLINE void igl::per_vertex_point_to_plane_quadrics(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXi & EMAP,
+  const Eigen::MatrixXi & EF,
+  const Eigen::MatrixXi & EI,
+  std::vector<
+    std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> > & quadrics)
+{
+  using namespace std;
+  typedef std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> Quadric;
+  const int dim = V.cols();
+  //// Quadrics per face
+  //std::vector<Quadric> face_quadrics(F.rows());
+  // Initialize each vertex quadric to zeros
+  quadrics.resize(
+    V.rows(),{Eigen::MatrixXd::Zero(dim,dim),Eigen::RowVectorXd::Zero(dim),0});
+  Eigen::MatrixXd I = Eigen::MatrixXd::Identity(dim,dim);
+  // Rather initial with zeros, initial with a small amount of energy pull
+  // toward original vertex position
+  const double w = 1e-10;
+  for(int v = 0;v<V.rows();v++)
+  {
+    std::get<0>(quadrics[v]) = w*I;
+    Eigen::RowVectorXd Vv = V.row(v);
+    std::get<1>(quadrics[v]) = w*-Vv;
+    std::get<2>(quadrics[v]) = w*Vv.dot(Vv);
+  }
+  // Generic nD qslim from "Simplifying Surfaces with Color and Texture
+  // using Quadric Error Metric" (follow up to original QSlim)
+  for(int f = 0;f<F.rows();f++)
+  {
+    int infinite_corner = -1;
+    for(int c = 0;c<3;c++)
+    {
+      if(isinf(V(F(f,c),0)) || isinf(V(F(f,c),1)) || isinf(V(F(f,c),2)))
+      {
+        assert(infinite_corner == -1 && "Should only be one infinite corner");
+        infinite_corner = c;
+      }
+    }
+    // Inputs:
+    //   p  1 by n row point on the subspace 
+    //   S  m by n matrix where rows coorespond to orthonormal spanning
+    //     vectors of the subspace to which we're measuring distance (usually
+    //     a plane, m=2)
+    //   weight  scalar weight
+    // Returns quadric triple {A,b,c} so that A-2*b+c measures the quadric
+    const auto subspace_quadric = [&I](
+      const Eigen::RowVectorXd & p,
+      const Eigen::MatrixXd & S,
+      const double  weight)->Quadric
+    {
+      // Dimension of subspace
+      const int m = S.rows();
+      // Weight face's quadric (v'*A*v + 2*b'*v + c) by area
+      // e1 and e2 should be perpendicular
+      Eigen::MatrixXd A = I;
+      Eigen::RowVectorXd b = -p;
+      double c = p.dot(p);
+      for(int i = 0;i<m;i++)
+      {
+        Eigen::RowVectorXd ei = S.row(i);
+        for(int j = 0;j<i;j++) assert(std::abs(S.row(j).dot(ei)) < 1e-10);
+        A += -ei.transpose()*ei;
+        b += p.dot(ei)*ei;
+        c += -pow(p.dot(ei),2);
+      }
+      return { weight*A, weight*b, weight*c };
+    };
+    if(infinite_corner == -1)
+    {
+      // Finite (non-boundary) face
+      Eigen::RowVectorXd p = V.row(F(f,0));
+      Eigen::RowVectorXd q = V.row(F(f,1));
+      Eigen::RowVectorXd r = V.row(F(f,2));
+      Eigen::RowVectorXd pq = q-p;
+      Eigen::RowVectorXd pr = r-p;
+      // Gram Determinant = squared area of parallelogram 
+      double area = sqrt(pq.squaredNorm()*pr.squaredNorm()-pow(pr.dot(pq),2));
+      Eigen::RowVectorXd e1 = pq.normalized();
+      Eigen::RowVectorXd e2 = (pr-e1.dot(pr)*e1).normalized();
+      Eigen::MatrixXd S(2,V.cols());
+      S<<e1,e2;
+      Quadric face_quadric = subspace_quadric(p,S,area);
+      // Throw at each corner
+      for(int c = 0;c<3;c++)
+      {
+        quadrics[F(f,c)] = quadrics[F(f,c)] + face_quadric;
+      }
+    }else
+    {
+      // cth corner is infinite --> edge opposite cth corner is boundary
+      // Boundary edge vector
+      const Eigen::RowVectorXd p = V.row(F(f,(infinite_corner+1)%3));
+      Eigen::RowVectorXd ev = V.row(F(f,(infinite_corner+2)%3)) - p;
+      const double length = ev.norm();
+      ev /= length;
+      // Face neighbor across boundary edge
+      int e = EMAP(f+F.rows()*infinite_corner);
+      int opp = EF(e,0) == f ? 1 : 0;
+      int n =  EF(e,opp);
+      int nc = EI(e,opp);
+      assert(
+        ((F(f,(infinite_corner+1)%3) == F(n,(nc+1)%3) && 
+          F(f,(infinite_corner+2)%3) == F(n,(nc+2)%3)) || 
+          (F(f,(infinite_corner+1)%3) == F(n,(nc+2)%3) 
+          && F(f,(infinite_corner+2)%3) == F(n,(nc+1)%3))) && 
+        "Edge flaps not agreeing on shared edge");
+      // Edge vector on opposite face
+      const Eigen::RowVectorXd eu = V.row(F(n,nc)) - p;
+      assert(!isinf(eu(0)));
+      // Matrix with vectors spanning plane as columns
+      Eigen::MatrixXd A(ev.size(),2);
+      A<<ev.transpose(),eu.transpose();
+      // Use QR decomposition to find basis for orthogonal space
+      Eigen::HouseholderQR<Eigen::MatrixXd> qr(A);
+      const Eigen::MatrixXd Q = qr.householderQ();
+      const Eigen::MatrixXd N = 
+        Q.topRightCorner(ev.size(),ev.size()-2).transpose();
+      assert(N.cols() == ev.size());
+      assert(N.rows() == ev.size()-2);
+      Eigen::MatrixXd S(N.rows()+1,ev.size());
+      S<<ev,N;
+      Quadric boundary_edge_quadric = subspace_quadric(p,S,length);
+      for(int c = 0;c<3;c++)
+      {
+        if(c != infinite_corner)
+        {
+          quadrics[F(f,c)] = quadrics[F(f,c)] + boundary_edge_quadric;
+        }
+      }
+    }
+  }
+}
+

include/igl/per_vertex_point_to_plane_quadrics.h

@@ -0,0 +1,46 @@
+#ifndef IGL_PER_VERTEX_POINT_TO_PLANE_QUADRICS_H
+#define IGL_PER_VERTEX_POINT_TO_PLANE_QUADRICS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+#include <tuple>
+namespace igl
+{
+  // Compute quadrics per vertex of a "closed" triangle mesh (V,F). Rather than
+  // follow the qslim paper, this implements the lesser-known _follow up_
+  // "Simplifying Surfaces with Color and Texture using Quadric Error Metrics".
+  // This allows V to be n-dimensional (where the extra coordiantes store
+  // texture UVs, color RGBs, etc.
+  //
+  // Inputs:
+  //   V  #V by n list of vertex positions. Assumes that vertices with
+  //     infinite coordinates are "points at infinity" being used to close up
+  //     boundary edges with faces. This allows special subspace quadrice for
+  //     boundary edges: There should never be more than one "point at
+  //     infinity" in a single triangle.
+  //   F  #F by 3 list of triangle indices into V
+  //   E  #E by 2 list of edge indices into V.
+  //   EMAP #F*3 list of indices into E, mapping each directed edge to unique
+  //     unique edge in E
+  //   EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  //     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).
+  // Outputs:
+  //   quadrics  #V list of quadrics, where a quadric is a tuple {A,b,c} such
+  //     that the quadratic energy of moving this vertex to position x is
+  //     given by x'Ax - 2b + c
+  //
+  IGL_INLINE void per_vertex_point_to_plane_quadrics(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
+    std::vector<
+      std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> > & quadrics);
+};
+#ifndef IGL_STATIC_LIBRAY
+#  include "per_vertex_point_to_plane_quadrics.cpp"
+#endif
+#endif

include/igl/qslim.cpp

@@ -0,0 +1,76 @@
+#include "qslim.h"
+
+#include "collapse_edge.h"
+#include "connect_boundary_to_infinity.h"
+#include "decimate.h"
+#include "edge_flaps.h"
+#include "max_faces_stopping_condition.h"
+#include "per_vertex_point_to_plane_quadrics.h"
+#include "qslim_optimal_collapse_edge_callbacks.h"
+#include "quadric_binary_plus_operator.h"
+#include "remove_unreferenced.h"
+#include "slice.h"
+#include "slice_mask.h"
+
+IGL_INLINE bool igl::qslim(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const size_t max_m,
+  Eigen::MatrixXd & U,
+  Eigen::MatrixXi & G,
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & I)
+{
+  using namespace igl;
+  // Original number of faces
+  const int orig_m = F.rows();
+  // Tracking number of faces
+  int m = F.rows();
+  typedef Eigen::MatrixXd DerivedV;
+  typedef Eigen::MatrixXi DerivedF;
+  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);
+  // Quadrics per vertex
+  typedef std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> Quadric;
+  std::vector<Quadric> quadrics;
+  per_vertex_point_to_plane_quadrics(VO,FO,EMAP,EF,EI,quadrics);
+  // State variables keeping track of edge we just collapsed
+  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 int)> pre_collapse;
+  std::function<void(const int,const bool)> post_collapse;
+  qslim_optimal_collapse_edge_callbacks(
+    E,quadrics,v1,v2, cost_and_placement, pre_collapse,post_collapse);
+  // Call to greedy decimator
+  bool ret = decimate(
+    VO, FO,
+    cost_and_placement,
+    max_faces_stopping_condition(m,orig_m,max_m),
+    pre_collapse,
+    post_collapse,
+    E, EMAP, EF, EI,
+    U, G, J, I);
+  // Remove phony boundary faces and clean up
+  const Eigen::Array<bool,Eigen::Dynamic,1> keep = (J.array()<orig_m);
+  igl::slice_mask(Eigen::MatrixXi(G),keep,1,G);
+  igl::slice_mask(Eigen::VectorXi(J),keep,1,J);
+  Eigen::VectorXi _1,I2;
+  igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1,I2);
+  igl::slice(Eigen::VectorXi(I),I2,1,I);
+  return ret;
+}

include/igl/qslim.h

@@ -0,0 +1,34 @@
+#ifndef IGL_QSLIM_H
+#define IGL_QSLIM_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+
+  // Decimate (simplify) a triangle mesh in nD according to the paper
+  // "Simplifying Surfaces with Color and Texture using Quadric Error Metrics"
+  // by [Garland and Heckbert, 1987] (technically a followup to qslim). The
+  // mesh can have open boundaries but should be edge-manifold.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions. Assumes that vertices w
+  //   F  #F by 3 list of triangle indices into V
+  //   max_m  desired number of output faces
+  // Outputs:
+  //   U  #U by dim list of output vertex posistions (can be same ref as V)
+  //   G  #G by 3 list of output face indices into U (can be same ref as G)
+  //   J  #G list of indices into F of birth face
+  //   I  #U list of indices into V of birth vertices
+  IGL_INLINE bool qslim(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const size_t max_m,
+    Eigen::MatrixXd & U,
+    Eigen::MatrixXi & G,
+    Eigen::VectorXi & J,
+    Eigen::VectorXi & I);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "qslim.cpp"
+#endif
+#endif

include/igl/qslim_optimal_collapse_edge_callbacks.cpp

@@ -0,0 +1,68 @@
+#include "qslim_optimal_collapse_edge_callbacks.h"
+
+IGL_INLINE void igl::qslim_optimal_collapse_edge_callbacks(
+  Eigen::MatrixXi & E,
+  std::vector<std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> > & 
+    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 int)> & pre_collapse,
+  std::function<void(const int,const bool)> & post_collapse)
+{
+  typedef std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> Quadric;
+  cost_and_placement = [&quadrics,&v1,&v2](
+    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)
+  {
+    // Combined quadric
+    Quadric quadric_p;
+    quadric_p = quadrics[E(e,0)] + quadrics[E(e,1)];
+    // Quadric: p'Ap + 2b'p + c
+    // optimal point: Ap = -b, or rather because we have row vectors: pA=-b
+    const auto & A = std::get<0>(quadric_p);
+    const auto & b = std::get<1>(quadric_p);
+    const auto & c = std::get<2>(quadric_p);
+    p = -b*A.inverse();
+    cost = p.dot(p*A) + 2*p.dot(b) + c;
+    // Force infs and nans to infinity
+    if(isinf(cost) || cost!=cost)
+    {
+      cost = std::numeric_limits<double>::infinity();
+    }
+  };
+  // Remember endpoints
+  pre_collapse = [&v1,&v2,&E](const int e)->bool
+  {
+    v1 = E(e,0);
+    v2 = E(e,1);
+    return true;
+  };
+  // update quadric
+  post_collapse = [&v1,&v2,&quadrics](
+    const int e, 
+    const bool collapsed)
+  {
+    if(collapsed)
+    {
+      quadrics[v1<v2?v1:v2] = quadrics[v1] + quadrics[v2];
+    }
+  };
+}
+

include/igl/qslim_optimal_collapse_edge_callbacks.h

@@ -0,0 +1,46 @@
+#ifndef IGL_QSLIM_OPTIMAL_COLLAPSE_EDGE_CALLBACKS_H
+#define IGL_QSLIM_OPTIMAL_COLLAPSE_EDGE_CALLBACKS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <functional>
+#include <vector>
+#include <tuple>
+namespace igl
+{
+
+  // Prepare callbacks for decimating edges using the qslim optimal placement
+  // metric.
+  //
+  // Inputs:
+  //   E  #E by 2 list of working edges
+  //   quadrics  reference to list of working per vertex quadrics 
+  //   v1  working variable to maintain end point of collapsed edge
+  //   v2  working variable to maintain end point of collapsed edge
+  // Outputs
+  //   cost_and_placement  callback for evaluating cost of edge collapse and
+  //     determining placement of vertex (see collapse_edge)
+  //   pre_collapse  callback before edge collapse (see collapse_edge)
+  //   post_collapse  callback after edge collapse (see collapse_edge)
+  IGL_INLINE void qslim_optimal_collapse_edge_callbacks(
+    Eigen::MatrixXi & E,
+    std::vector<std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> > & 
+      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 int)> & pre_collapse,
+    std::function<void(const int,const bool)> & post_collapse);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "qslim_optimal_collapse_edge_callbacks.cpp"
+#endif
+#endif

include/igl/quadric_binary_plus_operator.cpp

@@ -0,0 +1,17 @@
+#include "quadric_binary_plus_operator.h"
+
+IGL_INLINE std::tuple< Eigen::MatrixXd, Eigen::RowVectorXd, double> 
+  igl::operator+(
+    const std::tuple< Eigen::MatrixXd, Eigen::RowVectorXd, double>  & a, 
+    const std::tuple< Eigen::MatrixXd, Eigen::RowVectorXd, double>  & b)
+{
+  std::tuple<
+    Eigen::MatrixXd,
+    Eigen::RowVectorXd,
+    double>  c;
+  std::get<0>(c) = (std::get<0>(a) + std::get<0>(b)).eval();
+  std::get<1>(c) = (std::get<1>(a) + std::get<1>(b)).eval();
+  std::get<2>(c) = (std::get<2>(a) + std::get<2>(b));
+  return c;
+}
+

include/igl/quadric_binary_plus_operator.h

@@ -0,0 +1,28 @@
+#ifndef IGL_QUADRIC_BINARY_PLUS_OPERATOR_H
+#define IGL_QUADRIC_BINARY_PLUS_OPERATOR_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <tuple>
+
+namespace igl
+{
+  // A binary addition operator for Quadric tuples compatible with qslim,
+  // computing c = a+b
+  //
+  // Inputs:
+  //   a  QSlim quadric
+  //   b  QSlim quadric
+  // Output
+  //   c  QSlim quadric
+  //
+  IGL_INLINE std::tuple< Eigen::MatrixXd, Eigen::RowVectorXd, double> 
+    operator+(
+      const std::tuple< Eigen::MatrixXd, Eigen::RowVectorXd, double>  & a, 
+      const std::tuple< Eigen::MatrixXd, Eigen::RowVectorXd, double>  & b);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "quadric_binary_plus_operator.cpp"
+#endif
+
+#endif

tutorial/703_Decimation/main.cpp

@@ -75,7 +75,11 @@ 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))
+        const auto & always_try = [](const int e)->bool{return true;};
+        const auto & never_care = [](const int e, const bool collapsed){};
+        if(!collapse_edge(
+          shortest_edge_and_midpoint,always_try,never_care,
+          V,F,E,EMAP,EF,EI,Q,Qit,C))
         {
           break;
         }