Posts Tagged ‘decimation’

Implementing QSlim mesh simplification using libigl

Friday, November 4th, 2016

One of my research projects is looking at mesh simplification and the perennial comparison is always to Garland and Heckbert’s “Surface simplification using quadric error metrics” (aka QSlim). And relevant to us, the lesser-known followup for “Surfaces with Color and Texture using Quadric Error” (aka QSlim in nD)

Currently, libigl has built in edge-collapse based mesh simplification with an exposed function handle: cost_and_placement which determines the cost of collapsing an edge and the placement. This function is used for two purposes: to initialize the heap of edges to collapse and to update edges in that heap after a collapse.

Currently I had only implemented a very simple metric: cost is the length of the edge and collapsed vertices are placed at edge midpoints.

So a bit out of research necessity and a bit just for a challenge, I timed how long it would take for me to implement QSlim in nD.

cactus decimated using qslim implemented in libigl

So, 1 hour and 20 minutes later, here is my QSlim in nD demo based off of the decimation entry in the libigl tutorial:

#include <igl/collapse_edge.h>
#include <igl/edge_flaps.h>
#include <igl/read_triangle_mesh.h>
#include <igl/viewer/Viewer.h>
#include <Eigen/Core>
#include <iostream>
#include <set>

int main(int argc, char * argv[])
{
  using namespace std;
  using namespace Eigen;
  using namespace igl;
  cout<<"Usage: ./703_Decimation_bin [filename.(off|obj|ply)]"<<endl;
  cout<<"  [space]  toggle animation."<<endl;
  cout<<"  'r'  reset."<<endl;
  // Load a closed manifold mesh
  string filename(argv[1]);
  if(argc>=2)
  {
    filename = argv[1];
  }
  MatrixXd V,OV;
  MatrixXi F,OF;
  read_triangle_mesh(filename,OV,OF);

  igl::viewer::Viewer viewer;

  // 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;
  // If an edge were collapsed, we'd collapse it to these points:
  MatrixXd C;
  int num_collapsed;


  // Function for computing cost of collapsing edge (lenght) and placement
  // (midpoint)
  const auto & shortest_edge_and_midpoint = [](
    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,
    RowVectorXd & p)
  {
    cost = (V.row(E(e,0))-V.row(E(e,1))).norm();
    p = 0.5*(V.row(E(e,0))+V.row(E(e,1)));
  };

  // Quadrics per vertex
  typedef std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> Quadric;
  std::vector<Quadric> quadrics;

  // c **is** allowed to be a or b.
  const auto & plus = [](const Quadric & a, const Quadric & b, Quadric & 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));
  };
  // State variables keeping track of whether we've just collpased edge (v1,v2)
  int v1 = -1;
  int v2 = -1;
  const auto & qslim_optimal = [&quadrics,&plus,&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,
    RowVectorXd & p)
  {
    // Then we just collapsed (v1,v2)
    if(v1>=0 && v2>=0)
    {
      plus(quadrics[v1],quadrics[v2],quadrics[v1<v2?v1:v2]);
      v1 = -1;
      v2 = -1;
    }
    // Combined quadric
    Quadric quadric_p;
    plus(quadrics[E(e,0)],quadrics[E(e,1)],quadric_p);
    // 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;
  };

  // Function to reset original mesh and data structures
  const auto & reset = [&]()
  {
    F = OF;
    V = OV;
    edge_flaps(F,E,EMAP,EF,EI);
    Qit.resize(E.rows());
    Q.clear();

    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++)
    {
      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();
      // e1 and e2 be perpendicular
      assert(std::abs(e1.dot(e2)) < 1e-10);
      // Weight face's quadric (v'*A*v + 2*b'*v + c) by area
      const Eigen::MatrixXd A = I-e1.transpose()*e1-e2.transpose()*e2;
      const Eigen::RowVectorXd b = p.dot(e1)*e1 + p.dot(e2)*e2 - p;
      const double c = (p.dot(p) - pow(p.dot(e1),2) - pow(p.dot(e2),2));
      face_quadrics[f] = { area*A, area*b, area*c };
      // Throw at each corner
      for(int c = 0;c<3;c++)
      {
        plus(
          quadrics[F(f,c)],
          face_quadrics[f],
          quadrics[F(f,c)]);
      }
    }




    C.resize(E.rows(),V.cols());
    VectorXd costs(E.rows());
    v1 = -1;
    v2 = -1;
    for(int e = 0;e<E.rows();e++)
    {
      double cost = e;
      RowVectorXd p(1,3);
      qslim_optimal(e,V,F,E,EMAP,EF,EI,cost,p);
      C.row(e) = p;
      Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
    }
    num_collapsed = 0;
    viewer.data.clear();
    viewer.data.set_mesh(V,F);
    viewer.data.set_face_based(true);
  };

  const auto &pre_draw = [&](igl::viewer::Viewer & viewer)->bool
  {
    // If animating then collapse 10% of edges
    if(viewer.core.is_animating && !Q.empty())
    {
      bool something_collapsed = false;
      // collapse edge
      const int max_iter = std::ceil(0.01*Q.size());
      for(int j = 0;j<max_iter;j++)
      {
        std::pair<double,int> p = *(Q.begin());
        bool shouldnt_collapse = false;
        if((! Q.empty()) && 
            (p.first != std::numeric_limits<double>::infinity()))
        {
          v1 = E(p.second,0);
          v2 = E(p.second,1);
        }else
        {
          v1 = -1;
          v2 = -1;
          shouldnt_collapse = true;
        }
        if(!collapse_edge(qslim_optimal,V,F,E,EMAP,EF,EI,Q,Qit,C))
        {
          break;
        }else
        {
          // we just collapsed. 
          assert(!shouldnt_collapse);
        }
        something_collapsed = true;
        num_collapsed++;
      }

      if(something_collapsed)
      {
        viewer.data.clear();
        viewer.data.set_mesh(V,F);
        viewer.data.set_face_based(true);
      }
    }
    return false;
  };

  const auto &key_down =
    [&](igl::viewer::Viewer &viewer,unsigned char key,int mod)->bool
  {
    switch(key)
    {
      case ' ':
        viewer.core.is_animating ^= 1;
        break;
      case 'R':
      case 'r':
        reset();
        break;
      default:
        return false;
    }
    return true;
  };

  reset();
  viewer.core.is_animating = true;
  viewer.callback_key_down = key_down;
  viewer.callback_pre_draw = pre_draw;
  return viewer.launch();
}

I’m pretty happy that this didn’t take me all day, but admittedly it should have been faster. One confusing part was that the current API in libigl for edge collapsing doesn’t have an explicit function handle that’s called when an edge is successfully collapsed. Detecting this is necessary to propagate the quadric to the collapsed vertex (but isn’t necessary for simple costs like edge-length). I faked this for now using v1 and v2 as state variables. Probably, I’ll change the libigl API to accommodate this. The other time sink was worry about degenerate Quadrics (when the system matrix is singular). Rather than try to detect this categorically, I opted to sprinkle in a little energy pulling vertices toward their (ancestors’) original positions. This seems to work fairly well.

Nested Cages project page

Friday, October 2nd, 2015

nested cages bunny

We’ve posted a project page for our upcoming SIGGRAPH Asia paper Nested Cages, a collaboration between Leonardo Sacht, Etienne Vouga and myself.

Abstract: Many tasks in geometry processing and physical simulation benefit from multiresolution hierarchies. One important characteristic across a variety of applications is that coarser layers strictly encage finer layers, nesting one another. Existing techniques such as surface mesh decimation, voxelization, or contouring distance level sets do not provide sufficient control over the quality of the output surfaces while maintaining strict nesting. We propose a solution that enables use of application-specific decimation and quality metrics. The method constructs each next-coarsest level of the hierarchy, using a sequence of decimation, flow, and contact-aware optimization steps. From coarse to fine, each layer then fully encages the next while retaining a snug fit. The method is applicable to a wide variety of shapes of complex geometry and topology. We demonstrate the effectiveness of our nested cages not only for multigrid solvers, but also for conservative collision detection, domain discretization for elastic simulation, and cage-based geometric modeling.