## Archive for May, 2015

### Code for web-based perceptual study on GitHub

Wednesday, 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.

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

Tuesday, 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”

Thursday, 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”

Wednesday, May 6th, 2015

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

Wednesday, 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);


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

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


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


### Depth peeling mini-app

Monday, May 4th, 2015

I’ve been lamenting poor transparency handling in opengl for as long as I’ve used opengl. From a graphics programmer’s perspective, it’s just so frustrating that you can’t just:

glEnable(GL_ALPHA_BLEND_CORRECTLY);


Finally I’ve gotten around to implementing a correct solution. The idea is simple, and, by now, well-known. Render the scene k times. For each pass keep track of the color and depth of the fragment with the smallest depth value which is less than the previous pass’s. In this way, each pass peels off a layer of pixels from front-to-back. Finally render all of these as composite image from back-to-front with the usual alpha blending.

If k is more than the maximum number of overlapping transparent objects then this will give the correct result. In practice if the alpha values aren’t all tiny (everything’s only a bit transparent). Then only a few passes are needed.

Here’s a little 280-line glut app demonstrating this:

// make sure the modern opengl headers are included before any others
#include <OpenGL/gl3.h>
#define __gl_h_
#include <igl/frustum.h>
#include <Eigen/Core>
#include <GLUT/glut.h>
#include <string>

void init_render_to_texture(
const size_t w, const size_t h, GLuint & tex, GLuint & dtex, GLuint & fbo)
{
const auto & gen_tex = [](GLuint & tex)
{
// http://www.opengl.org/wiki/Framebuffer_Object_Examples#Quick_example.2C_render_to_texture_.282D.29
glGenTextures(1, &tex);
glBindTexture(GL_TEXTURE_2D, tex);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
};
// Generate texture for colors and attached to color component of framebuffer
gen_tex(tex);
glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, w, h, 0, GL_BGRA, GL_FLOAT, NULL);
glBindTexture(GL_TEXTURE_2D, 0);
glGenFramebuffers(1, &fbo);
glBindFramebuffer(GL_FRAMEBUFFER, fbo);
// Generate texture for depth and attached to depth component of framebuffer
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, tex, 0);
gen_tex(dtex);
glTexImage2D(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT32, w, h, 0, GL_DEPTH_COMPONENT, GL_FLOAT, NULL);
glFramebufferTexture2D(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_TEXTURE_2D, dtex, 0);
// Clean up
glBindFramebuffer(GL_FRAMEBUFFER, 0);
glBindTexture(GL_TEXTURE_2D,0);
}

// For rendering a full-viewport quad, set tex-coord from position
#version 330 core
in vec3 position;
out vec2 tex_coord;
void main()
{
gl_Position = vec4(position,1.);
tex_coord = vec2(0.5*(position.x+1), 0.5*(position.y+1));
}
)";
// Render directly from color or depth texture
#version 330 core
in vec2 tex_coord;
out vec4 color;
uniform sampler2D color_texture;
uniform sampler2D depth_texture;
uniform bool show_depth;
void main()
{
vec4 depth = texture(depth_texture,tex_coord);
// Mask out background which is set to 1
if(depth.r<1)
{
color = texture(color_texture,tex_coord);
if(show_depth)
{
// Depth of background seems to be set to exactly 1.
color.rgb = vec3(1,1,1)*(1.-depth.r)/0.006125;
}
}else
{
}
}
)";

// Pass-through vertex shader with projection and model matrices
#version 330 core
uniform mat4 proj;
uniform mat4 model;
in vec3 position;
void main()
{
gl_Position = proj * model * vec4(position,1.);
}
)";
// Render if first pass or farther than closest frag on last pass
#version 330 core
out vec4 color;
uniform bool first_pass;
uniform float width;
uniform float height;
uniform sampler2D depth_texture;
void main()
{
color = vec4(0.8,0.4,0.0,0.75);
color.rgb *= (1.-gl_FragCoord.z)/0.006125;
if(!first_pass)
{
vec2 tex_coord = vec2(float(gl_FragCoord.x)/width,float(gl_FragCoord.y)/height);
float max_depth = texture(depth_texture,tex_coord).r;
if(gl_FragCoord.z <= max_depth)
{
}
}
}
)";

// shader id, vertex array object
GLuint scene_p_id=0,tex_p_id;
GLuint VAO,QVAO;
// Number of passes
#define k 4
GLuint tex_id[k],dtex_id[k],fbo_id[k];
// full width/height of window, width/height of viewports
int full_w=1440,full_h=480,w=full_w/(k+2),h=full_h/1;
// Mesh data: RowMajor is important to directly use in OpenGL
typedef Eigen::Matrix< float,Eigen::Dynamic,3,Eigen::RowMajor> MatrixV;
typedef Eigen::Matrix<GLuint,Eigen::Dynamic,3,Eigen::RowMajor> MatrixF;
MatrixV V,QV;
MatrixF F,QF;
int main(int argc, char * argv[])
{
// Init glut and create window + OpenGL context
glutInit(&argc,argv);
glutInitDisplayMode(GLUT_3_2_CORE_PROFILE|GLUT_RGBA|GLUT_DOUBLE|GLUT_DEPTH);
glutInitWindowSize(full_w,full_h);
glutCreateWindow("test");
// Prepare VAOs
const auto & vao = [](const MatrixV & V, const MatrixF & F, GLuint & VAO)
{
// 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);
};

// Read input mesh from file
V.rowwise() -= V.colwise().mean();
V /= (V.colwise().maxCoeff()-V.colwise().minCoeff()).maxCoeff();
vao(V,F,VAO);
// square
const MatrixV QV = (MatrixV(4,3)<<-1,-1,0,1,-1,0,1,1,0,-1,1,0).finished();
const MatrixF QF = (MatrixF(2,3)<< 0,1,2, 0,2,3).finished();
vao(QV,QF,QVAO);

// Main display routine
glutDisplayFunc(
[]()
{
// Projection and modelview matrices
Eigen::Matrix4f proj = Eigen::Matrix4f::Identity();
float near = 0.01;
float far = 3;
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(M_PI/180.*count++,Eigen::Vector3f(0,1,0)));

glEnable(GL_DEPTH_TEST);
glViewport(0,0,w,h);
// select program and attach uniforms
glUseProgram(scene_p_id);
GLint proj_loc = glGetUniformLocation(scene_p_id,"proj");
glUniformMatrix4fv(proj_loc,1,GL_FALSE,proj.data());
GLint model_loc = glGetUniformLocation(scene_p_id,"model");
glUniformMatrix4fv(model_loc,1,GL_FALSE,model.matrix().data());
glUniform1f(glGetUniformLocation(scene_p_id,"width"),w);
glUniform1f(glGetUniformLocation(scene_p_id,"height"),h);
glBindVertexArray(VAO);
glDisable(GL_BLEND);
for(int pass = 0;pass<k;pass++)
{
const bool first_pass = pass == 0;
glUniform1i(glGetUniformLocation(scene_p_id,"first_pass"),first_pass);
if(!first_pass)
{
glUniform1i(glGetUniformLocation(scene_p_id,"depth_texture"),0);
glActiveTexture(GL_TEXTURE0 + 0);
glBindTexture(GL_TEXTURE_2D, dtex_id[pass-1]);
}
glBindFramebuffer(GL_FRAMEBUFFER, fbo_id[pass]);
glClearColor(0.0,0.4,0.7,0.);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDrawElements(GL_TRIANGLES, F.size(), GL_UNSIGNED_INT, 0);
}
// clean up and set to render to screen
glBindVertexArray(0);
glBindFramebuffer(GL_FRAMEBUFFER, 0);
glActiveTexture(GL_TEXTURE0 + 0);
glBindTexture(GL_TEXTURE_2D,0);

glBindVertexArray(QVAO);
glClearColor(0.0,0.4,0.7,0.);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glUseProgram(tex_p_id);
// Draw result of each peel
for(int pass = 0;pass<k;pass++)
{
GLint color_tex_loc = glGetUniformLocation(tex_p_id,"color_texture");
glUniform1i(color_tex_loc, 0);
glActiveTexture(GL_TEXTURE0 + 0);
glBindTexture(GL_TEXTURE_2D, tex_id[pass]);
GLint depth_tex_loc = glGetUniformLocation(tex_p_id,"depth_texture");
glUniform1i(depth_tex_loc, 1);
glActiveTexture(GL_TEXTURE0 + 1);
glBindTexture(GL_TEXTURE_2D, dtex_id[pass]);
glViewport(pass*w,0*h,w,h);
glUniform1i(glGetUniformLocation(tex_p_id,"show_depth"),0);
glDrawElements(GL_TRIANGLES,6, GL_UNSIGNED_INT, 0);
}

// Render final result as composite of all textures
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glEnable(GL_BLEND);
glDepthFunc(GL_ALWAYS);
glViewport(k*w,0*h,w,h);
glUniform1i(glGetUniformLocation(tex_p_id,"show_depth"),0);
GLint color_tex_loc = glGetUniformLocation(tex_p_id,"color_texture");
GLint depth_tex_loc = glGetUniformLocation(tex_p_id,"depth_texture");
for(int pass = k-1;pass>=0;pass--)
{
glUniform1i(color_tex_loc, 0);
glActiveTexture(GL_TEXTURE0 + 0);
glBindTexture(GL_TEXTURE_2D, tex_id[pass]);
glUniform1i(depth_tex_loc, 1);
glActiveTexture(GL_TEXTURE0 + 1);
glBindTexture(GL_TEXTURE_2D, dtex_id[pass]);
glDrawElements(GL_TRIANGLES,6, GL_UNSIGNED_INT, 0);
}
glDepthFunc(GL_LESS);
// Render scene using naive GL_BLEND transparency
glUseProgram(scene_p_id);
glBindVertexArray(VAO);
glViewport((k+1)*w,0*h,w,h);
glDrawElements(GL_TRIANGLES, F.size(), GL_UNSIGNED_INT, 0);
glBindVertexArray(0);
glDisable(GL_BLEND);

glutSwapBuffers();
glutPostRedisplay();
}
);
glutReshapeFunc(
[](int w,int h)
{
full_h=h;
full_w=w;
::w=full_w/(k+2);
::h=full_h/(1);
// (re)-initialize textures and buffers
for(size_t i = 0;i<k;i++)
{
init_render_to_texture(::w,::h,tex_id[i],dtex_id[i],fbo_id[i]);
}
});
glutMainLoop();
}


Running this for an elephant mesh you get:

The left 4 images are the individual peeled layers, rendered only for demonstration. The next image is the composite (the main result) and for comparison the nasty mess you get if you just try to render your scene with:

glBlendFunc(GL_SRC, GL_ONE_MINUS_SRC_ALPHA);
glEnable(GL_BLEND);


**Update: ** Turns out to get the right blending you should be careful to turn off GL_BLEND when doing the peeling pass and only turn it on for the texture compositing.

### Tiny mesh viewer example using modern OpenGL core profile + Eigen + libigl + GLUT

Saturday, May 2nd, 2015

I was pleased to find out that the glut built in on recent macs does indeed support modern opengl. You can enable it via the initialization with:

glutInitDisplayMode(GLUT_3_2_CORE_PROFILE|GLUT_RGBA|GLUT_DOUBLE|GLUT_DEPTH);


If you’re using my patched glut you can also use the string version:

glutInitDisplayString("rgba depth double samples>=8 stencil core");


Together with libigl and Eigen (and made easy with C++11 features) you can write a little mesh viewer in modern opengl code in under 130 lines:

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

#include <igl/frustum.h>
#include <Eigen/Core>
#include <GLUT/glut.h>
#include <string>
#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;
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;
int main(int argc, char * argv[])
{
// Init glut and create window + OpenGL context
glutInit(&argc,argv);
//glutInitDisplayString("rgba depth double samples>=8 stencil core");
glutInitDisplayMode(GLUT_3_2_CORE_PROFILE|GLUT_RGBA|GLUT_DOUBLE|GLUT_DEPTH);
glutInitWindowSize(w,h);
glutCreateWindow("test");

const auto & compile_shader = [](const GLint type,const char * str) -> GLuint
{
return id;
};
prog_id = glCreateProgram();
GLint status;

// Read input mesh from file
V.rowwise() -= V.colwise().mean();
V /= (V.colwise().maxCoeff()-V.colwise().minCoeff()).maxCoeff();

// 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
glutDisplayFunc(
[]()
{
// 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.01*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);

glutSwapBuffers();
glutPostRedisplay();
}
);
glutReshapeFunc([](int w,int h){::w=w,::h=h;});
std::cout<<"OpenGL version: "<<glGetString(GL_VERSION)<<std::endl;
glutMainLoop();
}


This produces:

and on my mac outputs:

OpenGL version: 4.1 INTEL-10.6.20


### Generate list of k distant unique colors

Saturday, May 2nd, 2015

As far as I know, the best way to do picking in OpenGL is still to render each pickable object (pickle?) using a unique color. Here’s what I use to get k distant and unique RGB colors. The idea is very straightforward. Sample the RGB color cube at points on the k lattice. Here’s what that looks like for k = 103

This is easy to generate in matlab with:

k = 10*10*10;
s = ceil(k^(1/3));
[R,G,B] = ind2sub([s s s],1:k);
rgb = uint8([R;G;B]'-1)/(s-1)*255);


Here’s a plot of the colors in one image:

If you’re object ids have spacial coherency then this will still place similar colors near each other (shouldn’t really be a problem if you use a decent precision for your picking buffer). But you can always randomize the order with:

I = randperm(k); rgb = uint8([R(I);G(I);B(I)]’-1)/(s-1)*255);

to give you something like

Here’s a C++ version using Eigen and libigl:

  size_t k = 1024;
size_t s = ceil(cbrt(k));
Matrix<unsigned char,Dynamic,3> rgb(k,3);
ArrayXi I,re;
igl::randperm(k,I);
igl::mod(I,s*s,re);
rgb.col(2) = (((I-re)/(s*s)*255)/(s-1)).cast<unsigned char>();
I = re;
igl::mod(I,s,re);
rgb.col(1) = (((I-re)/(s)*255)/(s-1)).cast<unsigned char>();
rgb.col(0) = ((re*255)/(s-1)).cast<unsigned char>();


Or for the faint of heart:

  igl::randperm(k,I);
{
size_t i = 0;
for(size_t b = 0;b<s&&i<k;b++)
{
for(size_t g = 0;g<s&&i<k;g++)
{
for(size_t r = 0;r<s&&i<k;r++)
{
rgb(I(i),0) = r*255/(s-1);
rgb(I(i),1) = g*255/(s-1);
rgb(I(i),2) = b*255/(s-1);
i++;
}
}
}
}


Update: Here’s another way which requires less computation and doesn’t require knowing how many unique colors you want ahead of time (though it sort of assumes a smallish number). This won’t promise distant colors like above, but will still work for picking if you use unsigned byte precision:

[R,G,B] = ind2sub([256 256 256],mod((0:k-1)*(103811),256^3)+1);
rgb = uint8(permute(([R;G;B]'-1),[1 3 2]));


The idea is to sample the 2563 cube with a large prime stride. I choose 103811 because it produced a nice balance of colors for k==100.

;