Physical simulation “skinning” with biharmonic coordinates

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);
VV=bsxfun(@minus,VV,off);
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], ...
  'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
%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), ...
  fphong,'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
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];
  set(gca,'Visible','off');
  set(gcf,'Color','w');
end
hold off;
axis equal;
axis([-0.8 1.8 -0.5 0.5 0 1.5]);
axis manual;
camproj('persp');
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, ...
    'Vm1',UU-vel);
  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));
  end

  drawnow;
end

Hausdorff distance between two triangle meshes

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;

drawnow; 
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;
apply_ambient_occlusion(t);
set(gcf,'Color','w');
set(gca,'Visible','off');
camproj('persp');
view(-27,6);

GLFW c++ mesh viewer from matlab with Core OpenGL

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>
#define GLFW_INCLUDE_GLU
#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
  if(!glfwInit())
  {
    mexErrMsgTxt("Could not initialize glfw");
  }
  const auto & error = []
    (int error, const char* description){
      mexErrMsgTxt(description);
    };
  glfwSetErrorCallback(error);
  glfwWindowHint(GLFW_SAMPLES, 4);
  glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
  glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 2);
  glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
  glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
  GLFWwindow* window = glfwCreateWindow(w, h, "OpenGL in matlab? ! ? !", NULL, NULL);
  if(!window)
  {
    glfwTerminate();
    mexErrMsgTxt("Could not create glfw window");
  }

  glfwMakeContextCurrent(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));

  glfwSetInputMode(window,GLFW_CURSOR,GLFW_CURSOR_NORMAL);
  const auto & reshape = [] (GLFWwindow* window, int w, int h)
  {
    ::w=w,::h=h;
  };
  glfwSetWindowSizeCallback(window,reshape);

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


  // Compile each shader
  const auto & compile_shader = [](const GLint type,const char * str) -> GLuint
  {
    GLuint id = glCreateShader(type);
    glShaderSource(id,1,&str,NULL);
    glCompileShader(id);
    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();
  glAttachShader(prog_id,vid);
  glAttachShader(prog_id,fid);
  glLinkProgram(prog_id);
  GLint status;
  glGetProgramiv(prog_id, GL_LINK_STATUS, &status);
  glDeleteShader(vid);
  glDeleteShader(fid);

  // 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");
  igl::parse_rhs_double(prhs+1,V);
  igl::parse_rhs_index(prhs+0,F);

  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);
  glBindVertexArray(VAO);
  glBindBuffer(GL_ARRAY_BUFFER, VBO);
  glBufferData(GL_ARRAY_BUFFER, sizeof(float)*V.size(), V.data(), GL_STATIC_DRAW);
  glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, EBO);
  glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(GLuint)*F.size(), F.data(), GL_STATIC_DRAW);
  glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(GLfloat), (GLvoid*)0);
  glEnableVertexAttribArray(0);
  glBindBuffer(GL_ARRAY_BUFFER, 0); 
  glBindVertexArray(0);

  // Main display routine
  while (!glfwWindowShouldClose(window))
  {
      double tic = igl::get_seconds();
      // clear screen and set viewport
      glClearColor(0.0,0.4,0.7,0.);
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
      glViewport(0,0,w,h);

      // 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;
      igl::frustum(-right,right,-top,top,near,far,proj);
      Eigen::Affine3f model = Eigen::Affine3f::Identity();
      model.translate(Eigen::Vector3f(0,0,-1.5));
      // spin around
      static size_t count = 0;
      model.rotate(Eigen::AngleAxisf(0.005*count++,Eigen::Vector3f(0,1,0)));

      // select program and attach uniforms
      glUseProgram(prog_id);
      GLint proj_loc = glGetUniformLocation(prog_id,"proj");
      glUniformMatrix4fv(proj_loc,1,GL_FALSE,proj.data());
      GLint model_loc = glGetUniformLocation(prog_id,"model");
      glUniformMatrix4fv(model_loc,1,GL_FALSE,model.matrix().data());

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

      glfwSwapBuffers(window);

      {
        glfwPollEvents();
        // In microseconds
        double duration = 1000000.*(igl::get_seconds()-tic);
        const double min_duration = 1000000./60.;
        if(duration<min_duration)
        {
          std::this_thread::sleep_for(std::chrono::microseconds((int)(min_duration-duration)));
        }
      }
  }
  glfwDestroyWindow(window);
  glfwTerminate();
  // Restore the std stream buffer Important!
  std::cout.rdbuf(outbuf);
}

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);
tsurf_opengl(F,V);

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

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

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

Code for web-based perceptual study on GitHub

May 20th, 2015

I’ve released the code for our web-based perceptual study as a github repo. It’s a combination of client-side javascript to serve up the survey questions and server-side php scripts to save the response to plain text files. This is a good set up for smallish surveys ~100/200 subjects. I found that it was much easier to work with organized plain text files than setting up an SQL database.

Hopefully it’s easy to see how to edit the php scripts/files/etc. to make use of this for your projects.

screenshot

Implementing “Linear Subspace Design for Real-Time Shape Deformation”

May 12th, 2015

Here’s a way to implement the linearly precise smoothness energy using gptoolbox in manner consistent with the paper.

First construct the usual cotangent laplacian L:

C = cotangent(V,F);
L = sparse( ...
  F(:,[2 3 1 3 1 2 2 3 1 3 1 2]), ...
  F(:,[3 1 2 2 3 1 2 3 1 3 1 2]), ...
  [C C -C -C], ...
  size(V,1),size(V,1));

This throws the cotangent of the angle across from each half-edge (i,j) to the off-diagonal entries L(i,j) and L(j,i) and minus that value to the diagonal entries L(i,i) and L(j,j).

You could also just use:

L = cotmatrix(V,F);

To compute the normal derivative matrix N, first find all of the half-edges on the boundary and the subscript indices into F so that if on_b_f(1) = f and on_b_c(1) = c then we know the half-edge across from vertex F(f,c) is on the boundary.

[~,on_b] = on_boundary(F);
[on_b_f,on_b_c] = find(on_b);

Now build lists of triplets (p,i,q) so that the half-edge i–>p is on the boundary and q is opposite it:

P = F(sub2ind(size(F),on_b_f,mod(on_b_c+1,3)+1));
I = F(sub2ind(size(F),on_b_f,mod(on_b_c,3)+1));
Q = F(sub2ind(size(F),on_b_f,on_b_c));

Collect the cotangents across from vertices i and j and throw accordingly at entries N(i,j), N(i,k), N(i,i), etc.:

CP = C(sub2ind(size(F),on_b_f,mod(on_b_c+1,3)+1));
CI = C(sub2ind(size(F),on_b_f,mod(on_b_c,3)+1));
N = sparse( ...
  [P P P P I I I I], ...
  [Q P Q I Q I Q P], ...
  [-CI CI -CP CP -CP CP -CI CI], ...
  size(V,1),size(V,1));

Now we can compute the “linearly precise Laplacian” K:

K = L+N;

Build the per-vertex mass matrix, and square our Laplacian to result in the bilaplacian smoothness energies quadratic form A:

M = massmatrix(V,F);
A = K' * (M\ K);

Accompanying video for “Linear Subspace Design for Real-Time Shape Deformation”

May 7th, 2015

Here’s the accompanying video for SIGGRAPH 2015 paper “Linear Subspace Design for Real-Time Shape Deformation”.

Preprint for “Linear Subspace Design for Real-Time Shape Deformation”

May 6th, 2015

linear subspace design thumbnail

Our upcoming SIGGRAPH 2015 paper “Linear Subspace Design for Real-Time Shape Deformation” is online now. This is joint work with Yu Wang, jernej Barbič and Ladislav Kavan. One way to interpret our contribution is a method for computing shape-aware (rather than cage-aware) generalized barycentric coordinates without a cage. Another way to interpret it is as a generalization of linear blend skinning to incorporate transformations of region regions/bones and translations at point handles. The distinguishing characteristic of our weight computation is its speed. The choice of handles (the choice of linear subspace) is now interactive and becomes part of the design process. It absorbs user-creativity just as manipulating the handles does.

Closed mesh of piece-wise constant height field surface from an image

May 6th, 2015

I pushed a little function box_height_field.m to gptoolbox. This function creates a height field from an image where each pixel becomes an extruded square (rather than just a point/vertex). There are two modes, one that’s fast, vectorized version which doesn’t add vertices so that the mesh is closed (though the underlying surface will still be “water-tight”). And a slower version which really creates a perfectly closed height field. Heres’ the result of the second one on the red-channel of the Hans Hass image:

im = im2double(imresize(rgb2gray(imread('hans-hass.jpg')),0.1));
[V,F] = box_height_field(im);

box height field

Here’s the same result but computed after quantizing the colors:

im = round(im/0.25)*0.25;
[V,F] = box_height_field(im);

box height field

You can clearly see the piecewise-constant regions. Using remesh_planar_patches you can reduce the number of version with losing any information:

[W,G] = remesh_planar_patches(V,F);

box height field