A Simple Method for Correcting Facet Orientations in Polygon Meshes Based on Ray Casting

November 23rd, 2014

correct facet orientations on a bike, cell phone, chair, and truck

We’ve finally published our paper A Simple Method for Correcting Facet Orientations in Polygon Meshes Based on Ray Casting in the Journal of Computer Graphics Tools. The paper was written by Kenshi Takayama, Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung.

Abstract: We present a method for fixing incorrect orientations of facets in an input polygon mesh, a problem often seen in popular 3D model repositories, such that the front side of facets is visible from viewpoints outside of a solid shape represented or implied by the mesh. As opposed to previously proposed methods which are rather complex and hard to reproduce, our method is very simple, only requiring sampling visibilities by shooting many rays. We also propose a simple heuristic to handle interior facets that are invisible from exterior viewpoints. Our method is evaluated extensively with the SHREC Generic 3D Warehouse dataset containing 3168 manually designed meshes, and is demonstrated to be very effective.

You can find an implementation of this method in the embree extension of libigl: igl/embree/reorient_facets_raycast.h. If you store your mesh by its vertices in rows of a #V by 3 matrix V and triangle indices in the rows of a #F by 3 matrix F. Then you can quickly reorient your triangles to all point outward consistently with a single function call:

#include <igl/embree/reorient_facets_raycast.h>
...
Eigen::MatrixXI FF; // #F by 3 list of output triangle indices, some rows potentially reversed
Eigen::VectorXi I; // #F list of booleans revealing whether facet was reversed
igl::reorient_facets_raycast(V,F,FF,I);

As a preprocessor to our generalized winding numbers, this forms a powerful tool for determining the inside from outside for arbitrary meshes. This is especially important for creating volumetric tetrahedral meshes.

Tangible and modular input device on Swiss TV

November 17th, 2014

input device on tv

Our input device was featured on the Swiss television show “Gadgets” on the SRF channel. They’re using it to control a robot made out of a Ken Barbie doll.

Drop-in replacement header for switching from MathJax to KaTeX in markdown

November 7th, 2014

I wanted to try out switching to KaTeX from MathJax for our libigl tutorial. KaTeX needs a little bit of javascript besides the including the script to know to iterate over all class=math tags and render the latex code.

Also, it seems that markdown translates \\(x+y\\) into <span class=math>\(x+y\)</span> and the \( parenthesis throw off the KaTeX. You get an error like:

ParseError: KaTeX parse error: Expected 'EOF', got '\(' at position 2: \(x+y

So my header to switch from mathjax simply iterates over all math tags on document load using jquery and strips these parenthesis (which I expect on all tags):

<link rel="stylesheet" href="http://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.1.1/katex.min.css">
<script src="http://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.1.1/katex.min.js"></script>
<script src="http://code.jquery.com/jquery-1.11.0.min.js"></script>
<script src="http://code.jquery.com/jquery-migrate-1.2.1.min.js"></script>
<script>
$(function(){
  $(".math").each(function() {
    var texTxt = $(this).text();
    texTxt =  texTxt.slice(2,-2);
    el = $(this).get(0);
    if(el.tagName == "DIV"){
        addDisp = "\\displaystyle";
    } else {
        addDisp = "";
    }
    try {
        katex.render(addDisp+texTxt, el);
    }
    catch(err) {
        $(this).html("<span class='err'>"+err);
    }
  }); 
}); 
</script>

I didn’t invest much more into bulletproofing this because it seems like KaTeX is still in its infancy in terms of LaTeX support. It doesn’t support unicode characters (e.g. ∑, ∆, π) or common tags like \mathbf or \begin{cases}. Hopefully it will eventually, it is much faster than MathJax for the simple equations.

Using two arrays to store undirected edges instead of a map

November 7th, 2014

I profiled some of the libigl code the other day and found a performance leak corresponding to the way I sometimes dealt with undirected edges. In a triangle mesh, each oriented triangle has three directed edges. It’s sometimes useful to find other triangles with the same edges. So if my first triangle has an edge {i,j} I may want to find all other triangles who have either also an edge {i,j} or have the reverse edge {j,i}. The slow way of storing such a relationship is with a std::map or std::unordered_map. For example you might write,

std::map<std::pair<int,int>,std::vector<int>> edge_to_faces;
for each face f in F
{
  for each edge {i,j} in F(f)
  {
    if(i<j)
    {
      edge_to_faces[make_pair<int,int>(i,j)].push_back(f);
    }else
    {
      edge_to_faces[make_pair<int,int>(j,i)].push_back(f);
    }
  }
}

You can do slightly better using an std::unordered_map (aka hash map) and writing a compress(i,j) function which maps your {i,j} or {j,i} pair to a single number.

Turns out it’s much faster to build and to use two arrays. One which maps each directed edge instance to a unique undirected edge and then an array for the data stored at each undirected edge. Libigl provides a way of getting this first array via the igl::unique_edge_map function:

igl::unique_edge_map(F,E,uE,EMAP,uE2E);

Then the previous example looks like:

std::vector<int,std::vector<int>> uedge_to_faces;
for each face f in F
{
  for each cth edge {i,j} in F(f)
  {
    uedge_to_faces[EMAP(f+c*m)].push_back(f);
  }
}

Storing EMAP and uedge_to_faces is technically more memory than just edge_to_faces, but not asymptotically more if the data we’re storing for each undirected is O(#F), as is the case here. I think this boils down to memory access. The maps access is scattered but the array is contiguous. The price to build EMAP is a sort: O(m log m), but it can be down in place once and for all, where as the map needs to maintain a sort while building and conduct O(log m) searches when accessing.

Robust mesh boolean operations in libigl, gptoolbox

November 4th, 2014

I’ve added robust mesh boolean operations to libigl and a mex wrappers for matlab in gptoolbox. For comparison and as an alternative, I also included new wrappers cork’s boolean operations.

Check out the boolean entry in the libigl tutorial.

mesh booleans on cheburashka and knight

Linker order error with matlab dynamic libraries and gcc on mac os x

November 3rd, 2014

I had a nasty time tracking down a runtime error in our libigl tutorials. It seems that the newest version of matlab comes with some dynamic libraries which do not agree with the stdc++ library of my gcc compiler. I had a link command that looked something like this:

/opt/local/bin/g++ -std=c++11
  /Applications/MATLAB_R2014b.app/bin/maci64/libmex.dylib \
  /Applications/MATLAB_R2014b.app/bin/maci64/libmx.dylib \
  ...

Here’s the debugger output of the error I got at runtime

* thread #1: tid = 0x153aa3, 0x00007fff86d08866 libsystem_kernel.dylib`__pthread_kill + 10, queue = 'com.apple.main-thread', stop reason = signal SIGABRT
  * frame #0: 0x00007fff86d08866 libsystem_kernel.dylib`__pthread_kill + 10
    frame #1: 0x00007fff8341f35c libsystem_pthread.dylib`pthread_kill + 92
    frame #2: 0x00007fff8571fb1a libsystem_c.dylib`abort + 125
    frame #3: 0x00007fff8b1c007f libsystem_malloc.dylib`free + 411
    frame #4: 0x00000001001d47a0 libmex.dylib`std::basic_stringbuf<char, std::char_traits<char>, std::allocator<char> >::~basic_stringbuf() + 96
    frame #5: 0x0000000100451da2 libstdc++.6.dylib`std::basic_ostringstream<char, std::char_traits<char>, std::allocator<char> >::~basic_ostringstream() + 36

Seems like libmex and libmx were the culprits. I fixed this problem by telling the linker to find the stdc++ library before the matlab dynamic libraries:

/opt/local/bin/g++ -std=c++11
  -lstdc++ \
  /Applications/MATLAB_R2014b.app/bin/maci64/libmex.dylib \
  /Applications/MATLAB_R2014b.app/bin/maci64/libmx.dylib \
  ...

Triangle mesh for image as surface

October 26th, 2014

Here’s a little chuck of matlab code I use to make a surface mesh out of an image:

% load image as grayscale
im = rgb2gray(imresize(im2double(imread('hans-hass.jpg')),0.75));
% create triangle mesh grid
[V,F] = create_regular_grid(size(im,2),size(im,1),0,0);
% Scale (X,Y) to fit image
V(:,1) = V(:,1)*size(im,2);
V(:,2) = (1-V(:,2))*size(im,1);
V(:,3) = im(:);

A very simple method for approximate blue noise sampling on a triangle mesh

October 25th, 2014

This morning I was looking over the recent methods for generating blue-noise random samples on an arbitrary surface. Many of them involve dart throwing with complicated rejection schemes to achieve a Poisson-disk sampling (e.g. “Feature-preserving direct blue noise sampling for surface meshes” [Peyrot et al. 2014]), some involve hefty optimization procedures (e.g. “Blue Noise through Optimal Transport” [de Goes et al. 12]). A lot of the effort is spent trying to preserve features. I wasn’t particularly interested in this at the moment, and I just wanted something simple I could quickly code up and play with.

My idea stemmed from the observation that if my surface is a sphere, and I have a enough samples then the distance between samples will approximate their distance in Euclidean space. So an algorithm for Poisson-disk sampling the sphere looks like:

  1. Generate n uniformly random samples on the sphere
  2. For each sample i find distance dij to closest other sample j, in usual R3 sense
  3. Repeat steps 1 and 2 recursively for all samples i such that i < j and dij > re.

where re is the expected radius between samples on the surface. On the unit sphere, this is easy to compute: something like re = sqrt(4/n).

This works quite well:

blue noise sampling on sphere

on the left is the original “white noise” uniform sampling, and on the right I show the iterations of this algorithm converging to a Poisson-disk-like sampling.

Now, for an arbitrary surface my original assumption is no longer true. There may be parts very close to each other in R3 but very far away in terms of the geodesic distance on the surface.

However, there exist a family of methods that lift a surface to an (usually higher dimensional) embedding Rn so that Euclidean distance in Rn approximates geodesic distance (e.g. “Biharmonic Distance” [Lipman et al. 2010]). Using this I can alter my method only slightly to work with arbitrary surfaces:

  1. Lift surface to Rn using “Biharmonic embedding”
  2. Generate n uniformly random samples on the original surface (the one still in R3)
  3. For each sample i find distance dij to closest other sample j, after lifting to Rn
  4. Repeat steps 2 and 3 recursively for all samples i such that i < j and dij > re,

but this time re is not so obviously computed. It would be easy to compute the expected geodesic distance radius, but unfortunately the biharmonic embedding I’m using does not preserve units (i.e. dij is only relative to the geodesic distance between samples i and j on the mesh). I could probably compute the expected distortion during the lift to Rn and derive re from that. But I found it much easier to just take the average closest distance between the original n samples after being lifted to Rn. If I recall correctly, this should converge to the correct expected radius as n goes to infinity. So for large n this should be a reasonable thing to do.

It’s working quite well:

blue noise sampling on cheburashka

again, on the left is the original “white noise” uniform sampling, and on the right I show the iterations of this algorithm converging to a Poisson-disk-like sampling.

I’ve put up matlab code for this in the mesh/random_points_on_mesh.m in my gptoolbox on git. Here’s a use example:

[V,F] = load_mesh('/usr/local/igl/libigl/examples/shared/cheburashka.off');
n = 3000;
Pwhite = random_points_on_mesh(V,F,n);
Pblue = random_points_on_mesh(V,F,n,'Color','blue');
t = tsurf([F;size(V,1)+F],[V;V(:,1)+1.1*(max(V(:,1))-min(V(:,1))) V(:,2:3)],'FaceColor','w','EdgeColor','none','FaceAlpha',0.8);
apply_ambient_occlusion(t);
hold on;
scatter3(Pblue(:,1),Pblue(:,2),Pblue(:,3));
scatter3(Pwhite(:,1)+1.1*(max(V(:,1))-min(V(:,1))),Pwhite(:,2),Pwhite(:,3));
hold off;
axis equal;

Note: The current code is rebuilding a k-nearest neighbor acceleration data structure for all samples every iteration. Surely, this could be optimized.

The Freedom Tower is a coarse triangle mesh

October 19th, 2014

Let’s apply Loop subdivision to make it smoother.

freedom tower loop subdivision

Extruding a Bezier curve into a triangle mesh in maya

October 17th, 2014

Today I struggled to convince Maya to let me extrude a Bezier Curve into a solid shape (sweep a closed curve and finish with planar end caps). I could used the Surface > Extrude tool to extrude the curve and then select the boundary edges and use the Surface > Planar tool to close the endcaps, but this just creates a group of 3 surfaces which are not topologically connected.

My end goal today was to create something to send to the 3D printer. So in this case I eventually wanted a triangle mesh. Here’re the steps I took to convert a bezier curve to a polygonal mesh:

  1. draw bezier curves
    maya bezier curve to polygon mesh
  2. find the Polygons > Plane tool
    maya bezier curve to polygon mesh
  3. draw a plane behind the curves
    maya bezier curve to polygon mesh
  4. Select one of the curves and the plane
    maya bezier curve to polygon mesh
  5. Choose Edit Mesh > Project curve onto mesh
    maya bezier curve to polygon mesh
  6. Select the new projected curve and the plane
    maya bezier curve to polygon mesh
  7. Choose Edit Mesh > Split Mesh with projected curve
    maya bezier curve to polygon mesh
  8. Right click and hold and drag to select “Face” selection mode
    maya bezier curve to polygon mesh
  9. Mouse over the plane until just the filled curve appears (there seem to be many overlapping faces.
    maya bezier curve to polygon mesh
  10. Choose Edit > Invert Selection
    maya bezier curve to polygon mesh
  11. Then choose Edit > Delete
    maya bezier curve to polygon mesh
  12. Select just the original Bezier curve and delete it.
    maya bezier curve to polygon mesh
  13. Repeat steps 2-12 for the other curves (why can’t we do all curves at once?)
    maya bezier curve to polygon mesh
  14. Select both filled curves and choose the Polygons > Extrude tool
    maya bezier curve to polygon mesh
  15. Pull down on the widget’s arrow to extrude.
    maya bezier curve to polygon mesh
  16. Select both extruded surfaces and choose Mesh > Cleanup...
    maya bezier curve to polygon mesh
  17. Make sure 4-sided faces is unchecked
    maya bezier curve to polygon mesh
  18. Make sure Faces with zero geometry area is checked with very small Area tolerance (e.g. 0.00001)
    maya bezier curve to polygon mesh
  19. Hit Cleanup
    maya bezier curve to polygon mesh
  20. The choose Mesh > Triangulate
    maya bezier curve to polygon mesh
  21. The surface is now triangulated. Select everything.
    maya bezier curve to polygon mesh
  22. File > Export Selection and save as an obj
    maya bezier curve to polygon mesh

Wow. So 21 steps. Not particularly easy for a task I thought would be dead simple. I must be missing some faster way to do this.