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

Depth peeling mini-app

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 <igl/read_triangle_mesh.h>
#include <igl/create_shader_program.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
std::string tex_v_shader = R"(
#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
std::string tex_f_shader = R"(
#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
  {
    discard;
  }
}
)";

// Pass-through vertex shader with projection and model matrices
std::string scene_v_shader = R"(
#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
std::string scene_f_shader = R"(
#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)
    {
      discard;
    }
  }
}
)";

// 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");
  // Compile shaders
  igl::create_shader_program(scene_v_shader,scene_f_shader,{},scene_p_id);
  igl::create_shader_program(tex_v_shader,tex_f_shader,{},tex_p_id);
  // 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
  igl::read_triangle_mesh(argv[1],V,F);
  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);
      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);

      // Get read to draw quads
      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_ONE, 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:

depth peeling vs gl-blend

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_ONE, GL_ONE_MINUS_SRC_ALPHA);
glEnable(GL_BLEND);

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

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 <igl/read_triangle_mesh.h>
#include <Eigen/Core>
#include <GLUT/glut.h>
#include <string>
#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;
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");

  // 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 file
  igl::read_triangle_mesh(argv[1],V,F);
  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:

cheburashka modern opengl

and on my mac outputs:

OpenGL version: 4.1 INTEL-10.6.20

Generate list of k distant unique colors

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

rgb cube

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:

k unique colors in order

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

k unique colors in order

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.

k unique colors alternative;

Writing a mesh to an .obj file in a single line with Eigen

April 29th, 2015

I was recently revisiting our igl::writeOBJ code and realized that for simple meshes, the code for writing an .obj file can be really simple:

ofstream("test.obj")<<
  V.format(IOFormat(FullPrecision,0," ","\n","v ","","","\n"))<<
  (F.array()+1).format(IOFormat(FullPrecision,0," ","\n","f ","","","\n"));

Update: And because why not, here’s a .off writer:

ofstream("test.off")<<
  "OFF\n"<<V.rows()<<" "<<F.rows()<<" 0\n"<<
  V.format(IOFormat(FullPrecision,0," ","\n","","","","\n"))<<
  (F.array()).format(IOFormat(FullPrecision,0," ","\n","3 ","","","\n"));

Quick and dirty Eigen Matrix/Vector as std::map key

April 28th, 2015

I tried to use an Eigen matrix type as the key to a std::map with something like:

std::map<Eigen::VectorXd,bool> m;

But got compile-time errors of the sort:

error: no viable conversion from 'const CwiseBinaryOp<std::less<Scalar>, const Eigen::Array<double, -1, 1, 0, -1, 1>, const Eigen::Array<double, -1, 1, 0, -1, 1> >' to 'bool'
    {return __x < __y;}
            ^~~~~~~~~

Seems that map wants a proper less than operator and Eigen has overloaded that with the coefficient-wise operator. A reasonable way to sort vectors would be lexicographically. Fortunately stl has a function for that, so I define my map like this:

std::map<
  Eigen::VectorXd,
  bool,
  std::function<bool(const Eigen::VectorXd&,const Eigen::VectorXd&)> >
  m([](const VectorXd & a, const VectorXd & b)->bool
  {
    return std::lexicographical_compare(
      a.data(),a.data()+a.size(),
      b.data(),b.data()+b.size());
  });

This will ignore the internal ordering of the matrix elements (i.e. ColMajor vs RowMajor) but if you’re just using the map for a uniqueness check this is good enough.