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

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
[~,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);
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.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);

light_pos = [l.Position strcmp(l.Style,'local')];
ground = [0 0 -1 0];
d = ground * light_pos';
st.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
su.Vertices = bsxfun(@rdivide,H(:,1:3),H(:,4));
end

drawnow;
end


### Hausdorff distance between two triangle meshes

Tuesday, June 23rd, 2015

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


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.

### 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/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 <iostream>

#version 330 core
uniform mat4 proj;
uniform mat4 model;
in vec3 position;
void main()
{
gl_Position = proj * model * vec4(position,1.);
}
)";
#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);
}

const auto & compile_shader = [](const GLint type,const char * str) -> GLuint
{
return id;
};
prog_id = glCreateProgram();
GLint 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");
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)
{
}
}
}
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:

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.

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

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