Archive for April, 2015

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

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

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

CGAL fails to compile intersection code only in debug mode.

Tuesday, April 28th, 2015

I tried to compile:

#include <CGAL/assertions.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/intersections.h>

typedef CGAL::Epeck Kernel;
typedef Kernel::FT CGALScalar;
typedef CGAL::Triangle_3<Kernel> Triangle_3; 

int main(int argc, char * argv[])
{
  Triangle_3 A,B;
  CGAL::do_intersect(A,B);
  CGAL::intersection(A,B);
}

But got stream of CGAL template error diarrhea:


In file included from ./bar.cpp.cpp:44:
In file included from /opt/local/include/CGAL/Exact_predicates_exact_constructions_kernel.h:28:
In file included from /opt/local/include/CGAL/Simple_cartesian.h:28:
In file included from /opt/local/include/CGAL/Cartesian/Cartesian_base.h:68:
In file included from /opt/local/include/CGAL/Cartesian/function_objects.h:28:
In file included from /opt/local/include/CGAL/Kernel/function_objects.h:34:
In file included from /opt/local/include/CGAL/intersection_3.h:34:
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:72:7: error: no matching function for call to 'possibly'
      CGAL_kernel_assertion(obj);
      ^~~~~~~~~~~~~~~~~~~~~~~~~~
/opt/local/include/CGAL/kernel_assertions.h:50:5: note: expanded from macro 'CGAL_kernel_assertion'
   (CGAL::possibly(EX)?(static_cast<void>(0)): ::CGAL::assertion_fail( # EX , __FILE__, __LINE__))
    ^~~~~~~~~~~~~~
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:108:3: note: in instantiation of function template
      specialization
      'CGAL::internal::intersection_coplanar_triangles_cutoff<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >'
      requested here
  intersection_coplanar_triangles_cutoff(p,q,r,k,inter_pts); //line pq
  ^
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:190:12: note: in instantiation of function template
      specialization 'CGAL::internal::intersection_coplanar_triangles<CGAL::Simple_cartesian<CGAL::Interval_nt<false> >
      >' requested here
    return intersection_coplanar_triangles(t1,t2,k);
           ^
/opt/local/include/CGAL/Kernel/function_objects.h:2577:24: note: in instantiation of function template specialization
      'CGAL::internal::intersection<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >' requested here
    { return internal::intersection(t1, t2, K() ); }
                       ^
/opt/local/include/CGAL/Lazy.h:440:31: note: in instantiation of function template specialization
      'CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> >
      >::operator()<CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > > >' requested here
BOOST_PP_REPEAT_FROM_TO(2, 9, CGAL_LAZY_REP, _)
                              ^
/opt/local/include/boost/preprocessor/repetition/repeat_from_to.hpp:36:125: note: expanded from macro
      'BOOST_PP_REPEAT_FROM_TO_1'
  ...l, m, dt) BOOST_PP_REPEAT_FROM_TO_D_1(BOOST_PP_AUTO_REC(BOOST_PP_WHILE_P, 256), f, l, m, dt)
                                                                                           ^
/opt/local/include/boost/preprocessor/repetition/repeat_from_to.hpp:54:136: note: expanded from macro
      'BOOST_PP_REPEAT_FROM_TO_D_1'
  ...f, l, m, dt) BOOST_PP_REPEAT_1(BOOST_PP_SUB_D(d, l, f), BOOST_PP_REPEAT_FROM_TO_M_1, (d, f, m, dt))
                                                                                                 ^
/opt/local/include/boost/preprocessor/repetition/repeat.hpp:38:63: note: expanded from macro 'BOOST_PP_REPEAT_1'
# define BOOST_PP_REPEAT_1(c, m, d) BOOST_PP_REPEAT_1_I(c, m, d)
                                                              ^
note: (skipping 11 expansions in backtrace; use -fmacro-backtrace-limit=0 to see all)
/opt/local/include/boost/preprocessor/repetition/repeat_from_to.hpp:79:118: note: expanded from macro
      'BOOST_PP_REPEAT_FROM_TO_M_1_I'
  ...n, d, f, m, dt) BOOST_PP_REPEAT_FROM_TO_M_1_II(z, BOOST_PP_ADD_D(d, n, f), m, dt)
                                                                                ^
/opt/local/include/boost/preprocessor/repetition/repeat_from_to.hpp:83:54: note: expanded from macro
      'BOOST_PP_REPEAT_FROM_TO_M_1_II'
# define BOOST_PP_REPEAT_FROM_TO_M_1_II(z, n, m, dt) m(z, n, dt)
                                                     ^
/opt/local/include/CGAL/Lazy.h:434:29: note: expanded from macro 'CGAL_LAZY_REP'
    : Lazy_rep<AT, ET, E2A>(ac( BOOST_PP_ENUM(n, CGAL_LN, CGAL::approx) )), BOOST_PP_ENUM(n, CGAL_LINIT, _) \
                            ^
/opt/local/include/CGAL/Lazy.h:1552:39: note: in instantiation of member function
      'CGAL::Lazy_rep_2<boost::optional<boost::variant<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> >
      >, CGAL::Segment_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      std::__1::vector<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      std::__1::allocator<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > > > >,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_> >, boost::optional<boost::variant<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Gmpq>
      >, CGAL::Segment_3<CGAL::Simple_cartesian<CGAL::Gmpq> >, CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> >,
      std::__1::vector<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Gmpq> >,
      std::__1::allocator<CGAL::Point_3<CGAL::Simple_cartesian<CGAL::Gmpq> > > >, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_,
      boost::detail::variant::void_, boost::detail::variant::void_, boost::detail::variant::void_> >,
      CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq> >,
      CGAL::Cartesian_converter<CGAL::Simple_cartesian<CGAL::Gmpq>, CGAL::Simple_cartesian<CGAL::Interval_nt<false> >,
      CGAL::NT_converter<CGAL::Gmpq, CGAL::Interval_nt<false> > >, CGAL::Triangle_3<CGAL::Epeck>,
      CGAL::Triangle_3<CGAL::Epeck> >::Lazy_rep_2' requested here
      Lazy<AT, ET, EFT, E2A> lazy(new Lazy_rep_2<AT, ET, AC, EC, E2A, L1, L2>(AC(), EC(), l1, l2));
                                      ^
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:230:1: note: in instantiation of function template
      specialization 'CGAL::Lazy_construction_variant<CGAL::Epeck,
      CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq> >
      >::operator()<CGAL::Triangle_3<CGAL::Epeck>, CGAL::Triangle_3<CGAL::Epeck> >' requested here
CGAL_INTERSECTION_FUNCTION_SELF(Triangle_3, 3)
^
/opt/local/include/CGAL/Intersection_traits.h:88:25: note: expanded from macro 'CGAL_INTERSECTION_FUNCTION_SELF'
    return BOOST_PP_CAT(K().intersect_, BOOST_PP_CAT(DIM, _object()(a, b))); \
                        ^
/opt/local/include/boost/preprocessor/cat.hpp:22:47: note: expanded from macro 'BOOST_PP_CAT'
#    define BOOST_PP_CAT(a, b) BOOST_PP_CAT_I(a, b)
                                              ^
/opt/local/include/boost/preprocessor/cat.hpp:29:34: note: expanded from macro 'BOOST_PP_CAT_I'
#    define BOOST_PP_CAT_I(a, b) a ## b
                                 ^
./bar.cpp.cpp:56:9: note: in instantiation of function template specialization 'CGAL::intersection<CGAL::Epeck>'
      requested here
  CGAL::intersection(A,B);
        ^
/opt/local/include/CGAL/Uncertain.h:263:13: note: candidate function not viable: no known conversion from 'typename
      Intersection_traits<Simple_cartesian<Interval_nt<false> >, typename Simple_cartesian<Interval_nt<false>
      >::Line_3, typename Simple_cartesian<Interval_nt<false> >::Line_3>::result_type' (aka 'optional<variant_type>')
      to 'bool' for 1st argument
inline bool possibly(bool b) { return b; }
            ^
/opt/local/include/CGAL/Uncertain.h:275:6: note: candidate function not viable: no known conversion from 'typename
      Intersection_traits<Simple_cartesian<Interval_nt<false> >, typename Simple_cartesian<Interval_nt<false>
      >::Line_3, typename Simple_cartesian<Interval_nt<false> >::Line_3>::result_type' (aka 'optional<variant_type>')
      to 'Uncertain<bool>' for 1st argument
bool possibly(Uncertain<bool> c)
     ^
In file included from ./bar.cpp.cpp:44:
In file included from /opt/local/include/CGAL/Exact_predicates_exact_constructions_kernel.h:28:
In file included from /opt/local/include/CGAL/Simple_cartesian.h:28:
In file included from /opt/local/include/CGAL/Cartesian/Cartesian_base.h:68:
In file included from /opt/local/include/CGAL/Cartesian/function_objects.h:28:
In file included from /opt/local/include/CGAL/Kernel/function_objects.h:34:
In file included from /opt/local/include/CGAL/intersection_3.h:34:
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:72:7: error: no matching function for call to 'possibly'
      CGAL_kernel_assertion(obj);
      ^~~~~~~~~~~~~~~~~~~~~~~~~~
/opt/local/include/CGAL/kernel_assertions.h:50:5: note: expanded from macro 'CGAL_kernel_assertion'
   (CGAL::possibly(EX)?(static_cast<void>(0)): ::CGAL::assertion_fail( # EX , __FILE__, __LINE__))
    ^~~~~~~~~~~~~~
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:108:3: note: in instantiation of function template
      specialization 'CGAL::internal::intersection_coplanar_triangles_cutoff<CGAL::Simple_cartesian<CGAL::Gmpq> >'
      requested here
  intersection_coplanar_triangles_cutoff(p,q,r,k,inter_pts); //line pq
  ^
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:190:12: note: in instantiation of function template
      specialization 'CGAL::internal::intersection_coplanar_triangles<CGAL::Simple_cartesian<CGAL::Gmpq> >' requested
      here
    return intersection_coplanar_triangles(t1,t2,k);
           ^
/opt/local/include/CGAL/Kernel/function_objects.h:2577:24: note: in instantiation of function template specialization
      'CGAL::internal::intersection<CGAL::Simple_cartesian<CGAL::Gmpq> >' requested here
    { return internal::intersection(t1, t2, K() ); }
                       ^
/opt/local/include/CGAL/Lazy.h:1573:20: note: in instantiation of function template specialization
      'CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq>
      >::operator()<CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> >,
      CGAL::Triangle_3<CGAL::Simple_cartesian<CGAL::Gmpq> > >' requested here
      ET exact_v = EC()(CGAL::exact(l1), CGAL::exact(l2));
                   ^
/opt/local/include/CGAL/Triangle_3_Triangle_3_intersection.h:230:1: note: in instantiation of function template
      specialization 'CGAL::Lazy_construction_variant<CGAL::Epeck,
      CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Interval_nt<false> > >,
      CGAL::CommonKernelFunctors::Intersect_3<CGAL::Simple_cartesian<CGAL::Gmpq> >
      >::operator()<CGAL::Triangle_3<CGAL::Epeck>, CGAL::Triangle_3<CGAL::Epeck> >' requested here
CGAL_INTERSECTION_FUNCTION_SELF(Triangle_3, 3)
^
/opt/local/include/CGAL/Intersection_traits.h:88:25: note: expanded from macro 'CGAL_INTERSECTION_FUNCTION_SELF'
    return BOOST_PP_CAT(K().intersect_, BOOST_PP_CAT(DIM, _object()(a, b))); \
                        ^
/opt/local/include/boost/preprocessor/cat.hpp:22:47: note: expanded from macro 'BOOST_PP_CAT'
#    define BOOST_PP_CAT(a, b) BOOST_PP_CAT_I(a, b)
                                              ^
/opt/local/include/boost/preprocessor/cat.hpp:29:34: note: expanded from macro 'BOOST_PP_CAT_I'
#    define BOOST_PP_CAT_I(a, b) a ## b
                                 ^
./bar.cpp.cpp:56:9: note: in instantiation of function template specialization 'CGAL::intersection<CGAL::Epeck>'
      requested here
  CGAL::intersection(A,B);
        ^
/opt/local/include/CGAL/Uncertain.h:263:13: note: candidate function not viable: no known conversion from 'typename
      Intersection_traits<Simple_cartesian<Gmpq>, typename Simple_cartesian<Gmpq>::Line_3, typename
      Simple_cartesian<Gmpq>::Line_3>::result_type' (aka 'optional<variant_type>') to 'bool' for 1st argument
inline bool possibly(bool b) { return b; }
            ^
/opt/local/include/CGAL/Uncertain.h:275:6: note: candidate function not viable: no known conversion from 'typename
      Intersection_traits<Simple_cartesian<Gmpq>, typename Simple_cartesian<Gmpq>::Line_3, typename
      Simple_cartesian<Gmpq>::Line_3>::result_type' (aka 'optional<variant_type>') to 'Uncertain<bool>' for 1st
      argument
bool possibly(Uncertain<bool> c)

Compiling with either -DNDEBUG -DCGAL_NDEBUG hides these errors. After a lot of digging, I found a small note in the CGAL doc:

The CGAL::intersection() function used to return an Object, but starting with CGAL 4.2 the return type is determined by a metafunction defined by the kernel. To preserve backward compatibility Object can be constructed from the new return types implicitly, but switching to the new style is recommended. To enable the old style without any overhead, the macro CGAL_INTERSECTION_VERSION must be defined to 1 before any CGAL header is included.

Indeed adding

#define CGAL_INTERSECTION_VERSION 1

Fixes these errors. The doc is pretty confusing and this example is riddled with errors. Just adding auto result = in front of CGAL::intersection(A,B); doesn’t make the errors above go away.

Update: This will also avoid the problem above:

#define CGAL_KERNEL_NO_ASSERTIONS

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

Edge collapse and mesh decimation in libigl

Thursday, April 16th, 2015

Libigl purposefully does not build itself around complicated mesh data-structures. This is freeing for many reasons: it’s very easy to include our code in an arbitrary project, functions are not artificially limited to manifold meshes, array based data structures are easily serialized and fast, matrix operations on vertex positions are directly exposed, etc.

mesh decimation in libigl

But our choice is also limiting. In particular, combinatorial changes to the mesh are potentially expensive and difficult to implement. Recently, I took a first stab at implementing an efficient edge-collapse routine for libigl. Indeed I’m seeing good performance: O(1) constant time for a single edge collapse and O(log m) time for cost-priority-queue-based collapses. The data structures are a little “messy” in the sense that when edges or faces disappear their rows are not deleted, but rather just redacted: entries are replaced with NULL values. This is because my approach is array-based and constant resizing would ruin performance. Luckily for edge-collapse, the number of elements only decreases. For edge-split I’ll have to think about amortized costs…

For now, check out the updated tutorial entry for libigl’s new mesh decimation and edge collapse features.

The code is “programmable” in the sense I expose function handles for computing edge costs and merged vertex placements. Though the default functions I currently provide are quite naive, this should support rather advanced “Progressive Meshes” style simplification with creases, sharp features, adaptivity, etc.

Port and Modification of Pascal Frey’s medit software for tet-mesh visualization on github

Wednesday, April 15th, 2015

I’m moving my patched and improved version of Pascal Frey’s medit software from libigl/external to its own github repository. The new version still depends on libigl and AntTweakbar. There’s a matlab wrapper for it in my gptoolbox.

This tool has been essential to my research on tetrahedral meshes. There are more advanced tet-mesh, 3d-slicing visualization tools around, but this one is straightforward and open-source. I’ve made some modifications to make it easier to visualization data and make more complicated slices (hot-dog view!). The input is the same .mesh tet mesh format used by TetGen.

medit screen capture

qslim for Mac OS X on github

Wednesday, April 15th, 2015

A while ago I patched Michael Garland’s qslim to compile on modern Mac OS X. I’ve finally wrangled in the instructions and modified files to a github repo. There are also wrappers for this tool in my gptoolbox. It’d be cool to have a proper C++ API to this tool, but that’s future work.

Qslim is one of the seminal mesh decimation/edge-collapse/simplification techniques. It’s fast, simple and predictable. It’s a good starting place if you’re trying to simplify your mesh, especially in an automated way.

DEC using gptoolbox

Tuesday, April 14th, 2015

Speaking with the Caltech Discrete Exterior Calculus (DEC) promulgators, it seems that despite using DEC notation to explain their math they advocate for building the final algebraic entities (discrete cotangent Laplacian) directly. That is, rather than by multiplying the core DEC stars and d’s.

Nonetheless, when re-implementing these papers I find it useful to follow the texts step-by-step. To do that, I build DEC operators (differentials) using functions in my gptoolbox for MATLAB. I’ll use D0 for d1 and S0 for ★0 and so on:

% List of all "half"-edges: 3*#F by 2
allE = [F(:,[2 3]); F(:,[3 1]); F(:,[1 2])];
% Sort each edge and remember if order changed
[sortallE,I] = sort(allE,2);
% 1 if edge matches half-edge orientation, -1 otherwise
flip = 3-2*I(:,1);
% find unique edges. EMAP(i) tells us where to
% find sortallE(i,:) in uE: so that sortallE(i,:) = uE(EMAP(i),:)
[E,~,EMAP] = unique(sortallE,'rows');
% number of edges
ne = size(E,1);
% number of vertices
nv = size(V,1);
% number of faces
nf = size(F,1);
% S0  nv by nv lumped diagonal mass matrix
S0 = massmatrix(V,F);
% S1  ne by ne diagonal matrix of edge dual-primal ratios (cotangents)
C = cotangent(V,F);
S1 = sparse(EMAP,EMAP,-C(:),ne,ne);
% S2  nf by nf diagonal matrix of inverted face areas
S2 = diag(sparse(doublearea(V,F).^-1))/2;
% D0  ne by nv edge to vertex signed incidence matrix
D0 = sparse(repmat(1:ne,2,1)',E,repmat([1 -1],ne,1),ne,nv);
% D1  nf by ne edge to face signed incidence matrix
D1 = sparse(repmat(1:nf,1,3)',EMAP,flip,nf,ne);

Then you can build familiar matrices like the discrete Laplacian:

L = D0' * S1 * D0;

The analogous code for building these in C++ libigl should hopefully be obvious.

Battle with acmsiggraph.cls’s copyright space

Monday, April 13th, 2015

Switching to the final version of the acm siggraph latex template always messes up the first page, but lately I’ve had a lot of trouble with the copyright space. Somehow I end up with blank space all other the first two columns rather than just a small portion at the bottom of the first column. To fix this, in acmsiggraph.cls I replaced:

\footnotetext[0]{\rule[\ACMcopyrightspace]{2.71828in}{0in}}%

with

\footnotetext[0]{\vspace*{1in}}%

MeshFix on github

Monday, April 13th, 2015

I’ve moved the meshfix port I made for mac from its buried location inside libigl to its own github repo.

If you’re commonly dealing with bad meshes (self-intersections, non-manifold bits, holes, etc.). This is a great fully automatic tool. I often use it before sending meshes to tetgen. It works particularly for single component meshes with lots of little problems.