Archive for June, 2015

Physical simulation “skinning” with biharmonic coordinates

Thursday, June 25th, 2015

This is a cute little example we didn’t have time to put into our upcoming SIGGRAPH paper Linear Subspace Design for Real-Time Shape Deformation.

octopus upsampling

In this example I first decimate the high resolution (blue) octopus to create the low-res orange octopus. I’m careful to use a decimate that keeps vertices at (or at least near) their original positions. Then I compute “biharmonic coordinates” for vertices of the blue mesh with respect to the vertices of the coarse orange mesh. Notice that the orange mesh does not need to be an enclosing cage. Rather its points just need to lie in or on the blue octopus (there’re even ways to deal with points outside the blue octopus, but right now I’m skipping that). These weights are computed in gptoolbox with:

b = knnsearch(orange_V,blue_tet_V);
W = biharmonic_coordinates(blue_tet_V,blue_tet_T,b);

Then I run a sort of co-rotational linear-elasticity simulation on the orange mesh with collisions activated for the floor. The blue mesh tracks the orange mesh by simply interpolating the orange vertex positions via the coordinates:

blue_V = W * orange_V;

Simple as that. Smooth physics “skinning” with a matrix-matrix multiplication. This is due to the affine precision of our coordinates and their squared Laplacian energy minimization. The “standard” approach to physics “skinning” would be to create the orange tet-mesh so that it envelops the blue mesh and then use piecewise-linear (i.e. non-smooth) interpolation to move the blue mesh vertices.

Update: I should make it super clear that this is not a subspace (aka model reduction) simulation. The simulation of the orange mesh is not aware of the blue mesh at all. Sometimes subspace physics would be more appropriate, though other times pushing the subspace through all parts of the simulation pipeline is impractical/impossible. Hence this prototypical result.

For posterity here’s the entire example code:

% load hi-res mesh
[uV,uF] = load_mesh('~/Dropbox/models/octopus-med.ply');
% load low-res mesh
[~,CV,CF] = qslim(uV,uF,900,'QSlimFlags','-O 0');
[CV,CF] = meshfix(CV,CF);
% Find closest points
b = knnsearch(uV,CV);
CV = uV(b,:);

% Tet-mesh fine resolution mesh and compute weights
[dV,dT,dF] = tetgen(uV,uF,'Flags','-q2');
W = biharmonic_coordinates(dV,dT,b);
W = W(1:size(uV,1),:);

[VV,VT,VF] = tetgen(CV,CF,'Flags','-q100 -Y');
off = mean(VV);
UU = bsxfun(@plus,VV,[0 0 1]);
vel = zeros(size(UU));
data = [];
M = massmatrix(VV,VT);
% This normalization is very important
M = M/max(M(:));
dt = 1e-2;

t = tsurf(VF,UU, ...
  'FaceColor',[0.8 0.5 0.2], ...
%t.Visible = 'off';
l = light('Style','local','Position',[5 0 20]);
l2 = light('Style','infinite','Position',[-1 -1 2]);
hold on;
up = @(UU) bsxfun(@plus,[1 0 0],W*UU(1:numel(b),:));
u = tsurf(uF,up(UU),'EdgeColor','none', ...
  'FaceVertexCData',repmat([0.3 0.4 1.0],size(uV,1),1), ...
shadows = false;
if shadows
  e = -1e-4;
  QV = [-0.8,-0.5 e; -0.8  0.5 e; 1.5  0.5 e;1.5,-0.5 e];
  QQ = [1 2 3 4];
  trisurf(QQ,QV(:,1),QV(:,2),QV(:,3),'EdgeColor','none','FaceColor',[0.7 0.7 0.7]);
  st = add_shadow(t,l,'Ground',[0 0 -1 0]);
  su = add_shadow(u,l,'Ground',[0 0 -1 0]);
  st.FaceColor = [0.5 0.5 0.5];
  su.FaceColor = [0.5 0.5 0.5];
hold off;
axis equal;
axis([-0.8 1.8 -0.5 0.5 0 1.5]);
axis manual;
T = get(gca,'tightinset');
set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);

while true
  dur = 0;
  UU0 = UU;
  [UU,data] = arap(VV,VT,[],[], ...
    'Energy','elements', ...
    'V0',UU0, ...
    'Data',data, ...
    'Dynamic',M*repmat([0 0 -9.8],size(VV,1),1), ...
    'MaxIter',100, ...
    'TimeStep',dt, ...
  vel = UU-UU0;
  C = UU(:,3)<0;
  UU(C,3) = -UU(C,3);
  vel(C,3) = -vel(C,3);
  dur = dur+dt;
  t.Vertices = UU;
  u.Vertices = up(UU);

  if shadows
    light_pos = [l.Position strcmp(l.Style,'local')];
    ground = [0 0 -1 0];
    d = ground * light_pos';
    shadow_mat = d*eye(4) - light_pos'*ground;
    H = [t.Vertices ones(size(t.Vertices,1),1)]*shadow_mat';
    st.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
    H = [u.Vertices ones(size(u.Vertices,1),1)]*shadow_mat';
    su.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));


Hausdorff distance between two triangle meshes

Tuesday, June 23rd, 2015

hausdorff distance

Unless I’m missing something the Hausdorff distance between two triangle meshes will be the maximum over the maximum minimum distance from the vertices of mesh A to the surface of mesh B and the maximum minimum distance from the vertices of mesh B to the surface of mesh A. That is, if p on A and q on B determine the Hausdorff distance between A and B, then it follows that p and/or q is a vertex of A or B respectively. This means you can compute Hausdorff distances easily if you already have a function to compute distances of points (vertices) to a triangle mesh.

So, I’ve added Hausdorff distance calculation between two meshes to my gptoolbox as the function hausdorff. Above I’m drawing the pair of points that determine the Hausdorff distance between the Cheburashka and the knight. I’m also drawing the intersection curves of these two meshes. Here’s how I generate this image:

% compute Hausdorff distance and "determining" pair
[d,pair] = hausdorff(VA,FA,VB,FB);
% Mesh intersection, reindex as a single mesh
[SV,SF,~,J,IM] = selfintersect([VA;VB],[FA;size(VA,1)+FB]);
SF = IM(SF);
% Extract edges along intersection
allE = [SF(:,[2 3]); SF(:,[3 1]); SF(:,[1 2])];
sortallE = sort(allE,2);
[uE,~,IC] = unique(sortallE,'rows');
uE2F = sparse(IC(:),repmat(1:size(SF,1),1,3)',1);
BE = uE(any(uE2F(:,J<=size(FA,1)),2)&any(uE2F(:,J>size(FA,1)),2),:);

% Draw meshes, intersection lines and Hausdorff distance "determiner"
t = tsurf(SF,SV,'CData',1*(J<=size(FA,1)),'EdgeColor','none','FaceAlpha',0.6);
colormap([0.3 0.4 1.0;0.75 0.8 1.0]);axis equal;view(2);
t.SpecularStrength = 0;
t.DiffuseStrength = 0.7;
t.AmbientStrength = 0.3;
hold on;
plot_edges(SV,BE,'LineWidth',2,'Color',[0.7 0.65 0.55]);
p = tsurf([1 2 2],pair,'EdgeColor',[0.9 0.5 0.1],'LineWidth',6);
hold off;

l = light('Position',[1 -1.5 1.8],'Style','infinite');
ht = add_shadow(t,l);
hp = add_shadow(p,l,'Ground',[0 0 -1 min(SV(:,3))]);
ht.FaceColor = [0.9 0.9 0.9];
hp.FaceColor = ht.FaceColor;
hp.EdgeColor = hp.FaceColor;
hp.LineWidth = p.LineWidth;

Update: I was missing something. The “determiner” points might (at least) be lying on the edges. Consider the two different triangulations of a non-planar quad. For now, consider the libigl and gptoolbox versions of hausdorff distance to be broken. A hopeful patch would be to do in-plane subdivision once, so that there are vertices on all of the edge midpoints, but I’m hesitant to claim that this will fix the problem.

non planar quad

GLFW c++ mesh viewer from matlab with Core OpenGL

Monday, June 15th, 2015

Here’s a proof of concept mex function for matlab that passes a triangle mesh and renders it in a new window using glfw and native opengl:

// make sure the modern opengl headers are included before any others
#include <OpenGL/gl3.h>
#define __gl_h_

#include <igl/frustum.h>
#include <igl/read_triangle_mesh.h>
#include <igl/get_seconds.h>
#include <igl/matlab/mexStream.h>
#include <igl/matlab/parse_rhs.h>
#include <igl/matlab/mexErrMsgTxt.h>
#include <Eigen/Core>
#include <GLFW/glfw3.h>
#include <mex.h>

#include <string>
#include <chrono>
#include <thread>
#include <iostream>

std::string vertex_shader = R"(
#version 330 core
uniform mat4 proj;
uniform mat4 model;
in vec3 position;
void main()
  gl_Position = proj * model * vec4(position,1.);
std::string fragment_shader = R"(
#version 330 core
out vec3 color;
void main()
  color = vec3(0.8,0.4,0.0);

// width, height, shader id, vertex array object
int w=800,h=600;
double highdpi=1;
GLuint prog_id=0;
GLuint VAO;
// Mesh data: RowMajor is important to directly use in OpenGL
Eigen::Matrix< float,Eigen::Dynamic,3,Eigen::RowMajor> V;
Eigen::Matrix<GLuint,Eigen::Dynamic,3,Eigen::RowMajor> F;
void mexFunction(int nlhs, mxArray *plhs[], 
   int nrhs, const mxArray *prhs[])
//#define mexErrMsgTxt(X) std::cerr<<X<<std::endl;exit(EXIT_FAILURE);
//int main()
  using namespace std;
  igl::MexStream mout;        
  std::streambuf *outbuf = cout.rdbuf(&mout);
  // Init glfw and create window + OpenGL context
    mexErrMsgTxt("Could not initialize glfw");
  const auto & error = []
    (int error, const char* description){
  glfwWindowHint(GLFW_SAMPLES, 4);
  GLFWwindow* window = glfwCreateWindow(w, h, "OpenGL in matlab? ! ? !", NULL, NULL);
    mexErrMsgTxt("Could not create glfw window");


      int major, minor, rev;
      major = glfwGetWindowAttrib(window, GLFW_CONTEXT_VERSION_MAJOR);
      minor = glfwGetWindowAttrib(window, GLFW_CONTEXT_VERSION_MINOR);
      rev = glfwGetWindowAttrib(window, GLFW_CONTEXT_REVISION);
      printf("OpenGL version recieved: %d.%d.%d\n", major, minor, rev);
      printf("Supported OpenGL is %s\n", (const char*)glGetString(GL_VERSION));
      printf("Supported GLSL is %s\n", (const char*)glGetString(GL_SHADING_LANGUAGE_VERSION));

  const auto & reshape = [] (GLFWwindow* window, int w, int h)

    int width, height;
    glfwGetFramebufferSize(window, &width, &height);
    int width_window, height_window;
    glfwGetWindowSize(window, &width_window, &height_window);
    highdpi = width/width_window;

  // Compile each shader
  const auto & compile_shader = [](const GLint type,const char * str) -> GLuint
    GLuint id = glCreateShader(type);
    return id;
  GLuint vid = compile_shader(GL_VERTEX_SHADER,vertex_shader.c_str());
  GLuint fid = compile_shader(GL_FRAGMENT_SHADER,fragment_shader.c_str());
  // attach shaders and link
  prog_id = glCreateProgram();
  GLint status;
  glGetProgramiv(prog_id, GL_LINK_STATUS, &status);

  // Read input mesh from rhs
  const int dim = mxGetN(prhs[1]);
  igl::mexErrMsgTxt(dim == 3,
    "Mesh vertex list must be #V by 3 list of vertex positions");
  igl::mexErrMsgTxt(dim == mxGetN(prhs[0]),
    "Mesh \"face\" simplex size must equal dimension");

  V.rowwise() -= V.colwise().mean();
  V /= (V.colwise().maxCoeff()-V.colwise().minCoeff()).maxCoeff();
  V /= 1.2;

  // Generate and attach buffers to vertex array
  glGenVertexArrays(1, &VAO);
  GLuint VBO, EBO;
  glGenBuffers(1, &VBO);
  glGenBuffers(1, &EBO);
  glBindBuffer(GL_ARRAY_BUFFER, VBO);
  glBufferData(GL_ARRAY_BUFFER, sizeof(float)*V.size(),, GL_STATIC_DRAW);
  glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(GLuint)*F.size(),, GL_STATIC_DRAW);
  glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), (GLvoid*)0);
  glBindBuffer(GL_ARRAY_BUFFER, 0); 

  // Main display routine
  while (!glfwWindowShouldClose(window))
      double tic = igl::get_seconds();
      // clear screen and set viewport

      // Projection and modelview matrices
      Eigen::Matrix4f proj = Eigen::Matrix4f::Identity();
      float near = 0.01;
      float far = 100;
      float top = tan(35./360.*M_PI)*near;
      float right = top * (double)::w/(double)::h;
      Eigen::Affine3f model = Eigen::Affine3f::Identity();
      // spin around
      static size_t count = 0;

      // select program and attach uniforms
      GLint proj_loc = glGetUniformLocation(prog_id,"proj");
      GLint model_loc = glGetUniformLocation(prog_id,"model");

      // Draw mesh as wireframe
      glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
      glDrawElements(GL_TRIANGLES, F.size(), GL_UNSIGNED_INT, 0);


        // In microseconds
        double duration = 1000000.*(igl::get_seconds()-tic);
        const double min_duration = 1000000./60.;
  // Restore the std stream buffer Important!

Compile with something like:

mex('LDOPTIMFLAGS=','CXX=clang++','LD=clang++','-DMEX -v -largeArrayDims','CXXFLAGS=$CXXFLAGS -std=c++11','-I/usr/local/igl/libigl/include','-I/usr/local/igl/libigl/external/glfw/include','-L/usr/local/igl/libigl/external/glfw/lib -lglfw3_clang','-I/opt/local/include/eigen3/','LDFLAGS=\$LDFLAGS -framework Foundation -framework AppKit -framework Accelerate -framework GLUT -framework OpenGL -framework IOKit -framework Cocoa -framework CoreVideo', 'tsurf_opengl.cpp');

And then if you issue:

[F,V] = surf2patch(surface(16*membrane(1,10)),'triangles');
V = V*axisangle2matrix([1 0 0],pi*0.5)*axisangle2matrix([0 1 0],-0.25*pi)*axisangle2matrix([1 0 0],-pi*0.1);

You should see something like:

opengl matlab glfw There are still quite a few kinks, but it’s at least a start.

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.

Improved distance queries in libigl and gptoolbox

Sunday, June 7th, 2015

tiny torus wire mesh, distance isosurface

I’ve stripped out the CGAL dependency from the libigl igl::point_mesh_squared_distance and igl::signed_distance functions. I’ve been meaning to do this for a while since it wasn’t really using anything that needs CGAL. Just conveniently wrapped around its axis-aligned bounding box hierarchy for a list of triangles.

Unfortunately CGAL tends to rely on callers to handle degeneracies, so zero-area triangles threw assertions. And it’s not particularly easy to use CGAL’s AABB hierarchy to hold a mixture of faces, edges and points.

Implementing a AABB-tree is really quite simple and libigl already had one for the igl::in_element function. I added distance query functions to this implementation and cleaned it up a bit.

The matlab wrappers (e.g. signed_distance_isosurface) in gptoolbox now use these distances. For example, the image above was created by using

 signed_distance_isosurface(V,E,'SignedDistanceType','unsigned', ...)

where E are the edges of the “Tiny Torus mesh” (grey with black wireframe). The blue surface is the watertight “inflated wire mesh”.