Posts Tagged ‘qslim’

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.

Mesh decimation (aka simplification) in matlab

Thursday, October 15th, 2015

I’d forgotten that I’d discovered that matlab has a built-in function for triangle mesh simplification: reduce_patch. Here’s a little comparison between this method and the ones I’ve wrapped up in gptoolbox.

Input mesh

Given an input mesh with vertices V and faces F (this elephant has 14732 faces):

cartoon-elephant-input.gif

Built-in matlab reduce patch

We can decimate the mesh to 1000 faces using pure built-in matlab:

[mF,mV] = reducepatch(F,V,1000);

cartoon-elephant-decimated-matlab.gif

Notice that it’s a rather adaptive remeshing. There aren’t (m)any options for controlling reducepatch, just the 'fast' flag, but in this case it will produce an identical output.

Libigl’s decimate

Libigl has a framework for edge-collapse decimation on manifold meshes and a default setting for regular meshing:

[iV,iF] = decimate_libigl(V,F,1000);

cartoon elephant decimated libigl

CGAL’s decimation

I have a wrapper for CGAL’s decimation algorithms. These can result in a regular remeshing with 1000 faces:

[c1V,c1F] = decimate_cgal(V,F,1000/size(F,1));

cartoon-elephant-decimated-cgal-regular.gif

or an adaptive meshing:

[c2V,c2F] = decimate_cgal(V,F,1000/size(F,1),'Adaptive',true);

cartoon-elephant-decimated-cgal-adaptive.gif

CGAL has a fairly general interface for an edge-collapser, but these are the only metrics provided by default from CGAL.

QSlim

cartoon-elephant-decimated-qslim.gif

Timing

Suppose my original elephant had 3.7 million faces instead of just 14K, how do these methods measure up when reducing to 1000 faces:

Method Time
reducepatch 64.4s
libigl 74.56s
cgal (regular) 117.7s
cgal (adaptive) 214.8s
qslim 108.0s + 507.5s + 59.4s

My qslim wrapper is currently very silly. It’s calling qslim as an executable and then reading in the entire collapse log and conducting the collapse within matlab (the second time is for reading the log, the third for conducting the collapse).

None of these methods guarantee that the result is self-intersection free. I’ve noticed that qslim and matlab will create non-manifold meshes from manifold inputs. The current libigl implementation should not creative non-manifold output, but assumes the input is a closed surface. For casual decimation, it seems like matlab’s reducepatch and libigl’s decimate_libigl are the best bets for speedy adaptive and regular meshing respectively. Conveniently they also require the least in terms of dependencies.

Edge collapse and mesh decimation in libigl

Thursday, April 16th, 2015

Libigl purposefully does not build itself around complicated mesh data-structures. This is freeing for many reasons: it’s very easy to include our code in an arbitrary project, functions are not artificially limited to manifold meshes, array based data structures are easily serialized and fast, matrix operations on vertex positions are directly exposed, etc.

mesh decimation in libigl

But our choice is also limiting. In particular, combinatorial changes to the mesh are potentially expensive and difficult to implement. Recently, I took a first stab at implementing an efficient edge-collapse routine for libigl. Indeed I’m seeing good performance: O(1) constant time for a single edge collapse and O(log m) time for cost-priority-queue-based collapses. The data structures are a little “messy” in the sense that when edges or faces disappear their rows are not deleted, but rather just redacted: entries are replaced with NULL values. This is because my approach is array-based and constant resizing would ruin performance. Luckily for edge-collapse, the number of elements only decreases. For edge-split I’ll have to think about amortized costs…

For now, check out the updated tutorial entry for libigl’s new mesh decimation and edge collapse features.

The code is “programmable” in the sense I expose function handles for computing edge costs and merged vertex placements. Though the default functions I currently provide are quite naive, this should support rather advanced “Progressive Meshes” style simplification with creases, sharp features, adaptivity, etc.

qslim for Mac OS X on github

Wednesday, April 15th, 2015

A while ago I patched Michael Garland’s qslim to compile on modern Mac OS X. I’ve finally wrangled in the instructions and modified files to a github repo. There are also wrappers for this tool in my gptoolbox. It’d be cool to have a proper C++ API to this tool, but that’s future work.

Qslim is one of the seminal mesh decimation/edge-collapse/simplification techniques. It’s fast, simple and predictable. It’s a good starting place if you’re trying to simplify your mesh, especially in an automated way.

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