## Archive for the ‘code’ Category

### Eigen performance gotcha calling non-templated function from templated one

Tuesday, July 25th, 2017

I just spent a while tracking down a rather surprising performance bug in my code.

Here’s a minimal example:

#include <Eigen/Dense>
#include <iostream>

int simple_size(const Eigen::MatrixXi & Z)
{
return Z.size();
}

template <typename T> int templated_size(const Eigen::MatrixBase<T> & Y)
{
return simple_size(Y);
}

int main(int argc, const char * argv[])
{
const int s = 40000;
Eigen::MatrixXi X = Eigen::MatrixXi::Zero(40000,40000);
std::cout<<"Compare:"<<std::endl;
std::cout<<(X.size()         ?"done":"")<<std::endl;
std::cout<<(simple_size(X)          ?"done":"")<<std::endl;
std::cout<<(templated_size(X)?"done":"")<<std::endl;

}

Running this, it will show that the last call to templated_size is taking way too long. Inspection will show that a copy of Y is being created to create a Eigen::MatrixXi & reference.

Now, clearly it’s poor design to call a function expecting a Eigen::MatrixXi & reference with a generic templated type Eigen::MatrixBase<T> &, but unfortunately this happens quite often with legacy libigl functions. My expectation was that since T is Eigen::MatrixXi in this case a simple reference would be passed.

It’s worth noting that const is actually creating/hiding the problem. Because simple_size takes a const reference, the compiler is happy to construct a Eigen::MatrixXi on the fly to create a valid reference. Without the consts the compiler stops at an error.

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

### Project page for “Generalized Matryoshka: Computational Design of Nesting Objects”

Wednesday, June 14th, 2017

This April I had fun working on a little project that’s been tickling my mind for a while. Can you make any shape into a Matryoshka doll?

I’ll be presenting my paper on this at SGP 2017. It’s entitled
Generalized Matryoshka: Computational Design of Nesting Objects.

### Pause (and then resume) Battery-Guzzling programs

Wednesday, June 7th, 2017

My laptop battery dies quickly these days. Certain apps (cough, cough, Slack), have very high idle CPU-usage. You can pause these programs with

killall -STOP Slack

And later you can resume the application with

killall -CONT Slack

### Convincing LatexIt and Illustrator to use the new SIGGRAPH fonts

Saturday, May 20th, 2017

The SIGGRAPH Latex style changed to the Libertine font. Here’re the steps to convince Latexit to use the new stylesheet and then to convince Illustrator to use the libertine font for drag and drop math.

mkdir ~/Library/texmf/tex/latex/local/acmart.cls/
cp ~/Dropbox/boundary/Paper/acmart.cls ~/Library/texmf/tex/latex/local/acmart.cls

In Latexit, open up Preferences, add a new SIGGRAPH “Template” containing:

\documentclass[sigconf, review]{acmart}
\pagenumbering{gobble}

If you try to drag and drop these into illustrator you’ll see that illustrator has replaced the nice math font with Myriad or something silly.

Drag this into FontBook.app

Update: I also had to issue:

If you see boxes with X’s replacing symbols after dragging and dropping from LaTeXit, then drag into Finder instead (to create a .pdf file), then open this directly and Illustrator will give a warning and tell you which font it’s (still) missing.

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

### Using glfw background window in matlab mex thread

Monday, April 3rd, 2017

A while ago, I tried to get glfw playing nicely with matlab’s mex files. I didn’t quite succeed.

My current project requires “GP-GPU” programming. I’m currently using glfw’s “background window” feature to create an OpenGL context. My first tests show that matlab’s mex will play nicely if glfwTerminate() is never called:

#include <igl/opengl/glfw/background_window.h>
#include <mex.h>
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
GLFWwindow * window;
igl::opengl::glfw::background_window(window);
glfwDestroyWindow(window);
//glfwTerminate();
}

You can compile this (on mac os x) with:

mex( ...
'-largeArrayDims','-DMEX','CXXFLAGS=$CXXFLAGS -std=c++11', ... 'LDFLAGS=\$LDFLAGS -framework Foundation -framework AppKit -framework Accelerate -framework OpenGL -framework AppKit -framework Carbon -framework QuartzCore -framework IOKit ', ...
'-I/usr/local/igl/libigl//external/nanogui/ext/eigen/', ...
'-I/usr/local/igl/libigl/include', ...
'-I/usr/local/igl/libigl/external/nanogui/ext/glfw/include/', ...
'-L/usr/local/igl/libigl/lib/','-lglfw3', ...
'background_glfw.cpp');

I don’t know why including glfwTerminate() causes matlab to sometimes crash. The error report is impossible to read and seems to have something to do with Quartz.

### Offset surface of triangle mesh in matlab

Wednesday, March 22nd, 2017

Here’s a little demonstration of how to use gptoolbox and MATLAB to generate an offset surfaces of a triangle mesh. This takes a mesh in V,F and creates a mesh SV,SF of the isosurface at signed distance iso:

% Extract offset at minus 3% of bounind box diagonal length
iso = -0.03;
% Resolution grid → resolution of extracted offset surface
side = 60;
% Amount of smoothing to apply to distance field
sigma = 1.4;
bbd = norm(max(V)-min(V));
% Pad generously for positive iso values
[BC,side,r] = voxel_grid([V;max(V)+iso*1;min(V)-iso*1],side);
D = signed_distance(BC,V,F);
D = reshape(D,side([2 1 3]));
% Smooth signed distance field
D = imfilter(D,fspecial('gaussian',9,sigma),'replicate');
BC3 = reshape(BC,[side([2 1 3]) 3]);
% Use matlab's built-in marching cubes iso-surface mesher (works most of the time)
surf = isosurface(BC3(:,:,:,1),BC3(:,:,:,2),BC3(:,:,:,3),D,iso*bbd);
SV = surf.vertices;
SF = surf.faces;

Here’s a blue bunny with a positive offset surface, an orange “cage”:

Here’s a blue bunny with a negative offset surface. This is useful for hollowing out objects to 3d print:

Because the iso-surface extraction will over tesselate low curvature patches of the output surface, it would make a lot of sense to remesh/decimate this mesh.

(to create these fancy renderings:)

clf;
hold on;
t  = tsurf(F,V,'EdgeColor','none',fsoft,  'FaceVertexCData',repmat(blue,size(V,1),1),'FaceAlpha',1+(iso<0)*(0.35-1),fphong);
ts = tsurf(SF,SV,'EdgeAlpha',0.2+(iso<0)*(0-0.2),fsoft,'FaceVertexCData',repmat(orange,size(SV,1),1),fphong,'FaceAlpha',1+(iso>0)*(0.2-1));
apply_ambient_occlusion(ts);
hold off;
axis equal;
view(-20,20)
camlight;
t.SpecularStrength = 0.04;
l = light('Position',[5 -5 10],'Style','local');
set(gca,'pos',[0 0 1 1])
set(gca,'Visible','off');
set(gcf,'Color','w');
drawnow;

### Quad meshing in matlab using gptoolbox

Friday, February 17th, 2017

This morning I played around with a very simple way of quad meshing a triangle mesh. This is more or less the “ideal pipeline” for quad meshing via parameterization. It will likely only work for very simple meshes. I’m pleasantly surprised that it works at all.

% Load a mesh with boundary and no topological (donut) holes
V = V*axisangle2matrix([1 0 0],-pi/2);
% Construct the LSCM matrix
[~,Q] = lscm(V,F,[],[]);
M = massmatrix(V,F);
% Use Fiedler vector as parameterization
[EV,ED] = eigs(Q,repdiag(M,2),5,'sm');
U = reshape(EV(:,end-3),[],2);
% try to find "canonical" rotation
[sU,sS,sV] = svd(U'*U);
U = U*sU;
% Lay down a grid in 2D
h = ((max(U(:,1))-min(U(:,1)))/(160-2));
[X,Y] = meshgrid(min(U(:,1))-h:h:max(U(:,1))+h,min(U(:,2))-h:h:max(U(:,2))+h);
% Find all quads that are strictly inside the parameterized mesh
[QQ,QU] = surf2patch(X,Y,0*Y);QU = QU(:,1:2);
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
QQ = QQ(all(abs(sqrD(QQ))<1e-10,2),:);
[QU,I] = remove_unreferenced(QU,QQ);
QQ = fliplr(I(QQ));
% Snap boundary vertices of quad mesh to boundary of mesh, and
% solve little Dirichlet problem to diffuse the displacements smoothly
QF = [QQ(:,[1 2 3]);QQ(:,[1 3 4])];
L = cotmatrix(QU,QF);
b = unique(outline(QQ));
[~,~,QU(b,:)] = point_mesh_squared_distance(QU(b,:),U,outline(F));
% Locate each quad mesh vertex in the parameterized mesh
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
B = barycentric_coordinates(C(:,1:2),U(F(I,1),:),U(F(I,2),:),U(F(I,3),:));
% Interpolate 3D positions at quad mesh vertex locations
QV = V(F(I,1),:).*B(:,1) + V(F(I,2),:).*B(:,2) + V(F(I,3),:).*B(:,3);

To create the image above use:

clf;
hold on;
tq = tsurf(QQ,QV,'FaceVertexCData',repmat(blue,size(QV,1),1),'EdgeColor','none',fsoft,fphong);
tt = tsurf(F,V-[0.8 0 0],'FaceVertexCData',repmat(1-0.8*(1-blue),size(V,1),1),'EdgeAlpha',0.5,fsoft,fphong);
axis equal;
camproj('persp');
view(-38,17);
drawn;
camlight;
plot_edges(QV,[QQ(:,2:3);QQ(:,[4 1])],'Color',orange);
plot_edges(QV,[QQ(:,1:2);QQ(:,3:4)],'w');hold off;
l = light('Position',[-10 0 10],'Style','infinite');
set(gca,'Visible','off');
set(gcf,'Color','w');

There’s a lot to improve about the quad mesh: the boundaries are not aligned with the meshing direction, the meshing directions are not always aligned with principle curvature, etc. Not bad for 30 mins of coding using tools that are just lying around.

### Using git and github for course assignments

Tuesday, January 10th, 2017

Here’s my tested plan for using git and github to manage student assignments for my upcoming course.

## Initial setup as the instructor:

We’re going to create a student repo and an instructor fork for each assignment. The motivation for separating assignments is to ensure that the diffs that students are submitting for their solutions are tidy. In my case, there’s not much shared code between assignments

Create public* repo that students will see: test-gp-hw0-student.

If you try to fork this on Github, it will tell you that you already have a fork. So instead trick github into forking it via the github importer, import test-gp-hw0-student as a private test-gp-hw0-instructor. Now you have two repos. We’ll set up a branch on test-gp-hw0-instructor to track test-gp-hw0-student for easy merging:

INSTRUCTOR_REPO=test-gp-hw0-instructor
STUDENT_REPO=test-gp-hw0-student
git clone https://github.com/alecjacobson/$INSTRUCTOR_REPO.git cd$INSTRUCTOR_REPO
git remote add student https://github.com/alecjacobson/$STUDENT_REPO.git git fetch student git checkout -b student-tracker --track student/master git checkout master cd .. To test this out, let’s create a file visible to everyone (students and instructors) in the student repo: git clone https://github.com/alecjacobson/$STUDENT_REPO.git
cd $STUDENT_REPO echo "visible to everyone" >> hello-everyone.md git add hello-everyone.md git commit -m "add file visible to everyone" hello-everyone.md git push cd ../ Now, let’s merge change into the instructor repo cd$INSTRUCTOR_REPO
git checkout master
git pull student master
cat hello-everyone.md
cd ..

Just to be thorough, let’s try the reverse: creating a file just visible to instructors:

cd stewartdent
git remote add instructor https://stewartdent@github.com/alecjacobson/$INSTRUCTOR_REPO.git git fetch instructor You should see an error: remote: Repository not found. fatal: repository 'https://stewartdent@github.com/alecjacobson/test-gp-hw0-instructor.git/' not found Now, let’s act as the student to submit the homework # cd stewartdent echo "stewartdent was here" >> assignment-01.md git add assignment-01.md git commit -m "Stewart Dent's Assignment 01 submission" assignment-01.md git push cd .. The student has pushed their submission onto their own repo. They can continue to work on it using git in all its glory. When they’re down they “hand in” their homework via a pull request onto the public test-gp-hw0-student repo. The student does this by visiting the Github page for their own fork of test-gp-hw0-student and clicking “New Pull Request”. ## Grading as the instructor The instructor can now grade this. I’m imagining that the test-gp-hw0-instructor repo has a solution that can be compared (in some way or another) to the students work. To pull the students work into the instructor copy (but never merge): cd$INSTRUCTOR_REPO
git checkout -b stewartdent-master master
git pull https://github.com/stewartdent/test-gp-hw0-student.git master
cat empty-form.md