Posts Tagged ‘opengl’

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.

Unzip OBJ-style mesh into a per-vertex attribute mesh

Wednesday, January 6th, 2016

The .obj mesh file format allows corners of triangles to pull attributes from different sources. The triangle:

v 0 0 0
v 1 0 0
v 0 0 0
vn 1 0 0
vn 0 1 0
vt 0 0 0
f 1/1/1 2/1/1 3/1/2

pulls vertex positions from entries (1,2,3) in the v ... vertex position list, texture coordinates from entries (1,1,1) in the vt ... list, and normals from entries (1,2,2) in the vn normals list.

If we think of corners being described by all attributes there are potentially #F*3 distinct corners. Often information is shared, and in the best case the position/texture/normal indices are all the same, so #V distinct corners. Usually it’s some mixture in between.

I added igl::unzip_corners to libigl which “unzips” an OBJ-style mesh into per-vertex attribute mesh: each new vertex is a distinct corner, so the new face list indexes all attributes in lock step (ideal for OpenGL). I was careful to determine uniqueness combinatorially so, for example, combinatorially distinct input vertices happening to share the same attributes (two vertices at the same position, etc.) don’t get merged (you could always use igl::remove_duplicate_vertices if you wanted to do that).

Here’s a little demo program:

  Eigen::MatrixXd V, TC, N;
  Eigen::MatrixXi F,FTC,FN;
  igl::readOBJ(argv[1],V,TC,N,F,FTC,FN);
  if(FTC.size() == 0)
  {
    FTC = F;
  }
  Eigen::MatrixXi U,G,J;
  igl::unzip_corners<Eigen::MatrixXi>({F,FTC},U,G,J);
  // New mesh vertices and texture coordinates indexed by G
  GV = igl::slice(Eigen::MatrixXd(V),U.col(0),1);
  GTC = igl::slice(Eigen::MatrixXd(TC),U.col(1),1);

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

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 <igl/read_triangle_mesh.h>
#include <igl/opengl/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::opengl::create_shader_program(scene_v_shader,scene_f_shader,{},scene_p_id);
  igl::opengl::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);
      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);

      // 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_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:

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_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 <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

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

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;

Old-style GPGPU reduction, average pixel color

Tuesday, April 28th, 2015

Here’s a little demo which computes the average pixel value of an OpenGL rendering. As a sanity check I compute the average value on the cpu-side using glReadPixels and then compare to computing the average value using a “mip-map” style, ping-pong, texture-buffer GPGPU reduction. Finally I render the buffers for fun.

I’m using yimg just to read in a .png file to render as an example.

On my weak little laptop, the GPU code is about 30 times faster. Can’t shrug at that! Rookie mistake. Made timings in debug mode. In release there’s hardly a speed up : – (

I was careful to compute the average in a coherent way (the images get progressively blurred out, rather than averaging the image quadrants recursively). This would be useful if, say, computing the average of 100 renderings. You could render them into a 10s x 10s texture and run the GPU reduction until the result is just a 10×10 image containing the 100 results. That’d only require a single final call to glReadPixels (rather than calling glReadPixels to read the final single pixel result of each reduction).

#include <YImage.hpp>
#include <YImage.cpp>
#include <GLUT/glut.h>
#include <iostream>
#include <string>
#include <iomanip>

// Size of image rounded up to next power of 2
size_t s;
// Shader program for doing reduction
GLuint reduction_prog_id;
// Need two textures and buffers to ping-pong
GLuint tex_id[] = {0,0};
GLuint fbo_id[] = {0,0};
GLuint dfbo_id[] = {0,0};
// Image (something to render in the first place)
std::string filename("hockney-512.png");
YImage yimg;

void init_render_to_texture(
  const size_t width,
  const size_t height,
  GLuint & tex_id,
  GLuint & fbo_id,
  GLuint & dfbo_id)
{
  using namespace std;
  // Delete if already exists
  glDeleteTextures(1,&tex_id);
  glDeleteFramebuffersEXT(1,&fbo_id);
  glDeleteFramebuffersEXT(1,&dfbo_id);
  // http://www.opengl.org/wiki/Framebuffer_Object_Examples#Quick_example.2C_render_to_texture_.282D.29
  glGenTextures(1, &tex_id);
  glBindTexture(GL_TEXTURE_2D, tex_id);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
  //NULL means reserve texture memory, but texels are undefined
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F_ARB, width, height, 0, GL_BGRA, GL_FLOAT, NULL);
  glBindTexture(GL_TEXTURE_2D, 0);
  glGenFramebuffersEXT(1, &fbo_id);
  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, fbo_id);
  //Attach 2D texture to this FBO
  glFramebufferTexture2DEXT(GL_FRAMEBUFFER_EXT, GL_COLOR_ATTACHMENT0_EXT, GL_TEXTURE_2D, tex_id, 0);
  glGenRenderbuffersEXT(1, &dfbo_id);
  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, dfbo_id);
  glRenderbufferStorageEXT(GL_RENDERBUFFER_EXT, GL_DEPTH_COMPONENT24, width, height);
  //Attach depth buffer to FBO (for this example it's not really needed, but if
  //drawing a 3D scene it would be necessary to attach something)
  glFramebufferRenderbufferEXT(GL_FRAMEBUFFER_EXT, GL_DEPTH_ATTACHMENT_EXT, GL_RENDERBUFFER_EXT, dfbo_id);
  //Does the GPU support current FBO configuration?
  GLenum status;
  status = glCheckFramebufferStatusEXT(GL_FRAMEBUFFER_EXT);
  assert(status == GL_FRAMEBUFFER_COMPLETE_EXT);
  // Unbind to clean up
  glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, 0);
  glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, 0);
}

int main(int argc, char * argv[])
{
  glutInit(&argc,argv);
  if(argc>1)
  {
    filename = argv[1];
  }
  yimg.load(filename.c_str());
  s = std::max(yimg.width(),yimg.height());
  // http://stackoverflow.com/a/466278/148668
  s--;
  s |= s >> 1;
  s |= s >> 2;
  s |= s >> 4;
  s |= s >> 8;
  s |= s >> 16;
  s++;

  glutInitDisplayString("rgba depth double stencil");
  glutInitWindowSize(2*s,s);
  glutCreateWindow("texture-reduction");
  glutDisplayFunc(
    []()
    {
      using namespace std;
      // Initialize **both** buffers and set to render into first
      for(int i = 1;i>=0;i--)
      {
        glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, fbo_id[i]);
        glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, dfbo_id[i]);
        glClearColor(0,0,0,1);
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
      }

      // Render something
      yimg.flip();
      glMatrixMode(GL_PROJECTION);
      glPushMatrix();
      glLoadIdentity();
      glOrtho(0,s,0,s, -10000,10000);
      glViewport(0,0,s,s);
      glRasterPos2f(0,0);
      glDrawPixels(yimg.width(), yimg.height(), GL_RGBA, GL_UNSIGNED_BYTE, yimg.data());
      glPopMatrix();

      // Even the cpu code should use a buffer (rather than the screen)
      GLfloat * rgb = new GLfloat[yimg.width() * yimg.height() * 3];
      glReadPixels(0, 0, yimg.width(), yimg.height(), GL_RGB, GL_FLOAT, rgb);
      // Gather into double: sequential add is prone to numerical error
      double avg[] = {0,0,0};
      for(size_t i = 0;i<yimg.width()*yimg.height()*3;i+=3)
      {
        avg[0] += rgb[i + 0];
        avg[1] += rgb[i + 1];
        avg[2] += rgb[i + 2];
      }
      for_each(avg,avg+3,[](double & c){c/=yimg.width()*yimg.height();});
      delete[] rgb;

      // Size of square being rendered
      assert(((s != 0) && ((s & (~s + 1)) == s)) && "s should be power of 2");
      size_t h = s/2;
      // odd or even ping-pong iteration
      int odd = 1;
      // Tell shader about texel step size
      glUseProgram(reduction_prog_id);
      glUniform1f(glGetUniformLocation(reduction_prog_id,"dt"),0.5f/float(s));
      glEnable(GL_TEXTURE_2D);
      while(h)
      {
        // Select which texture to draw/compute into
        glBindFramebufferEXT(GL_FRAMEBUFFER_EXT, fbo_id[odd]);
        glBindRenderbufferEXT(GL_RENDERBUFFER_EXT, dfbo_id[odd]);
        // Select which texture to draw/compute from
        glBindTexture(GL_TEXTURE_2D,tex_id[1-odd]);
        // Scale to smaller square
        glViewport(0,0,h,h);
        const float f = 2.*(float)h/(float)s;
        // Draw quad filling viewport with shrinking texture coordinates
        glBegin(GL_QUADS);
        glTexCoord2f(0,0);
        glVertex2f  (-1,-1);
        glTexCoord2f(f,0);
        glVertex2f  (1,-1);
        glTexCoord2f(f,f);
        glVertex2f  (1,1);
        glTexCoord2f(0,f);
        glVertex2f  (-1,1);
        glEnd();

        // ping-pong
        odd = 1-odd;
        h = h/2;
      }
      // Read corner pixel of last render
      float px[3];
      glReadPixels(0, 0, 1, 1, GL_RGB, GL_FLOAT, px);
      // Correct for size not power of 2
      for_each(px,px+3,[](float& c){c*=(float)s*s/yimg.width()/yimg.height();});
      cout<<" gpu: "<< px[0]<<" "<< px[1]<<" "<< px[2]<<" "<<endl;
      cout<<" cpu: "<<avg[0]<<" "<<avg[1]<<" "<<avg[2]<<" "<<endl;

      // Purely for vanity, draw the buffers
      glUseProgram(0);
      glBindFramebufferEXT(GL_FRAMEBUFFER_EXT,0);
      glBindRenderbufferEXT(GL_RENDERBUFFER_EXT,0);
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
      for(odd = 0;odd<2;odd++)
      {
        glViewport(odd*s,0,s,s);
        glEnable(GL_TEXTURE_2D);
        glBindTexture(GL_TEXTURE_2D,tex_id[odd]);
        glBegin(GL_QUADS);
        glTexCoord2f(0,0);
        glVertex2f  (-1,-1);
        glTexCoord2f(1,0);
        glVertex2f  (1,-1);
        glTexCoord2f(1,1);
        glVertex2f  (1,1);
        glTexCoord2f(0,1);
        glVertex2f  (-1,1);
        glEnd();
      }
      glutSwapBuffers();
    }
   );
  init_render_to_texture(s,s,tex_id[0],fbo_id[0],dfbo_id[0]);
  init_render_to_texture(s,s,tex_id[1],fbo_id[1],dfbo_id[1]);


  // Vertex shader is a **true** pass-through
  const std::string vertex_shader = R"(
#version 120
void main()
{
  gl_Position = gl_Vertex;
  gl_TexCoord[0] = gl_MultiTexCoord0;
}
)";
  // fragment shader sums texture of this pixel and left/bottom neighbors
  const std::string fragment_shader = R"(
#version 120
// size of a half-texel in 1/pixels units: 0.5/(full size)
uniform float dt;
uniform sampler2D texture;
void main()
{
  gl_FragColor = 0.25*(
    texture2D(texture, gl_TexCoord[0].st - vec2( 0, 0)) +
    texture2D(texture, gl_TexCoord[0].st - vec2(dt, 0)) +
    texture2D(texture, gl_TexCoord[0].st - vec2(dt,dt)) +
    texture2D(texture, gl_TexCoord[0].st - vec2( 0,dt)));
}
)";

  // Compile and link reduction shaders into program
  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());
  reduction_prog_id = glCreateProgram();
  glAttachShader(reduction_prog_id,vid);
  glAttachShader(reduction_prog_id,fid);
  glLinkProgram(reduction_prog_id);

  glutMainLoop();
}

gluUnProject buggy for canonical views

Sunday, March 15th, 2015

I just ran into a very frustrating bug in gluUnProject. It seems that if the model view matrix is rotated to align the view axis with the y-axis then gluUnProject is unreliable. Here’s a little example:

Set up you model view, projection and viewport with:

MV = [1 0 0 0;0 0 -1 0;0 1 0 -4.3;0 0 0 1];
P = [1.81066 0 0 0;0 2.41421 0 0; 0 0 -1.00019 -0.0200019; 0 0 -1 0];
VP = [0 0 960 720];

Then given a point in the scene

p = [1.66667 0 0];

project it using gluProject (this seems to work OK):

proj_p = [816.867 360 0.99777];

move it one pixel up (as if with the mouse):

proj_p = [816.867 361 0.99777];

then unproject with gluUnProject:

wrong = [1.80114 -0.345327 -0.00534674];

Where did that movement in the y-direction come from? My guess is the implementation of gluUnProject is doing a bad job inverting the scene matrix in terms of floating point error. Luckily, libigl has an implementation of igl::unproject which uses Eigen’s robust inverse() function. This solved my problem and gives the correct unprojection:

correct = [1.66667 0 -0.00494755];

Understanding OpenGL screen z (depth) values

Saturday, February 8th, 2014

I was trying to unproject a depth image into an OpenGL and was bit surprised to find that my choice of near and far plane values affected the influence of my “window Z” values used to unproject. This really shouldn’t have been a surprise when considering how the perspective matrix is formed. If you’re using gluPersepctive then your projection matrix P is:

P =  / x 0 0           0         \
    |  0 y 0           0          |
    |  0 0 (f+n)/(n-f) 2fn/(n-f)  |
     \ 0 0 -1          0         /

where x and y control the field of view angle and aspect ratio, and f and n are the far and near clipping depths respectively.

The projection of a world point z = (0 0 z 1)T is then given by computing p = P * z then applying perspective division q = p / p.w. Despite the incorect documentation, gluProject then computes a window “depth” via winZ = (q.z+1)/2.

Following these formulae we can derive a relationship between winZ and the z-coordinate of our point on the z-axis. Namely,

winZ = (p.z/p.w+1)/2
p.z/p.w = 2*winZ-1
p.w = -z
p.z = z*(f+n)/(n-f) + 2fn/(n-f)
(z(f+n)/(n-f) + 2fn/(n-f))/(-z) = 2winZ-1
z = fn/(f*winZ-n*winZ-f)

or equivalently

winZ = f(n+z)/((f-n)z)

Notably, when z=-f then winZ = 1 and when z=-n then winZ = 0. In between, we have an inverse relationship:

winz gluproject

Speed up slow scatter and line plots in MATLAB

Tuesday, February 4th, 2014

Matlab seems to default to its CPU renderer 'painters' when creating new figures with scattered dots or lines. So

plot3(X,Y,Z);

where X and so on are very large can take a long time (say 1 minute).

Whereas first setting the renderer to the GPU renderer 'OpenGL' is much faster.

set(gcf,'Renderer','OpenGL');
plot3(X,Y,Z);

for the same large X takes 6 seconds.

Awkwardly, plotting a bunch of (disconnected) line segments from vertices in N and edge indices in E in matlab using something like:

plot3([N(E(:,1),1) N(E(:,2),1)]',[N(E(:,1),2) N(E(:,2),2)]',[N(E(:,1),3) N(E(:,2),3)]','k');

or

line([N(E(:,1),1) N(E(:,2),1)]',[N(E(:,1),2) N(E(:,2),2)]',[N(E(:,1),3) N(E(:,2),3)]');

is much slower to render and interact with than to just draw a bunch of degenerate triangles using:

trisurf([E E(:,1)],N(:,1),N(:,2),N(:,3))

This has the added bonus that it seems to throw the renderer automatically in to ‘OpenGL’ mode.