## 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.

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/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;

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)));
};

// c **is** allowed to be a or b.
{
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)
{
v1 = -1;
v2 = -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;
};

// 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();
// Initialize each vertex quadric to zeros
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++)
{
Eigen::RowVectorXd Vv = V.row(v);
}
// Generic nD qslim from "Simplifying Surfaces with Color and Texture
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(
}
}

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