Posts Tagged ‘eigen’

Trouble building opencv_contrib extras, can’t find unsupported/Eigen

Monday, May 23rd, 2016

I ran into a problem compiling the opencv_contrib modules.

Just running a vanilla cmake then make:

cmake -DOPENCV_EXTRA_MODULES_PATH=../opencv_contrib/modules ..
make

produced the error:

/Users/ajx/Downloads/opencv/opencv_contrib/modules/rgbd/src/odometry.cpp:41:10: fatal error: 'unsupported/Eigen/MatrixFunctions' file not found
#include <unsupported/Eigen/MatrixFunctions>

It seems the problem is that cmake was finding Eigen in /Library/Frameworks/Eigen.Framework/. I don’t even know how that Eigen got installed. Homebrew? Does it ship with Mac OS X now? In any case, I fixed the issue by pointing cmake directly to my homebrew’s Eigen installation:

cmake -DOPENCV_EXTRA_MODULES_PATH=../opencv_contrib/modules -DEIGEN_INCLUDE_PATH=/usr/local/include/eigen3/ ..

Using meshfix “from libigl”

Wednesday, March 30th, 2016

Here’s a little example of how to call meshfix using plain Eigen types to represent a mesh (i.e. in the libigl style. I’ve included this example in my github fork of the meshfix code.

#define MESHFIX_WITH_EIGEN
#include "meshfix.h"

#include <igl/read_triangle_mesh.h>
#include <igl/write_triangle_mesh.h>
#include <iostream>

int main(int argc, char * argv[])
{
  // Load in libigl's (V,F) format
  Eigen::MatrixXd V,W;
  Eigen::MatrixXi F,G;
  if(argc <= 2)
  {
    std::cout<<R"(Usage:
    ./meshfix-libigl [input](.obj|.ply|.stl|.off) [output](.obj|.ply|.stl|.off)
)";
    return EXIT_FAILURE;
  }
  igl::read_triangle_mesh(argv[1],V,F);
  meshfix(V,F,W,G);
  // Write to OBJ
  igl::write_triangle_mesh(argv[2],W,G);
}

Then you can run this by issuing something like:

./meshfix-libigl input.obj output.ply

Write a triangle mesh to xml file in libigl using Eigen matrix templated on CGAL’s Exact Kernel

Thursday, December 17th, 2015

Lately I’ve been working with CGAL and its arbitrary precision kernel. These are convenient, but writing to standard floating point precision mesh file formats requires rounding. Here’s a relatively straightforward way of using libigl‘s existing xml serialization routines to write an Eigen matrix tempalated on the CGAL Exact kernels number type (e.g. a matrix of vertex positions) to an xml file (using ASCII or base64 binary encoding).

Notice that there’s a bit of gnarly template specialization that needs to be done inside the igl::xml::serialization_xml namespace. The structure of the libigl xml serialization code is unfortunately not very flexible.

#include <igl/xml/serialize_xml.h>
#include <Eigen/Core>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <iostream>

typedef CGAL::Epeck::FT EScalar;
typedef Eigen::Matrix<EScalar,Eigen::Dynamic,Eigen::Dynamic> MatrixXE;

namespace igl
{
  namespace xml
  {
    namespace serialization_xml
    {
      template <> inline void serialize(
        const MatrixXE & obj,
        tinyxml2::XMLDocument* doc,
        tinyxml2::XMLElement* element,
        const std::string& name)
      {
        const std::function<std::string(const MatrixXE::Scalar &) > to_string = 
          [](const MatrixXE::Scalar & v)->std::string
          {
            return
              STR(CGAL::exact(v));
          };
        serialize(obj,name,to_string,doc,element);
      }
      template <> inline void deserialize(
        MatrixXE & obj,
        const tinyxml2::XMLDocument* doc,
        const tinyxml2::XMLElement* element,
        const std::string& name)
      {
        const std::function<void(const std::string &,MatrixXE::Scalar &)> & 
          from_string = 
          [](const std::string & s, MatrixXE::Scalar & v)
          {
            std::stringstream(s)>>v;
          };
        deserialize(doc,element,name,from_string,obj);
      }
    }
  }
}

int main(int argc, char * argv[])
{
  using namespace std;
  using namespace Eigen;
  // Little 4-vertex, 2-face mesh with rational coordinates
  MatrixXE V(4,3),W;
  V<<
    EScalar(1)/EScalar(3),  EScalar(1)/EScalar(13), EScalar(1)/EScalar(7),
    EScalar(1)/EScalar(7),  EScalar(1)/EScalar(3),  EScalar(1)/EScalar(13),
    EScalar(1)/EScalar(13), EScalar(1)/EScalar(7),  EScalar(1)/EScalar(3),
    EScalar(1)/EScalar(13), EScalar(1)/EScalar(7), -EScalar(1)/EScalar(3);
  MatrixXi F(2,3),G;
  F<<0,1,2,0,2,3;

  // Write mesh
  const bool binary = false;
  // Write vertices, overwriting file (true)
  igl::xml::serialize_xml(V,"vertices","exact.xml",binary,true);
  // Write faces to same file, appending (false)
  igl::xml::serialize_xml(F,"faces","exact.xml",binary,false);

  // Read mesh
  igl::xml::deserialize_xml(W,"vertices","exact.xml");
  igl::xml::deserialize_xml(G,"faces","exact.xml");

  // Verify 
  cout<<"V "<<(V.isApprox(W,0) ? "equals" : "does not equal")<<" W"<<endl;
}

Here’s a little feeling for the overhead of this format on a big single-precision mesh cast to the exact kernel (the expressions are simple and short, so this is sort of a best case scenario):

|           | .ply|.xml ascii|.xml binary| .obj|
|-----------|-----|----------|-----------|-----|
|File size: |172MB|     511MB|      330MB|288MB|
|Zip size:  |  68M|     134MB|       74MB| 81MB|
|Read time: |   8s|      270s|         9s|  32s|
|Write time:|  13s|       46s|        13s|  28s|

SSE illegal instruction run-time error in Eigen

Monday, November 30th, 2015

While debugging the other day we ran into a runtime error:


`
0x0000000000655679 in _mm_set_epi32 (__q0=0, __q1=0, __q2=0, __q3=0)
    at /usr/lib/gcc/x86_64-linux-gnu/4.9/include/emmintrin.h:595
595       return __extension__ (__m128i)(__v4si){ __q0, __q1, __q2, __q3 };
(gdb) bt
#0  0x0000000000655679 in _mm_set_epi32 (__q0=0, __q1=0, __q2=0, __q3=0)
    at /usr/lib/gcc/x86_64-linux-gnu/4.9/include/emmintrin.h:595
#1  _mm_set1_epi32 (__A=0) at /usr/lib/gcc/x86_64-linux-gnu/4.9/include/emmintrin.h:635
#2  Eigen::internal::pset1(Eigen::internal::unpacket_traits::type const&) (from=@0x7fffffffd61c: 0) at /usr/local/include/eigen3/Eigen/src/Core/arch/SSE/PacketMath.h:115
#3  0x00000000006a88ee in Eigen::internal::scalar_constant_op::packetOp (this=0x7fffffffd61c)
    at /usr/local/include/eigen3/Eigen/src/Core/Functors.h:522
#4  0x00000000006a0ea4 in Eigen::CwiseNullaryOp, Eigen::Matrix >::packet<0> (this=0x7fffffffd610, index=0)
    at /usr/local/include/eigen3/Eigen/src/Core/CwiseNullaryOp.h:88
#5  0x0000000000698582 in Eigen::DenseCoeffsBase, 1>::copyPacket, Eigen::Matrix >, 1, 0> (
    this=0x7fffffffd6a0, index=0, other=...) at /usr/local/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:537
#6  0x000000000068f403 in Eigen::internal::assign_impl, Eigen::CwiseNullaryOp, Eigen::Matrix >, 3, 0, 0>::run (dst=..., 
    src=...) at /usr/local/include/eigen3/Eigen/src/Core/Assign.h:410
#7  0x000000000068416e in Eigen::DenseBase >::lazyAssign, Eigen::Matrix > > (this=0x7fffffffd6a0, 
    other=...) at /usr/local/include/eigen3/Eigen/src/Core/Assign.h:499
#8  0x000000000067980d in Eigen::PlainObjectBase >::lazyAssign, Eigen::Matrix > > (
    this=0x7fffffffd6a0, other=...) at /usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h:414
#9  0x000000000066c1f8 in Eigen::internal::assign_selector, Eigen::CwiseNullaryOp, Eigen::Matrix >, false, false>::run (
    dst=..., other=...) at /usr/local/include/eigen3/Eigen/src/Core/Assign.h:520
#10 0x00000000006619b3 in Eigen::PlainObjectBase >::_set_noalias, Eigen::Matrix > > (
---Type  to continue, or q  to quit---
    this=0x7fffffffd6a0, other=...) at /usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h:621
#11 0x000000000069e5b7 in Eigen::PlainObjectBase >::_set_selector, Eigen::Matrix > > (
    this=0x7fffffffd6a0, other=...) at /usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h:606
#12 0x000000000069590c in Eigen::PlainObjectBase >::_set, Eigen::Matrix > > (this=0x7fffffffd6a0, 
    other=...) at /usr/local/include/eigen3/Eigen/src/Core/PlainObjectBase.h:598
#13 0x000000000068bc41 in Eigen::Matrix::operator=, Eigen::Matrix > > (this=0x7fffffffd6a0, other=...)
`

We were getting an “illegal instruction” trap. Finally we tracked down the issue to be that we’d compiled using the -msse4.2 flag on a machine that did not support SSE 4.2. In retrospect we should have suspected this from the header in which the error occurred (emmintrin.h). We had accidentally cross-compiled and then tried to run a bogus binary. To fix this issue we deleted the flag and replaced it with -march=native.

Memory gotcha in eigen

Wednesday, August 26th, 2015

I recently ran into a frustrating memory leak “gotcha” involving Eigen. The following code compiles

template <typename Derived>
Eigen::PlainObjectBase<Derived> leaky(
  const Eigen::PlainObjectBase<Derived> & A)
{
  Derived B = A;
  return B;
}

int main(int argc, char * argv[])
{
  Eigen::MatrixXd A(10000,1);
  Eigen::MatrixXd B = leaky(A);
}

but causes a memory leak (without any warnings/asserts from Eigen):

.main(33789,0x7fff7bb2d300) malloc: *** error for object 0x7f9ed188da00: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug

The problem seems to be the cast happening in the return. Changing the return type of leaky to Derived:

template <typename Derived>
Derived leaky(
  const Eigen::PlainObjectBase<Derived> & A)

fixes the issue. Or it could be fixed by changing the type of B inside leaky to Eigen::PlainObjectBase<Derived>.

  Eigen::PlainObjectBase<Derived> B = A;

Update:

The preferred answer from the Eigen developers is that leaky() should return Derived. This is a little disappointing if I want to have A and the return value be different types and I don’t want to expose all possible types for the return (i.e. I only want to template on top of Eigen matrices).

Sanity check that constructing a parameter inline circumvents aliasing functions with Eigen

Monday, July 27th, 2015

Suppose I have a function with Eigen types as inputs and outputs:

void aliasing_function(
  const Eigen::MatrixXd & A,
  Eigen::MatrixXd & B)
{
  for(int i = 0;i<B.rows();i++)
  {
    B.row(i) = A.row(A.rows()-i-1);
  }
}

This function is dangerous despite the const on the input parameter A because if A and B reference the same place in memory, writing to B will also write onto what this function thinks is A. Usually this is undesirable and it’s called aliasing. For example consider calling with the same input as output:

aliasing_function(C,C);

Supposing that I can’t/won’t/refuse to change the implementation of aliasing_function, here’s proof that you can avoid aliasing problems by using the copy constructor of the Eigen type directly in the parameter (avoiding the explicitly named temporary variable, but still costing a copy probably). The call above is replaced with:

aliasing_function(Eigen::MatrixXd(C),C);

To have a tangible example:

int main()
{
  using namespace std;
  Eigen::MatrixXd C = (Eigen::MatrixXd(3,3)<<1,2,3,4,5,6,7,8,9).finished();
  cout<<"C: "<<endl<<C<<endl;
  aliasing_function(Eigen::MatrixXd(C),C);
  cout<<"C: "<<endl<<C<<endl;
  aliasing_function(C,C);
  cout<<"C: "<<endl<<C<<endl;
}

This outputs:

C: 
1 2 3
4 5 6
7 8 9
C: 
7 8 9
4 5 6
1 2 3
C: 
1 2 3
4 5 6
1 2 3

The last print out shows the effect of aliasing.

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

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.

Counting non-zeros in boolean matrices

Monday, February 2nd, 2015

Sometimes I’ll use Eigen types with bool as the Scalar type. This invites “gotchas” when I assume the matrices will work like logicals do in MATLAB. Most immediately sum of logicals in matlab will return the integer number of trues (actually returns it as a double type). In Eigen, sum (not surprisingly) just calls operator+ and casts the result to bool again.

So rather than:

Matrix<bool,4,4> B;
B<<true,false,true,true;
B.sum(); // will just return true

if you want to compute the number of true entries use one of:

<del>B.nonZeros();</del>
B.count();
std::count(B.data(),B.data()+B.size(),true);

Things get trickier if your trying to use more advanced accessors. Consider trying to find the number of true entries per row. Instead of

B.rowwise().sum();

you could use

B.cast<int>().rowwise().sum();

but this casts the entire matrix first to int

B.rowwise().count();

or just roll up the loop over the rows.