Posts Tagged ‘mesh’

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.

Mex wrapper for graph segmentation

Thursday, May 4th, 2017

I wrote a small mex wrapper for the graph segmentation part of the “Graph Based Image Segmentation” code. Most of the previously matlab implementations/wrappers worked on images. I want to apply this to geometry so I needed access to the graph segmentation directly. Here’s the wrapper (soon to be part of gptoolbox):

// mexopts = gptoolbox_mexopts('Static',false,'Debug',true);
// mex('segment_graph.cpp',mexopts{:});
#ifdef MEX
#  include <mex.h>
#  include <igl/C_STR.h>
#  include <igl/matlab/mexErrMsgTxt.h>
#  undef assert
#  define assert( isOK ) ( (isOK) ? (void)0 : (void) ::mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
#endif
#include "segment-graph.h"
#include <igl/matlab/mexErrMsgTxt.h>
#include <igl/matlab/parse_rhs.h>
#include <igl/unique.h>
#include <igl/matlab/prepare_lhs.h>
#include <igl/matlab/requires_arg.h>
#include <igl/matlab/validate_arg.h>
#include <igl/matlab/MexStream.h>
#include <Eigen/Sparse>
void mexFunction(
 int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  using namespace igl::matlab;
  using namespace Eigen;
  using namespace std;
  igl::matlab::MexStream mout;        
  std::streambuf *outbuf = std::cout.rdbuf(&mout);


  mexErrMsgTxt(nrhs>0,"Too few inputs");
  mexErrMsgTxt(mxIsSparse(prhs[0]),"Matrix should be sparse");
  const mxArray * mx_data = prhs[0];
  const int m = mxGetM(mx_data);
  const int n = mxGetN(mx_data);
  mexErrMsgTxt(n == mxGetM(prhs[0]), "Matrix should be square");
  assert(mxIsSparse(mx_data));
  assert(mxGetNumberOfDimensions(mx_data) == 2);
  // TODO: It should be possible to directly load the data into the sparse
  // matrix without going through the triplets
  // Copy data immediately
  double * pr = mxGetPr(mx_data);
  mwIndex * ir = mxGetIr(mx_data);
  mwIndex * jc = mxGetJc(mx_data);
  const int num_edges = mxGetNzmax(mx_data);
  edge * edges = new edge[num_edges];
  int k = 0;
  for(int j=0; j<n;j++)
  {
    // Iterate over inside
    while(k<(int)jc[j+1])
    {
      //cout<<ir[k]<<" "<<j<<" "<<pr[k]<<endl;
      assert((int)ir[k]<m);
      assert((int)j<n);
      edges[k].a = ir[k];
      edges[k].b = j;
      edges[k].w = pr[k];
      k++;
    }
  }

  // defaults 
  int min_size = 0;
  // Threshold 
  int c = sqrt((double)n);
  {
    int i = 1;
    while(i<nrhs)
    {
      mexErrMsgTxt(mxIsChar(prhs[i]),"Parameter names should be strings");
      // Cast to char
      const char * name = mxArrayToString(prhs[i]);
      if(strcmp("Threshold",name) == 0)
      {
        requires_arg(i,nrhs,name);
        validate_arg_scalar(i,nrhs,prhs,name);
        validate_arg_double(i,nrhs,prhs,name);
        c = (double)*mxGetPr(prhs[++i]);
      }else if(strcmp("MinSize",name) == 0)
      {
        requires_arg(i,nrhs,name);
        validate_arg_scalar(i,nrhs,prhs,name);
        validate_arg_double(i,nrhs,prhs,name);
        min_size = (int)((double)*mxGetPr(prhs[++i]));
      }
      i++;
    }
  }

  universe *u = segment_graph(n, num_edges, edges, c);

  // post process small components
  for (int i = 0; i < num_edges; i++) {
    int a = u->find(edges[i].a);
    int b = u->find(edges[i].b);
    if ((a != b) && ((u->size(a) < min_size) || (u->size(b) < min_size)))
      u->join(a, b);
  }

  switch(nlhs)
  {
    case 1:
    {
      plhs[0] = mxCreateDoubleMatrix(m,1, mxREAL);
      Eigen::VectorXi C(m);
      for(int i = 0;i<m;i++)
      {
        C(i) = u->find(i);
      }
      Eigen::VectorXi uC,I,J;
      igl::unique(C,uC,I,J);
      prepare_lhs_index(J,plhs);
    }
    default: break;
  }

  delete[] edges;
  delete u;
  std::cout.rdbuf(outbuf);
}

It takes the graph as a sparse matrix and outputs the component ids:

C = segment_graph(A);

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/{cow.off,cheburashka.off,decimated-knight.off}
#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++)
  {
    igl::read_triangle_mesh(argv[i],VV[i-1],FF[i-1]);
    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)
  {
    while(continue_computing)
    {
      // Let's watch this at 500x: Real clocks are slow.
      std::this_thread::sleep_for(std::chrono::milliseconds((int)tic*5));
      V *= Eigen::AngleAxisd(deg/180.*igl::PI,
          Eigen::Vector3d(0,0,1)).toRotationMatrix();
      {
        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;
  switch(VV.size())
  {
    case 3:
      threads.emplace_back(std::bind(rot,60.*60.,30,std::ref(VV[2])));
    case 2:
      threads.emplace_back(std::bind(rot,60.,6,std::ref(VV[1])));
    case 1:
      threads.emplace_back(std::bind(rot,1.,6,std::ref(VV[0])));
    default: break;
  }
  igl::viewer::Viewer v;
  v.data.clear();
  // 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;
    igl::combine(VV,FF,V,F);
    v.data.set_mesh(V,F);
  };
  set_meshes(VV,FF,v);
  // 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
    { 
      if(redraw) 
      {
        set_meshes(VV,FF,v); 
        {
          std::lock_guard<std::mutex> lock(m);
          redraw = false;
        }
      }
      return false;
    };
  v.launch();
  // 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.

Thingi10K: A Dataset of 10,000 3D-Printing Models

Tuesday, May 17th, 2016

thingi dataset

Qingnan “James” Zhou and I have released a technical report detailing the contents and methodology behind our ten thousand model thingi10k dataset.

Abstract
Empirically validating new 3D-printing related algorithms and implementations requires testing data representative of inputs encountered in the wild. An ideal benchmarking dataset should not only draw from the same distribution of shapes people print in terms of class (e.g., toys, mechanisms, jewelry), representation type (e.g., triangle soup meshes) and complexity (e.g., number of facets), but should also capture problems and artifacts endemic to 3D printing models (e.g., self-intersections, non-manifoldness). We observe that the contextual and geometric characteristics of 3D printing models differ significantly from those used for computer graphics applications, not to mention standard models (e.g., Stanford bunny, Armadillo, Fertility). We present a new dataset of 10,000 models collected from an online 3D printing model-sharing database. Via analysis of both geometric (e.g., triangle aspect ratios, manifoldness) and contextual (e.g., licenses, tags, classes) characteristics, we demonstrate that this dataset represents a more concise summary of real-world models used for 3D printing compared to existing datasets. To facilitate future research endeavors, we also present an online query interface to select subsets of the dataset according to project-specific characteristics. The complete dataset and per-model statistical data are freely available to the public.

Boolean operations using generalized winding numbers

Tuesday, February 2nd, 2016

booleans using generalized winding number

I posted a little tech report about how to use the generalized winding number for constructive-solid geometry (CSG) style boolean operations (union, intersection, difference etc.) on nasty meshes.

Code implementing this has existed for a little while now in the mesh_boolean_winding_number function of gptoolbox.

Unzip OBJ-style mesh into a per-vertex attribute mesh

Wednesday, January 6th, 2016

The .obj mesh file format allows corners of triangles to pull attributes from different sources. The triangle:

v 0 0 0
v 1 0 0
v 0 0 0
vn 1 0 0
vn 0 1 0
vt 0 0 0
f 1/1/1 2/1/1 3/1/2

pulls vertex positions from entries (1,2,3) in the v ... vertex position list, texture coordinates from entries (1,1,1) in the vt ... list, and normals from entries (1,2,2) in the vn normals list.

If we think of corners being described by all attributes there are potentially #F*3 distinct corners. Often information is shared, and in the best case the position/texture/normal indices are all the same, so #V distinct corners. Usually it’s some mixture in between.

I added igl::unzip_corners to libigl which “unzips” an OBJ-style mesh into per-vertex attribute mesh: each new vertex is a distinct corner, so the new face list indexes all attributes in lock step (ideal for OpenGL). I was careful to determine uniqueness combinatorially so, for example, combinatorially distinct input vertices happening to share the same attributes (two vertices at the same position, etc.) don’t get merged (you could always use igl::remove_duplicate_vertices if you wanted to do that).

Here’s a little demo program:

  Eigen::MatrixXd V, TC, N;
  Eigen::MatrixXi F,FTC,FN;
  igl::readOBJ(argv[1],V,TC,N,F,FTC,FN);
  if(FTC.size() == 0)
  {
    FTC = F;
  }
  Eigen::MatrixXi U,G,J;
  igl::unzip_corners<Eigen::MatrixXi>({F,FTC},U,G,J);
  // New mesh vertices and texture coordinates indexed by G
  GV = igl::slice(Eigen::MatrixXd(V),U.col(0),1);
  GTC = igl::slice(Eigen::MatrixXd(TC),U.col(1),1);

Write a triangle mesh to xml file in libigl using Eigen matrix templated on CGAL’s Exact Kernel

Thursday, December 17th, 2015

Lately I’ve been working with CGAL and its arbitrary precision kernel. These are convenient, but writing to standard floating point precision mesh file formats requires rounding. Here’s a relatively straightforward way of using libigl‘s existing xml serialization routines to write an Eigen matrix tempalated on the CGAL Exact kernels number type (e.g. a matrix of vertex positions) to an xml file (using ASCII or base64 binary encoding).

Notice that there’s a bit of gnarly template specialization that needs to be done inside the igl::xml::serialization_xml namespace. The structure of the libigl xml serialization code is unfortunately not very flexible.

#include <igl/xml/serialize_xml.h>
#include <Eigen/Core>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <iostream>

typedef CGAL::Epeck::FT EScalar;
typedef Eigen::Matrix<EScalar,Eigen::Dynamic,Eigen::Dynamic> MatrixXE;

namespace igl
{
  namespace xml
  {
    namespace serialization_xml
    {
      template <> inline void serialize(
        const MatrixXE & obj,
        tinyxml2::XMLDocument* doc,
        tinyxml2::XMLElement* element,
        const std::string& name)
      {
        const std::function<std::string(const MatrixXE::Scalar &) > to_string = 
          [](const MatrixXE::Scalar & v)->std::string
          {
            return
              STR(CGAL::exact(v));
          };
        serialize(obj,name,to_string,doc,element);
      }
      template <> inline void deserialize(
        MatrixXE & obj,
        const tinyxml2::XMLDocument* doc,
        const tinyxml2::XMLElement* element,
        const std::string& name)
      {
        const std::function<void(const std::string &,MatrixXE::Scalar &)> & 
          from_string = 
          [](const std::string & s, MatrixXE::Scalar & v)
          {
            std::stringstream(s)>>v;
          };
        deserialize(doc,element,name,from_string,obj);
      }
    }
  }
}

int main(int argc, char * argv[])
{
  using namespace std;
  using namespace Eigen;
  // Little 4-vertex, 2-face mesh with rational coordinates
  MatrixXE V(4,3),W;
  V<<
    EScalar(1)/EScalar(3),  EScalar(1)/EScalar(13), EScalar(1)/EScalar(7),
    EScalar(1)/EScalar(7),  EScalar(1)/EScalar(3),  EScalar(1)/EScalar(13),
    EScalar(1)/EScalar(13), EScalar(1)/EScalar(7),  EScalar(1)/EScalar(3),
    EScalar(1)/EScalar(13), EScalar(1)/EScalar(7), -EScalar(1)/EScalar(3);
  MatrixXi F(2,3),G;
  F<<0,1,2,0,2,3;

  // Write mesh
  const bool binary = false;
  // Write vertices, overwriting file (true)
  igl::xml::serialize_xml(V,"vertices","exact.xml",binary,true);
  // Write faces to same file, appending (false)
  igl::xml::serialize_xml(F,"faces","exact.xml",binary,false);

  // Read mesh
  igl::xml::deserialize_xml(W,"vertices","exact.xml");
  igl::xml::deserialize_xml(G,"faces","exact.xml");

  // Verify 
  cout<<"V "<<(V.isApprox(W,0) ? "equals" : "does not equal")<<" W"<<endl;
}

Here’s a little feeling for the overhead of this format on a big single-precision mesh cast to the exact kernel (the expressions are simple and short, so this is sort of a best case scenario):

|           | .ply|.xml ascii|.xml binary| .obj|
|-----------|-----|----------|-----------|-----|
|File size: |172MB|     511MB|      330MB|288MB|
|Zip size:  |  68M|     134MB|       74MB| 81MB|
|Read time: |   8s|      270s|         9s|  32s|
|Write time:|  13s|       46s|        13s|  28s|

Nested Cages project page

Friday, October 2nd, 2015

nested cages bunny

We’ve posted a project page for our upcoming SIGGRAPH Asia paper Nested Cages, a collaboration between Leonardo Sacht, Etienne Vouga and myself.

Abstract: Many tasks in geometry processing and physical simulation benefit from multiresolution hierarchies. One important characteristic across a variety of applications is that coarser layers strictly encage finer layers, nesting one another. Existing techniques such as surface mesh decimation, voxelization, or contouring distance level sets do not provide sufficient control over the quality of the output surfaces while maintaining strict nesting. We propose a solution that enables use of application-specific decimation and quality metrics. The method constructs each next-coarsest level of the hierarchy, using a sequence of decimation, flow, and contact-aware optimization steps. From coarse to fine, each layer then fully encages the next while retaining a snug fit. The method is applicable to a wide variety of shapes of complex geometry and topology. We demonstrate the effectiveness of our nested cages not only for multigrid solvers, but also for conservative collision detection, domain discretization for elastic simulation, and cage-based geometric modeling.

Mesh untangling in gptoolbox

Thursday, June 11th, 2015

I’ve reimplemented the first half of our paper, “Consistent Volumetric Discretizations Inside Self-Intersecting Surfaces” [Sacht et al. 2013] as the untangle function in my gptoolbox. Our paper expands on the idea of “Interference-Aware Geometric Modeling” [Harmon et al. 2011]: given a self-intersecting mesh, run a mean-curvature flow until all self-intersections disappear and then reverse the flow but activate collision detection and response. My new implementation uses the conformalized mean curvature flow of “Can Mean-Curvature Flow Be Made Non-Singular?” [Kazhdan et al. 2012] and then reverse the flow using el topo, “Robust Topological Operations for Dynamic Explicit Surfaces” [Brochu & Bridson 2009]. For the reverse steps, I’m simply filtering the reversed flow “velocities”, then once returned to it’s original state, I run some iterations of ARAP gradient descent (with collision detection and response via el topo) to restore the intrinsic shape of the mesh.

Here’s the result on an example with a rather drastic self-intersection.

hand ok shrink inflate, untangle

For simpler models with lots of little self-intersections (like an otherwise nice subdivision surface), this method works very well. Unlike meshfix, this will not delete or remesh any part of the model.

Writing a mesh to an .obj file in a single line with Eigen

Wednesday, April 29th, 2015

I was recently revisiting our igl::writeOBJ code and realized that for simple meshes, the code for writing an .obj file can be really simple:

ofstream("test.obj")<<
  V.format(IOFormat(FullPrecision,0," ","\n","v ","","","\n"))<<
  (F.array()+1).format(IOFormat(FullPrecision,0," ","\n","f ","","","\n"));

Update: And because why not, here’s a .off writer:

ofstream("test.off")<<
  "OFF\n"<<V.rows()<<" "<<F.rows()<<" 0\n"<<
  V.format(IOFormat(FullPrecision,0," ","\n","","","","\n"))<<
  (F.array()).format(IOFormat(FullPrecision,0," ","\n","3 ","","","\n"));