Posts Tagged ‘mesh simplification’

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.

Compiling and running qslim on Mac OS x

Thursday, May 23rd, 2013

I recently successfully compiled and executed qslim on mac os x. It took a little tweaking but here’s how I did it.

fltk

Before anything else, be sure to install fltk. I did this using macports:


sudo port install fltk

libgfx

Then, following the instructions in qslim-2.1/README.txt, I first compiled libgfx:


cd libgfx
./configure

The libpng support seems out of date and I don’t really need it so in raster-png.cxx I replaced:


//#ifdef HAVE_LIBPNG
#if false

memcopy was shoing up undefined so I added the following include to raster.cxx:


#include <cstring>

Similarly there was a missing include in time.cxx:


#elif defined(HAVE_TIMES)
#include <sys/times.h>

Know, while inside libgfx I issue:


make -C src

mixkit

First travel inside mixkit and configure:


cd ../mixkit
./configure

Trying to make I got a bunch of errors of the form:


MxDynBlock.h:44:33: Error: 'resize' was not declared in this scope, and no declarations were found by argument-dependent lookup at the point of instantiation [-fpermissive]
MxDynBlock.h:44:33: Error: note: declarations in dependent base 'MxBlock<MxStdModel*>' are not found by unqualified lookup
MxDynBlock.h:44:33: Error: note: use 'this->resize' instead

There’re at least two ways to fix this. You can add this-> everywhere gcc is tell you to: edit MxStack.h and MxDynBlock.h or you can add -fpermissive to the CXXFLAGS. This is awkwardly done not in mix-config but in ../libgfx/gfx-config. Where you can edit into something like:


CXXFLAGS = -g -O2  -I/Users/ajx/Downloads/qslim-2.1/libgfx/include -DHAVE_CONFIG_H $(WIN_FLAGS) -fpermissive

Alternatively you could have done this when configuring libgfx with:


cd ../libgfx
./configure CXXFLAGS='-fpermissive'
cd ../mixkit/

If you do it this way then you’ll get warnings instead of errors.

Now you can build with:


make -C src

qslim and qvis

Travel to the qslim tools directory:


cd ../tools/qslim

Here I need to add libfltk_gl to the linker commands in Makefile:


        $(CXX) -o qslim $(OBJS) $(LDFLAGS) $(LIBMIX) -lm -lfltk_gl
...
        $(CXX) -o qvis $(QVIS_OBJS) $(LDFLAGS) $(LIBMIX) $(GUI_LIBS) -lm -lfltk_gl

Then build with:


make

(those fpermissive warnings may show up again, but no problem)

Finally I can run the command line program qslim or the gui qvis. For some reason the qui qvis.app immediately closes if I double click on it. But I can run it by issuing:


./qvis.app/Contents/MacOS/qvis 

where I’m then prompted to load a .smf file.

or


./qvis.app/Contents/MacOS/qvis input.smf

You can call the command line program with:


./qslim -t 1000 -M smf -q 2>/dev/null

which spills the .smf file to stdout

It seems .smf format is the same as .obj, at least if there are only vertices and triangular faces. So, you can compose—a rather long—bash oneliner that simplifies an input.obj file to an output.obj file:


grep "^[vf] " input.obj | sed -e "s/^f *\([0-9][0-9]*\)\/[^ ]*  *\([0-9][0-9]*\)\/[^ ]*  *\([0-9][0-9]*\)\/.*/f \1 \2 \3/g" | ./qslim -t 1000 -M smf -q 2>/dev/null | grep "^[vf]" > output.obj

qslim/filters

I also built the “filters” (converters?) with:


cd ../filters
make

Note that ply2smf only acts on stdin and stdout so call it with:


cat input.ply | ./ply2smf >output.smf

Update: I found an easier way to setup the configure script and ensure that all macports libraries are found:


env CPPFLAGS="-I/opt/local/include -fpermissive" LDFLAGS="-L/opt/local/lib" ./configure