Posts Tagged ‘libigl’

Eigen performance gotcha calling non-templated function from templated one

Tuesday, July 25th, 2017

I just spent a while tracking down a rather surprising performance bug in my code.

Here’s a minimal example:

#include <Eigen/Dense>
#include <iostream>

int simple_size(const Eigen::MatrixXi & Z)
  return Z.size();

template <typename T> int templated_size(const Eigen::MatrixBase<T> & Y)
  return simple_size(Y);

int main(int argc, const char * argv[])
  const int s = 40000;
  Eigen::MatrixXi X = Eigen::MatrixXi::Zero(40000,40000);
  std::cout<<(X.size()         ?"done":"")<<std::endl;
  std::cout<<(simple_size(X)          ?"done":"")<<std::endl;


Running this, it will show that the last call to templated_size is taking way too long. Inspection will show that a copy of Y is being created to create a Eigen::MatrixXi & reference.

Now, clearly it’s poor design to call a function expecting a Eigen::MatrixXi & reference with a generic templated type Eigen::MatrixBase<T> &, but unfortunately this happens quite often with legacy libigl functions. My expectation was that since T is Eigen::MatrixXi in this case a simple reference would be passed.

It’s worth noting that const is actually creating/hiding the problem. Because simple_size takes a const reference, the compiler is happy to construct a Eigen::MatrixXi on the fly to create a valid reference. Without the consts the compiler stops at an error.

Inflate Wire Mesh in libigl C++ or gptoolbox MATLAB

Wednesday, July 12th, 2017

For a visualization and 3D printing, it’s often useful to “inflate” a edge-network into a thickened surface mesh. One method to do this is described “Sculptural Forms from Hyperbolic Tessellations” by George C Hart. This method works by adding rotated polygons at the ends of each edge offset a bit from the vertices. Then for each vertex the convex hull of incident edges’ polygons is computed and unioned with the convex hull of the polygons at either end of each edge. Hart writes that polygons shared by “edge hulls” and “vertex hulls” can simply be discarded. This is unfortunately not true, in general. It’s not super easier to categorize which faces can be discarded (even in general position) since the answer depends on the thickness, the number of sides of the polygons, their rotations, their offsets, and the angle between neighbouring edges. Fortunately, libigl is very good at conducting unions. We can just conduct the union explicitly and exactly using libigl.

I’ve written a new function for libigl igl::wire_mesh that takes in a wire network and spits out a solid (i.e., closed, watertight, manifold) mesh of a the inflated surface.

I’ve also wrapped this up in a Matlab Mex function in gptooolbox wire_mesh.

Offset surface of triangle mesh in matlab

Wednesday, March 22nd, 2017

Here’s a little demonstration of how to use gptoolbox and MATLAB to generate an offset surfaces of a triangle mesh. This takes a mesh in V,F and creates a mesh SV,SF of the isosurface at signed distance iso:

% Extract offset at minus 3% of bounind box diagonal length
iso = -0.03;
% Resolution grid → resolution of extracted offset surface
side = 60;
% Amount of smoothing to apply to distance field
sigma = 1.4;
bbd = norm(max(V)-min(V));
% Pad generously for positive iso values 
[BC,side,r] = voxel_grid([V;max(V)+iso*1;min(V)-iso*1],side);
D = signed_distance(BC,V,F);
D = reshape(D,side([2 1 3]));
% Smooth signed distance field
D = imfilter(D,fspecial('gaussian',9,sigma),'replicate');
BC3 = reshape(BC,[side([2 1 3]) 3]);
% Use matlab's built-in marching cubes iso-surface mesher (works most of the time)
surf = isosurface(BC3(:,:,:,1),BC3(:,:,:,2),BC3(:,:,:,3),D,iso*bbd);
SV = surf.vertices;
SF = surf.faces;

Here’s a blue bunny with a positive offset surface, an orange “cage”:

bunny with positive offset surface

Here’s a blue bunny with a negative offset surface. This is useful for hollowing out objects to 3d print:

bunny with positive offset surface

Because the iso-surface extraction will over tesselate low curvature patches of the output surface, it would make a lot of sense to remesh/decimate this mesh.

(to create these fancy renderings:)

hold on;
t  = tsurf(F,V,'EdgeColor','none',fsoft,  'FaceVertexCData',repmat(blue,size(V,1),1),'FaceAlpha',1+(iso<0)*(0.35-1),fphong);
ts = tsurf(SF,SV,'EdgeAlpha',0.2+(iso<0)*(0-0.2),fsoft,'FaceVertexCData',repmat(orange,size(SV,1),1),fphong,'FaceAlpha',1+(iso>0)*(0.2-1));
hold off;
axis equal;
t.SpecularStrength = 0.04;
l = light('Position',[5 -5 10],'Style','local');
add_shadow(t,l,'Color',0.8*[1 1 1],'Fade','local','Ground',[0 0 -1 min([V(:,3);SV(:,3)])]);
set(gca,'pos',[0 0 1 1])

Command line program to view 3d meshes from files and piped input

Thursday, December 15th, 2016

Here’s a little C++ program to directly render meshes in 3D viewer from the command line.

This let’s you write little test programs without worrying about linking to a 3D viewer. You just need to output a mesh in a standard format. For example, here’s a tiny program that outputs a cube in an .off format:

#include <igl/read_triangle_mesh.h>
#include <Eigen/Core>
#include <iostream>

int main(int argc, char * argv[])
  using namespace Eigen;
  MatrixXd V(8,3);
  MatrixXi Q(6,4);
    "OFF\n"<<V.rows()<<" "<<Q.rows()<<" 0\n"<<
    V.format(IOFormat(FullPrecision,DontAlignCols," ","\n","","","","\n"))<<
    (Q.array()).format(IOFormat(FullPrecision,DontAlignCols," ","\n","4 ","","","\n"));
  return EXIT_SUCCESS;

Compile this into cube_off then issue:

./cube_off | view mesh

**Update: ** And here’s a funny, little one-liner you can call from matlab to display a mesh via the .obj format:

system(sprintf('echo \"%s%s\" | /usr/local/bin/viewmesh',sprintf('v %0.17f %0.17f %0.17f\n',V'),sprintf('f %d %d %d\n',F')))

Background computation threads with igl::viewer::Viewer

Sunday, December 4th, 2016

Here’s a minimal example showing how to launch background computation threads for each mesh in a list of meshes. Meanwhile the main thread runs a mesh viewer with all meshes concatenated into one huge multi-component mesh. Whenever a computation thread signals that an update to the mesh needs to be made, the main thread will re-concatenate the meshes and update the viewer. In this example, the “computation” is determining how much to move a clock “hand” (random mesh).

Here’s the program running on the cow, cheburashka and knight models:

// Tiny example to demonstrate spawning a background computation thread for
// each mesh in a list and update the viewer when computation results in a
// change to (one of) the meshes.
// In this example, three meshes are read in and interpreted as "hands" of a
// clock. The "computation" is just a busy-wait until the hand should move
// (after one second, one minute, one hour). This is, of course, a terrible way
// to implement a clock.
// ./test libigl/tutorial/shared/{,,}
#include <igl/read_triangle_mesh.h>
#include <igl/point_mesh_squared_distance.h>
#include <igl/combine.h>
#include <igl/viewer/Viewer.h>
#include <Eigen/Geometry>
#include <thread>
int main(int argc, char * argv[])
  using namespace std;
  using namespace Eigen;
  // Read in k meshes (<=3 will be used)
  std::vector<Eigen::MatrixXd> VV(argc-1);
  std::vector<Eigen::MatrixXi> FF(argc-1);
  for(int i = 1;i<argc;i++)
    VV[i-1].col(0).array() -= VV[i-1].col(0).mean();
    VV[i-1].col(1).array() -= VV[i-1].col(1).minCoeff();
  bool continue_computing = true;
  // Flag to redraw and mutex to guard it
  bool redraw = false;
  std::mutex m;
  // After waiting `tic` seconds, rotate `V` by `deg` degrees
  const auto rot = [&continue_computing,&redraw,&m](
    const double tic, const double deg, Eigen::MatrixXd& V)
      // Let's watch this at 500x: Real clocks are slow.
      V *= Eigen::AngleAxisd(deg/180.*igl::PI,
        std::lock_guard<std::mutex> lock(m);
        redraw = true;
  // Launch background "computation" threads for each "hand" of the clock
  // std::ref is important, otherwise std::bind passes by copy
  std::vector<std::thread> threads;
    case 3:
    case 2:
    case 1:
    default: break;
  igl::viewer::Viewer v;;
  // Helper function to view k meshes
  const auto & set_meshes = [](
    const std::vector<Eigen::MatrixXd> & VV,
    const std::vector<Eigen::MatrixXi> & FF,
    igl::viewer::Viewer & v)
    Eigen::MatrixXd V;
    Eigen::MatrixXi F;
  // Continuous draw loop. TODO: trigger draws using glfwPostEmptyEvent
  v.core.is_animating = true;
  // Before every draw check if the meshes have changed. 
  v.callback_pre_draw = 
    [&redraw,&m,&VV,&FF,&set_meshes](igl::viewer::Viewer & v)->bool
          std::lock_guard<std::mutex> lock(m);
          redraw = false;
      return false;
  // Tell computation threads to stop
  continue_computing = false;
  // Join with computation threads: return to serial main thread.
  for(auto & t : threads) if(t.joinable()) t.join();

mesh clock igl viewer threads

This is pretty hacky, but it will allow me to use the standard libigl viewer for making comparisons of mesh-editing methods whose performances are very different.

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

  // 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)
      v1 = -1;
      v2 = -1;
    // Combined quadric
    Quadric 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 =*A) + 2* + c;

  // Function to reset original mesh and data structures
  const auto & reset = [&]()
    F = OF;
    V = OV;

    const int dim = V.cols();
    // Quadrics per face
    std::vector<Quadric> face_quadrics(F.rows());
    // Initialize each vertex quadric to zeros
    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*;
    // 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(,2));
      Eigen::RowVectorXd e1 = pq.normalized();
      Eigen::RowVectorXd e2 = (*e1).normalized();
      // e1 and e2 be perpendicular
      assert(std::abs( < 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 =*e1 +*e2 - p;
      const double c = ( - pow(,2) - pow(,2));
      face_quadrics[f] = { area*A, area*b, area*c };
      // Throw at each corner
      for(int c = 0;c<3;c++)

    VectorXd costs(E.rows());
    v1 = -1;
    v2 = -1;
    for(int e = 0;e<E.rows();e++)
      double cost = e;
      RowVectorXd p(1,3);
      C.row(e) = p;
      Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
    num_collapsed = 0;;,F);;

  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);
          v1 = -1;
          v2 = -1;
          shouldnt_collapse = true;
          // we just collapsed. 
        something_collapsed = true;

    return false;

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

  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.

Adding libigl as a submodule

Thursday, June 9th, 2016

If you’re adding libigl as a submodule to your project, you should add it recursively so that all of its submodules are also pulled in:

git submodule add
git submodule update --init --recursive

Mesh Arrangements for Solid Geometry preprint

Thursday, April 21st, 2016

mesh booleans

I’m very excited to publish my work with Qingnan Zhou, Eitan Grinspun and Denis Zorin on mesh booleans at this year’s SIGGRAPH. The paper is titled Mesh Arrangements for Solid Geometry. The code is (and has been) in libigl and gptoolbox and pymesh.

Variadic mesh boolean operations in libigl and gptoolbox

Thursday, April 14th, 2016

I’ve just pushed some new changes to libigl and gptoolbox to expose the variadic implementation of our robust mesh boolean operations. By variadic, I mean that the boolean function can take one mesh as input or two meshes as input (usual binary case) or three or four and so on. This means you can easily take the union or intersection of n objects, but it also means that complex operations can be reduced to single call. For example, identifying regions inside at least k out of n objects. I first became aware of this variadic concept when reading “QuickCSG: Arbitrary and Faster Boolean Combinations of N Solids”.

In libigl, you can call igl::copyleft::cgal::mesh_boolean with an n-long list of mesh vertex arrays and an n-long list of mesh face arrays. Rather than the usual union, intersection, minus etc., you can also pass a C++ function handle. This function handle will filter the winding number vector at any point in space. For example, if you’d like to extract the region inside object 0 and not inside object 1 and either inside object 2 or not inside object 3. The filter would return w(0) > 0 && w(1)<=0 && (w(2)>0 || w(3)<= 0).

After a bit of pointer tweaking, I’ve also exposed this interface to the matlab wrapper in gptoolbox. You can pass a matlab function handle after the optional argument ‘WindingNumberFilter’. For example, assuming V and F are cell arrays containing vertex position arrays and face indices respectively, the following command will extract the region inside at least 3 of the inputs:

 [VV,FF,J] = mesh_boolean(V,F,'','WindingNumberFilter',@(w) sum(w>0)>=3);

You could also call it with

 [VV,FF,J] = mesh_boolean(V{1},F{1},V{2},F{2},V{3},F{3},V{4},F{4},'','WindingNumberFilter',@(w) sum(w>0)>=3);

Here’s an example of extracting min-k results from 4 spheres:

four spheres variadic boolean operation

Since extraction from our mesh arrangement is cheap, each operation takes the same amount of time. And fortunately it seems that the mex function handle overhead is not so bad.

Update: Variadic operations are also useful for condensing entire binary CSG trees: (A union B) minus (C union D) could be a single extraction function on (A,B,C,D) return (w(0)>0 || w(1)>0) && !(w(2)>0 || w(3)>0)

Using meshfix “from libigl”

Wednesday, March 30th, 2016

Here’s a little example of how to call meshfix using plain Eigen types to represent a mesh (i.e. in the libigl style. I’ve included this example in my github fork of the meshfix code.

#include "meshfix.h"

#include <igl/read_triangle_mesh.h>
#include <igl/write_triangle_mesh.h>
#include <iostream>

int main(int argc, char * argv[])
  // Load in libigl's (V,F) format
  Eigen::MatrixXd V,W;
  Eigen::MatrixXi F,G;
  if(argc <= 2)
    ./meshfix-libigl [input](.obj|.ply|.stl|.off) [output](.obj|.ply|.stl|.off)
    return EXIT_FAILURE;
  // Write to OBJ

Then you can run this by issuing something like:

./meshfix-libigl input.obj output.ply