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