Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: 58754475350b80f3567c9ae09f4913e1632ffaf3

include/igl/Camera.h

@@ -165,7 +165,6 @@ inline Eigen::Matrix4d igl::Camera::projection() const
 {
   Eigen::Matrix4d P;
   using namespace std;
-  using namespace igl;
   // http://stackoverflow.com/a/3738696/148668
   if(m_angle >= IGL_CAMERA_MIN_ANGLE)
   {
@@ -244,7 +243,6 @@ inline Eigen::Vector3d igl::Camera::up() const
 
 inline Eigen::Vector3d igl::Camera::unit_plane() const
 {
-  using namespace igl;
   // Distance of center pixel to eye
   const double d = 1.0;
   const double a = m_aspect;
@@ -263,7 +261,6 @@ inline void igl::Camera::dolly(const Eigen::Vector3d & dv)
 inline void igl::Camera::push_away(const double s)
 {
   using namespace Eigen;
-  using namespace igl;
 #ifndef NDEBUG
   Vector3d old_at = at();
 #endif
@@ -276,7 +273,6 @@ inline void igl::Camera::push_away(const double s)
 inline void igl::Camera::dolly_zoom(const double da)
 {
   using namespace std;
-  using namespace igl;
   using namespace Eigen;
 #ifndef NDEBUG
   Vector3d old_at = at();
@@ -316,7 +312,6 @@ inline void igl::Camera::dolly_zoom(const double da)
 inline void igl::Camera::turn_eye(const Eigen::Quaterniond & q)
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d old_eye = eye();
   // eye should be fixed
   //
@@ -330,7 +325,6 @@ inline void igl::Camera::turn_eye(const Eigen::Quaterniond & q)
 inline void igl::Camera::orbit(const Eigen::Quaterniond & q)
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d old_at = at();
   // at should be fixed
   //
@@ -351,7 +345,6 @@ inline void igl::Camera::look_at(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   // http://www.opengl.org/sdk/docs/man2/xhtml/gluLookAt.xml
   // Normalize vector from at to eye
   Vector3d F = eye-at;

include/igl/InElementAABB.h

@@ -111,7 +111,6 @@ inline void igl::InElementAABB::init(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   if(bb_mins.size() > 0)
   {
     assert(bb_mins.rows() == bb_maxs.rows() && "Serial tree arrays must match");
@@ -171,7 +170,6 @@ inline void igl::InElementAABB::init(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   const int dim = V.cols();
   const double inf = numeric_limits<double>::infinity();
   m_bb_min.setConstant(1,dim, inf);
@@ -274,7 +272,6 @@ inline std::vector<int> igl::InElementAABB::find(
     const bool first) const
 {
   using namespace std;
-  using namespace igl;
   using namespace Eigen;
   bool inside = true;
   const int dim = m_bb_max.size();

include/igl/MouseController.h

@@ -214,7 +214,6 @@ inline void igl::MouseController::reshape(const int w, const int h)
 inline bool igl::MouseController::down(const int x, const int y)
 {
   using namespace std;
-  using namespace igl;
   m_down_x = m_drag_x =x;
   m_down_y = m_drag_y =y;
   const bool widget_down = any_selection() && m_widget.down(x,m_height-y);

include/igl/RotateWidget.h

@@ -118,7 +118,6 @@ inline Eigen::Quaterniond igl::RotateWidget::axis_q(const int a)
 inline Eigen::Vector3d igl::RotateWidget::view_direction(const int x, const int y)
 {
   using namespace Eigen;
-  using namespace igl;
   const Vector3d win_s(x,y,0), win_d(x,y,1);
   const Vector3d s = unproject(win_s);
   const Vector3d d = unproject(win_d);
@@ -128,7 +127,6 @@ inline Eigen::Vector3d igl::RotateWidget::view_direction(const int x, const int
 inline Eigen::Vector3d igl::RotateWidget::view_direction(const Eigen::Vector3d & pos)
 {
   using namespace Eigen;
-  using namespace igl;
   const Vector3d ppos = project(pos);
   return view_direction(ppos(0),ppos(1));
 }
@@ -151,7 +149,6 @@ inline Eigen::Vector3d igl::RotateWidget::unproject_onto(
   const int y) const
 {
   using namespace Eigen;
-  using namespace igl;
   // KNOWN BUG: This projects to same depths as pos. I think what we actually
   // want is The intersection with the plane perpendicular to the view
   // direction at pos. If the field of view angle is small then this difference
@@ -179,7 +176,6 @@ inline bool igl::RotateWidget::intersect(
   Eigen::Vector3d & hit) const
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d view = view_direction(x,y);
   const Vector3d ppos = project(pos);
   Vector3d uxy = unproject(Vector3d(x,y,ppos(2)));
@@ -196,7 +192,6 @@ inline bool igl::RotateWidget::intersect(
 inline double igl::RotateWidget::unprojected_inner_radius() const
 {
   using namespace Eigen;
-  using namespace igl;
   Vector3d off,ppos,ppos_off,pos_off;
   project(pos,ppos);
   ppos_off = ppos;
@@ -207,7 +202,6 @@ inline double igl::RotateWidget::unprojected_inner_radius() const
 inline bool igl::RotateWidget::down(const int x, const int y)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   if(!m_is_enabled)
   {
@@ -289,7 +283,6 @@ inline bool igl::RotateWidget::down(const int x, const int y)
 
 inline bool igl::RotateWidget::drag(const int x, const int y)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   if(!m_is_enabled)
@@ -378,7 +371,6 @@ inline void igl::RotateWidget::draw() const
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   glPushAttrib(GL_ENABLE_BIT | GL_LIGHTING_BIT | GL_DEPTH_BUFFER_BIT | GL_LINE_BIT);
   glDisable(GL_CLIP_PLANE0);
 
@@ -474,7 +466,6 @@ inline void igl::RotateWidget::draw_guide() const
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   glPushAttrib(
     GL_DEPTH_BUFFER_BIT | 
     GL_ENABLE_BIT | 

include/igl/WindingNumberAABB.h

@@ -83,7 +83,6 @@ inline void igl::WindingNumberAABB<Point>::set_mesh(
 template <typename Point>
 inline void igl::WindingNumberAABB<Point>::init()
 {
-  using namespace igl;
   using namespace Eigen;
   assert(max_corner.size() == 3);
   assert(min_corner.size() == 3);
@@ -124,7 +123,6 @@ inline void igl::WindingNumberAABB<Point>::grow()
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   //cout<<"cap.rows(): "<<this->getcap().rows()<<endl;
   //cout<<"F.rows(): "<<this->getF().rows()<<endl;
 
@@ -300,7 +298,6 @@ inline double igl::WindingNumberAABB<Point>::max_simple_abs_winding_number(const
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // Only valid if not inside
   if(inside(p))
   {

include/igl/WindingNumberTree.h

@@ -406,7 +406,6 @@ inline double igl::WindingNumberTree<Point>::cached_winding_number(
   const Point & p) const
 {
   using namespace std;
-  using namespace igl;
   // Simple metric for "far".
   //   this             that
   //                   --------

include/igl/active_set.cpp

@@ -48,7 +48,6 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #ifdef ACTIVE_SET_CPP_DEBUG
 #  warning "ACTIVE_SET_CPP_DEBUG"
 #endif
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   SolverStatus ret = SOLVER_STATUS_ERROR;

include/igl/all_edges.cpp

@@ -49,4 +49,5 @@ template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matri
 template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::all_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
 #endif

include/igl/angular_distance.cpp

@@ -12,7 +12,6 @@ IGL_INLINE double igl::angular_distance(
   const Eigen::Quaterniond & A,
   const Eigen::Quaterniond & B)
 {
-  using namespace igl;
   assert(fabs(A.norm()-1)<FLOAT_EPS && "A should be unit norm");
   assert(fabs(B.norm()-1)<FLOAT_EPS && "B should be unit norm");
   //// acos is always in [0,2*pi)

include/igl/arap_linear_block.cpp

@@ -18,7 +18,6 @@ IGL_INLINE void igl::arap_linear_block(
   const igl::ARAPEnergyType energy,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   switch(energy)
   {
     case ARAP_ENERGY_TYPE_SPOKES:
@@ -44,7 +43,6 @@ IGL_INLINE void igl::arap_linear_block_spokes(
   const int d,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)
@@ -110,7 +108,6 @@ IGL_INLINE void igl::arap_linear_block_spokes_and_rims(
   const int d,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)
@@ -193,7 +190,6 @@ IGL_INLINE void igl::arap_linear_block_elements(
   const int d,
   Eigen::SparseMatrix<Scalar> & Kd)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)

include/igl/arap_rhs.cpp

@@ -19,7 +19,6 @@ IGL_INLINE void igl::arap_rhs(
   const igl::ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& K)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // Number of dimensions

include/igl/boolean/mesh_boolean.cpp

@@ -169,7 +169,7 @@ IGL_INLINE void igl::mesh_boolean(
     return;
   }
   MatrixX3S N,CN;
-  per_face_normals(V,F,N);
+  per_face_normals_stable(V,F,N);
   CN.resize(CF.rows(),3);
   for(size_t f = 0;f<(size_t)CN.rows();f++)
   {

include/igl/boundary_conditions.cpp

@@ -26,7 +26,6 @@ IGL_INLINE bool igl::boundary_conditions(
   Eigen::MatrixXd &       bc )
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
   if(P.size()+BE.rows() == 0)

include/igl/boundary_facets.cpp

@@ -21,7 +21,6 @@ IGL_INLINE void igl::boundary_facets(
   std::vector<std::vector<IntegerF> > & F)
 {
   using namespace std;
-  using namespace igl;
 
   if(T.size() == 0)
   {
@@ -111,7 +110,6 @@ IGL_INLINE void igl::boundary_facets(
   assert(T.cols() == 0 || T.cols() == 4 || T.cols() == 3);
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // Cop out: use vector of vectors version
   vector<vector<typename Eigen::PlainObjectBase<DerivedT>::Scalar> > vT;
   matrix_to_list(T,vT);

include/igl/boundary_loop.cpp

@@ -19,7 +19,6 @@ IGL_INLINE void igl::boundary_loop(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
 
   if(F.rows() == 0)
     return;

include/igl/bounding_box_diagonal.cpp

@@ -13,7 +13,6 @@
 IGL_INLINE double igl::bounding_box_diagonal(
   const Eigen::MatrixXd & V)
 {
-  using namespace igl;
   using namespace Eigen;
   VectorXd maxV,minV;
   VectorXi maxVI,minVI;

include/igl/cat.cpp

@@ -123,9 +123,8 @@ IGL_INLINE Mat igl::cat(const int dim, const Mat & A, const Mat & B)
 }
 
 template <class Mat>
-IGL_INLINE void cat(const std::vector<std::vector< Mat > > & A, Mat & C)
+IGL_INLINE void igl::cat(const std::vector<std::vector< Mat > > & A, Mat & C)
 {
-  using namespace igl;
   using namespace std;
   // Start with empty matrix
   C.resize(0,0);

include/igl/cgal/CGAL_includes.hpp

@@ -8,6 +8,8 @@
 #ifndef IGL_CGAL_INCLUDES_H
 #define IGL_CGAL_INCLUDES_H
 
+// http://www.alecjacobson.com/weblog/?p=4291
+#define CGAL_INTERSECTION_VERSION 1
 // Triangle triangle intersection
 #include <CGAL/intersections.h>
 // THIS CANNOT BE INCLUDED IN THE SAME FILE AS <CGAL/intersections.h>

include/igl/cgal/RemeshSelfIntersectionsParam.h

@@ -0,0 +1,23 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_REMESH_SELF_INTERSECTIONS_PARAM_H
+#define IGL_REMESH_SELF_INTERSECTIONS_PARAM_H
+
+namespace igl
+{
+  // Optional Parameters
+  //   DetectOnly  Only compute IF, leave VV and FF alone
+  struct RemeshSelfIntersectionsParam
+  {
+    bool detect_only;
+    bool first_only;
+    RemeshSelfIntersectionsParam():detect_only(false),first_only(false){};
+  };
+}
+
+#endif

include/igl/cgal/SelfIntersectMesh.h

@@ -9,13 +9,14 @@
 #define IGL_SELFINTERSECTMESH_H
 
 #include "CGAL_includes.hpp"
-#include "remesh_self_intersections.h"
+#include "RemeshSelfIntersectionsParam.h"
 
 #include <Eigen/Dense>
 #include <list>
 #include <map>
 #include <vector>
 
+//#define IGL_SELFINTERSECTMESH_DEBUG
 #ifndef IGL_FIRST_HIT_EXCEPTION
 #define IGL_FIRST_HIT_EXCEPTION 10
 #endif
@@ -325,18 +326,22 @@ inline igl::SelfIntersectMesh<
   using namespace std;
   using namespace Eigen;
 
-  //const auto & tictoc = []()
-  //{
-  //  static double t_start = igl::get_seconds();
-  //  double diff = igl::get_seconds()-t_start;
-  //  t_start += diff;
-  //  return diff;
-  //};
-  //tictoc();
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  const auto & tictoc = []()
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  tictoc();
+#endif
 
   // Compute and process self intersections
   mesh_to_cgal_triangle_list(V,F,T);
-  //cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
+#endif
   // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
   // Create the corresponding vector of bounding boxes
   std::vector<Box> boxes;
@@ -355,7 +360,9 @@ inline igl::SelfIntersectMesh<
       // _1 etc. in global namespace)
       std::placeholders::_1,
       std::placeholders::_2);
-  //cout<<"boxes and bind: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"boxes and bind: "<<tictoc()<<endl;
+#endif
   // Run the self intersection algorithm with all defaults
   try{
     CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
@@ -368,7 +375,9 @@ inline igl::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
-  //cout<<"box_self_intersection_d: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"box_self_intersection_d: "<<tictoc()<<endl;
+#endif
 
   // Convert lIF to Eigen matrix
   assert(lIF.size()%2 == 0);
@@ -387,7 +396,9 @@ inline igl::SelfIntersectMesh<
       i++;
     }
   }
-  //cout<<"IF: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"IF: "<<tictoc()<<endl;
+#endif
 
   if(params.detect_only)
   {
@@ -409,6 +420,9 @@ inline igl::SelfIntersectMesh<
   map<typename CDT_plus_2::Vertex_handle,Index> v2i;
   // Loop over offending triangles
   const size_t noff = offending.size();
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  double t_proj_del = 0;
+#endif
   // Unfortunately it looks like CGAL has trouble allocating memory when
   // multiple openmp threads are running. Crashes durring CDT...
 //# pragma omp parallel for if (noff>1000)
@@ -417,7 +431,13 @@ inline igl::SelfIntersectMesh<
     // index in F
     const Index f = offending[o];
     {
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+      const double t_before = get_seconds();
+#endif
       projected_delaunay(T[f],F_objects[f],cdt[o]);
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+      t_proj_del += (get_seconds()-t_before);
+#endif
     }
     // Q: Is this also delaunay in 3D?
     // A: No, because the projection is affine and delaunay is not affine
@@ -527,7 +547,10 @@ inline igl::SelfIntersectMesh<
       }
     }
   }
-  //cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
+#endif
+
   assert(NV_count == (Index)NV.size());
   // Build output
 #ifndef NDEBUG
@@ -565,7 +588,7 @@ inline igl::SelfIntersectMesh<
   }
   // Append vertices
   VV.resize(V.rows()+NV_count,3);
-  VV.block(0,0,V.rows(),3) = V;
+  VV.block(0,0,V.rows(),3) = V.template cast<typename DerivedVV::Scalar>();
   {
     Index i = 0;
     for(
@@ -576,13 +599,14 @@ inline igl::SelfIntersectMesh<
       for(Index d = 0;d<3;d++)
       {
         const Point_3 & p = *nvit;
-        VV(V.rows()+i,d) = CGAL::to_double(p[d]);
-        // This distinction does not seem necessary:
-//#ifdef INEXACT_CONSTRUCTION
-//        VV(V.rows()+i,d) = CGAL::to_double(p[d]);
-//#else
-//        VV(V.rows()+i,d) = CGAL::to_double(p[d].exact());
-//#endif
+        // Don't convert via double if output type is same as Kernel
+        if(std::is_same<typename DerivedVV::Scalar,typename Kernel::FT>::value)
+        {
+          VV(V.rows()+i,d) = p[d];
+        }else
+        {
+          VV(V.rows()+i,d) = CGAL::to_double(p[d]);
+        }
       }
       i++;
     }
@@ -619,7 +643,9 @@ inline igl::SelfIntersectMesh<
       v++;
     }
   }
-  //cout<<"Output + dupes: "<<tictoc()<<endl;
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"Output + dupes: "<<tictoc()<<endl;
+#endif
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.

include/igl/cgal/mesh_to_cgal_triangle_list.cpp

@@ -9,10 +9,13 @@
 
 #include <cassert>
 
-template <typename Kernel>
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Kernel>
 IGL_INLINE void igl::mesh_to_cgal_triangle_list(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   std::vector<CGAL::Triangle_3<Kernel> > & T)
 {
   typedef CGAL::Point_3<Kernel>    Point_3;
@@ -35,7 +38,4 @@ IGL_INLINE void igl::mesh_to_cgal_triangle_list(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::mesh_to_cgal_triangle_list<CGAL::Epeck>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);
-template void igl::mesh_to_cgal_triangle_list<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-template void igl::mesh_to_cgal_triangle_list<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Simple_cartesian<double> >, std::allocator<CGAL::Triangle_3<CGAL::Simple_cartesian<double> > > >&);
 #endif

include/igl/cgal/mesh_to_cgal_triangle_list.h

@@ -22,10 +22,13 @@ namespace igl
   //   F  #F by 3 list of triangle indices
   // Outputs:
   //   T  #F list of CGAL triangles
-  template <typename Kernel>
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Kernel>
   IGL_INLINE void mesh_to_cgal_triangle_list(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
     std::vector<CGAL::Triangle_3<Kernel> > & T);
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/cgal/remesh_self_intersections.h

@@ -8,6 +8,7 @@
 #ifndef IGL_REMESH_SELF_INTERSECTIONS_H
 #define IGL_REMESH_SELF_INTERSECTIONS_H
 #include <igl/igl_inline.h>
+#include "RemeshSelfIntersectionsParam.h"
 
 #include <Eigen/Dense>
 
@@ -17,18 +18,9 @@
 #  undef assert
 #  define assert( isOK ) ( (isOK) ? (void)0 : (void) mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
 #endif
-
+  
 namespace igl
 {
-  // Optional Parameters
-  //   DetectOnly  Only compute IF, leave VV and FF alone
-  struct RemeshSelfIntersectionsParam
-  {
-    bool detect_only;
-    bool first_only;
-    RemeshSelfIntersectionsParam():detect_only(false),first_only(false){};
-  };
-  
   // Given a triangle mesh (V,F) compute a new mesh (VV,FF) which is the same
   // as (V,F) except that any self-intersecting triangles in (V,F) have been
   // subdivided (new vertices and face created) so that the self-intersection

include/igl/collapse_edge.cpp

@@ -29,7 +29,6 @@ IGL_INLINE bool igl::collapse_edge(
   // never get collapsed to anything else since it is the smallest index)
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   const int eflip = E(e,0)>E(e,1);
   // source and destination
   const int s = eflip?E(e,1):E(e,0);
@@ -154,3 +153,117 @@ IGL_INLINE bool igl::collapse_edge(
   int e1,e2,f1,f2;
   return collapse_edge(e,p,V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
 }
+
+IGL_INLINE bool igl::collapse_edge(
+  const std::function<void(
+    const int,
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    double &,
+    Eigen::RowVectorXd &)> & cost_and_placement,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI,
+  std::set<std::pair<double,int> > & Q,
+  std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+  Eigen::MatrixXd & C)
+{
+  int e,e1,e2,f1,f2;
+  return 
+    collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+}
+
+
+IGL_INLINE bool igl::collapse_edge(
+  const std::function<void(
+    const int,
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    double &,
+    Eigen::RowVectorXd &)> & cost_and_placement,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI,
+  std::set<std::pair<double,int> > & Q,
+  std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+  Eigen::MatrixXd & C,
+  int & e,
+  int & e1,
+  int & e2,
+  int & f1,
+  int & f2)
+{
+  using namespace Eigen;
+  if(Q.empty())
+  {
+    // no edges to collapse
+    return false;
+  }
+  std::pair<double,int> p = *(Q.begin());
+  if(p.first == std::numeric_limits<double>::infinity())
+  {
+    // min cost edge is infinite cost
+    return false;
+  }
+  Q.erase(Q.begin());
+  e = p.second;
+  Qit[e] = Q.end();
+  std::vector<int> N  = circulation(e, true,F,E,EMAP,EF,EI);
+  std::vector<int> Nd = circulation(e,false,F,E,EMAP,EF,EI);
+  N.insert(N.begin(),Nd.begin(),Nd.end());
+  const bool collapsed =
+    collapse_edge(e,C.row(e),V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+  if(collapsed)
+  {
+    // Erase the two, other collapsed edges
+    Q.erase(Qit[e1]);
+    Qit[e1] = Q.end();
+    Q.erase(Qit[e2]);
+    Qit[e2] = Q.end();
+    // update local neighbors
+    // loop over original face neighbors
+    for(auto n : N)
+    {
+      if(F(n,0) != IGL_COLLAPSE_EDGE_NULL ||
+          F(n,1) != IGL_COLLAPSE_EDGE_NULL ||
+          F(n,2) != IGL_COLLAPSE_EDGE_NULL)
+      {
+        for(int v = 0;v<3;v++)
+        {
+          // get edge id
+          const int ei = EMAP(v*F.rows()+n);
+          // erase old entry
+          Q.erase(Qit[ei]);
+          // compute cost and potential placement
+          double cost;
+          RowVectorXd place;
+          cost_and_placement(ei,V,F,E,EMAP,EF,EI,cost,place);
+          // Replace in queue
+          Qit[ei] = Q.insert(std::pair<double,int>(cost,ei)).first;
+          C.row(ei) = place;
+        }
+      }
+    }
+  }else
+  {
+    // reinsert with infinite weight (the provided cost function must **not**
+    // have given this un-collapsable edge inf cost already)
+    p.first = std::numeric_limits<double>::infinity();
+    Qit[e] = Q.insert(p).first;
+  }
+  return collapsed;
+}

include/igl/collapse_edge.h

@@ -9,6 +9,8 @@
 #define IGL_COLLAPSE_EDGE_H
 #include "igl_inline.h"
 #include <Eigen/Core>
+#include <vector>
+#include <set>
 namespace igl
 {
   // Assumes (V,F) is a closed manifold mesh (except for previouslly collapsed
@@ -62,6 +64,60 @@ namespace igl
     Eigen::VectorXi & EMAP,
     Eigen::MatrixXi & EF,
     Eigen::MatrixXi & EI);
+  // Collapse least-cost edge from a priority queue and update queue 
+  //
+  // Inputs/Outputs:
+  //   cost_and_placement  function computing cost of collapsing an edge and 3d
+  //     position where it should be placed:
+  //     cost_and_placement(V,F,E,EMAP,EF,EI,cost,placement);
+  //   Q  queue containing pairs of costs and edge indices
+  //   Qit  list of iterators so that Qit[e] --> iterator of edge e in Q
+  //   C  #E by dim list of stored placements
+  IGL_INLINE bool collapse_edge(
+    const std::function<void(
+      const int,
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      double &,
+      Eigen::RowVectorXd &)> & cost_and_placement,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI,
+    std::set<std::pair<double,int> > & Q,
+    std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    Eigen::MatrixXd & C);
+  IGL_INLINE bool collapse_edge(
+    const std::function<void(
+      const int,
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      double &,
+      Eigen::RowVectorXd &)> & cost_and_placement,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI,
+    std::set<std::pair<double,int> > & Q,
+    std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    Eigen::MatrixXd & C,
+    int & e,
+    int & e1,
+    int & e2,
+    int & f1,
+    int & f2);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/compile_and_link_program.cpp

@@ -0,0 +1,33 @@
+#include "compile_and_link_program.h"
+#include "compile_shader.h"
+#include "report_gl_error.h"
+#include <iostream>
+#include <cassert>
+
+IGL_INLINE GLuint igl::compile_and_link_program(
+  const char * v_str, const char * f_str)
+{
+  GLuint vid = compile_shader(GL_VERTEX_SHADER,v_str);
+  GLuint fid = compile_shader(GL_FRAGMENT_SHADER,f_str);
+
+  GLuint prog_id = glCreateProgram();
+  assert(prog_id != 0 && "Failed to create shader.");
+  glAttachShader(prog_id,vid);
+  igl::report_gl_error("glAttachShader (vid): ");
+  glAttachShader(prog_id,fid);
+  igl::report_gl_error("glAttachShader (fid): ");
+
+  glLinkProgram(prog_id);
+  igl::report_gl_error("glLinkProgram: ");
+
+  GLint status;
+  glGetProgramiv(prog_id, GL_LINK_STATUS, &status);
+  if (status != GL_TRUE)
+  {
+    char buffer[512];
+    glGetProgramInfoLog(prog_id, 512, NULL, buffer);
+    std::cerr << "Linker error: " << std::endl << buffer << std::endl;
+    prog_id = 0;
+  }
+  return prog_id;
+}

include/igl/compile_and_link_program.h

@@ -0,0 +1,22 @@
+#ifndef IGL_COMPILE_AND_LINK_PROGRAM_H
+#define IGL_COMPILE_AND_LINK_PROGRAM_H
+#include "igl_inline.h"
+#include "OpenGL_convenience.h"
+namespace igl
+{
+  // Compile and link very simple vertex/fragment shaders
+  //
+  // Inputs:
+  //   v_str  string of vertex shader contents
+  //   f_str  string of fragment shader contents
+  // Returns id of program
+  //
+  // Known bugs: this seems to duplicate `create_shader_program` with less
+  // functionality.
+  IGL_INLINE GLuint compile_and_link_program(
+    const char * v_str, const char * f_str);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "compile_and_link_program.cpp"
+#endif
+#endif

include/igl/compile_shader.cpp

@@ -0,0 +1,28 @@
+#include "compile_shader.h"
+#include "report_gl_error.h"
+#include <iostream>
+
+IGL_INLINE GLuint igl::compile_shader(const GLint type, const char * str)
+{
+  GLuint id = glCreateShader(type);
+  igl::report_gl_error("glCreateShader: ");
+  glShaderSource(id,1,&str,NULL);
+  igl::report_gl_error("glShaderSource: ");
+  glCompileShader(id);
+  igl::report_gl_error("glCompileShader: ");
+
+  GLint status;
+  glGetShaderiv(id, GL_COMPILE_STATUS, &status);
+  if (status != GL_TRUE)
+  {
+    char buffer[512];
+    if (type == GL_VERTEX_SHADER)
+      std::cerr << "Vertex shader:" << std::endl;
+    else if (type == GL_FRAGMENT_SHADER)
+      std::cerr << "Fragment shader:" << std::endl;
+    std::cerr << str << std::endl << std::endl;
+    glGetShaderInfoLog(id, 512, NULL, buffer);
+    std::cerr << "Error: " << std::endl << buffer << std::endl;
+  }
+  return id;
+}

include/igl/compile_shader.h

@@ -0,0 +1,26 @@
+#ifndef IGL_COMPILE_SHADER_H
+#define IGL_COMPILE_SHADER_H
+#include "OpenGL_convenience.h"
+#include "igl_inline.h"
+namespace igl
+{
+  // Compile a shader given type and string of shader code
+  // 
+  // Inputs:
+  //   type  either GL_VERTEX_SHADER or GL_FRAGMENT_SHADER
+  //   str  contents of shader code
+  // Returns result of glCreateShader (id of shader)
+  //
+  // Example:
+  //     GLuint vid = compile_shader(GL_VERTEX_SHADER,vertex_shader.c_str());
+  //     GLuint fid = compile_shader(GL_FRAGMENT_SHADER,fragment_shader.c_str());
+  //     GLuint prog_id = glCreateProgram();
+  //     glAttachShader(prog_id,vid);
+  //     glAttachShader(prog_id,fid);
+  //     glLinkProgram(prog_id);
+  IGL_INLINE GLuint compile_shader(const GLint type, const char * str);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "compile_shader.cpp"
+#endif
+#endif

include/igl/cotmatrix.cpp

@@ -22,7 +22,6 @@ IGL_INLINE void igl::cotmatrix(
   const Eigen::PlainObjectBase<DerivedF> & F, 
   Eigen::SparseMatrix<Scalar>& L)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
 

include/igl/cotmatrix_entries.cpp

@@ -21,7 +21,6 @@ IGL_INLINE void igl::cotmatrix_entries(
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedC>& C)
 {
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
   // simplex size (3: triangles, 4: tetrahedra)

include/igl/covariance_scatter_matrix.cpp

@@ -20,7 +20,6 @@ IGL_INLINE void igl::covariance_scatter_matrix(
   const ARAPEnergyType energy,
   Eigen::SparseMatrix<double>& CSM)
 {
-  using namespace igl;
   using namespace Eigen;
   // number of mesh vertices
   int n = V.rows();

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -0,0 +1,68 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "crouzeix_raviart_massmatrix.h"
+#include "unique_simplices.h"
+#include "all_edges.h"
+
+#include "is_edge_manifold.h"
+#include "doublearea.h"
+
+#include <cassert>
+#include <vector>
+
+template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+void igl::crouzeix_raviart_massmatrix(
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    Eigen::SparseMatrix<MT> & M,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  // All occurances of directed edges
+  Eigen::MatrixXi allE;
+  all_edges(F,allE);
+  Eigen::VectorXi _1;
+  unique_simplices(allE,E,_1,EMAP);
+  return crouzeix_raviart_massmatrix(V,F,E,EMAP,M);
+}
+
+template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+void igl::crouzeix_raviart_massmatrix(
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
+    const Eigen::PlainObjectBase<DerivedE> & E,
+    const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    Eigen::SparseMatrix<MT> & M)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(F.cols() == 3);
+  // Mesh should be edge-manifold
+  assert(is_edge_manifold(V,F));
+  // number of elements (triangles)
+  int m = F.rows();
+  // Get triangle areas
+  VectorXd TA;
+  doublearea(V,F,TA);
+  vector<Triplet<MT> > MIJV(3*m);
+  for(int f = 0;f<m;f++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      MIJV[f+m*c] = Triplet<MT>(EMAP(f+m*c),EMAP(f+m*c),TA(f)*0.5);
+    }
+  }
+  M.resize(E.rows(),E.rows());
+  M.setFromTriplets(MIJV.begin(),MIJV.end());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit instanciation
+template void igl::crouzeix_raviart_massmatrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::crouzeix_raviart_massmatrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<float, 0, int>&);
+#endif

include/igl/crouzeix_raviart_massmatrix.h

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef CROUZEIX_RAVIART_MASSMATRIX_H
+#define CROUZEIX_RAVIART_MASSMATRIX_H
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // CROUZEIX_RAVIART_MASSMATRIX Compute the Crouzeix-Raviart mass matrix where
+  // M(e,e) is just the sum of the areas of the triangles on either side of an
+  // edge e.
+  //
+  // See for example "Discrete Quadratic Curvature Energies" [Wardetzky, Bergou,
+  // Harmon, Zorin, Grinspun 2007]
+  //
+  // Templates:
+  //   MT  type of eigen sparse matrix for M (e.g. double for
+  //     SparseMatrix<double>)
+  //   DerivedV  derived type of eigen matrix for V (e.g. derived from
+  //     MatrixXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  //   DerivedE  derived type of eigen matrix for E (e.g. derived from
+  //     MatrixXi)
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   M  #E by #E edge-based diagonal mass matrix
+  //   E  #E by 2 list of edges
+  //   EMAP  #F*3 list of indices mapping allE to E
+  //
+  //
+  template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+  void crouzeix_raviart_massmatrix(
+      const Eigen::PlainObjectBase<DerivedV> & V, 
+      const Eigen::PlainObjectBase<DerivedF> & F, 
+      Eigen::SparseMatrix<MT> & M,
+      Eigen::PlainObjectBase<DerivedE> & E,
+      Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+  // wrapper if E and EMAP are already computed (better match!)
+  template <typename MT, typename DerivedV, typename DerivedF, typename DerivedE, typename DerivedEMAP>
+  void crouzeix_raviart_massmatrix(
+      const Eigen::PlainObjectBase<DerivedV> & V, 
+      const Eigen::PlainObjectBase<DerivedF> & F, 
+      const Eigen::PlainObjectBase<DerivedE> & E,
+      const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+      Eigen::SparseMatrix<MT> & M);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "crouzeix_raviart_massmatrix.cpp"
+#endif
+  
+#endif

include/igl/dated_copy.cpp

@@ -19,7 +19,6 @@
 IGL_INLINE bool igl::dated_copy(const std::string & src_path, const std::string & dir)
 {
   using namespace std;
-  using namespace igl;
   // Get time and date as string
   char buffer[80];
   time_t rawtime;

include/igl/decimate.cpp

@@ -0,0 +1,164 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "decimate.h"
+#include "collapse_edge.h"
+#include "edge_flaps.h"
+#include "remove_unreferenced.h"
+#include <iostream>
+
+IGL_INLINE bool igl::decimate(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const size_t max_m,
+  Eigen::MatrixXd & U,
+  Eigen::MatrixXi & G)
+{
+  const auto & max_m_faces = 
+    [&max_m,&F](
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const std::set<std::pair<double,int> > &,
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,
+    const Eigen::MatrixXd &,
+    const int,
+    const int,
+    const int,
+    const int,
+    const int)->bool
+    {
+      using namespace std;
+      static int m = F.rows();
+      m-=2;
+      return m<=(int)max_m;
+    };
+  const auto & shortest_edge_and_midpoint = [](
+    const int e,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & /*F*/,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & /*EMAP*/,
+    const Eigen::MatrixXi & /*EF*/,
+    const Eigen::MatrixXi & /*EI*/,
+    double & cost,
+    Eigen::RowVectorXd & p)
+  {
+    cost = (V.row(E(e,0))-V.row(E(e,1))).norm();
+    p = 0.5*(V.row(E(e,0))+V.row(E(e,1)));
+  };
+  return decimate(V,F,shortest_edge_and_midpoint,max_m_faces,U,G);
+}
+
+IGL_INLINE bool igl::decimate(
+  const Eigen::MatrixXd & OV,
+  const Eigen::MatrixXi & OF,
+  const std::function<void(
+    const int,
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    double &,
+    Eigen::RowVectorXd &)> & cost_and_placement,
+  const std::function<bool(
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const std::set<std::pair<double,int> > &,
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,
+      const Eigen::MatrixXd &,
+      const int,
+      const int,
+      const int,
+      const int,
+      const int)> & stopping_condition,
+  Eigen::MatrixXd & U,
+  Eigen::MatrixXi & G)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Working copies
+  Eigen::MatrixXd V = OV;
+  Eigen::MatrixXi F = OF;
+  VectorXi EMAP;
+  MatrixXi E,EF,EI;
+  typedef std::set<std::pair<double,int> > PriorityQueue;
+  PriorityQueue Q;
+  std::vector<PriorityQueue::iterator > Qit;
+  edge_flaps(F,E,EMAP,EF,EI);
+  Qit.resize(E.rows());
+  // If an edge were collapsed, we'd collapse it to these points:
+  MatrixXd C(E.rows(),V.cols());
+  for(int e = 0;e<E.rows();e++)
+  {
+    double cost = e;
+    RowVectorXd p(1,3);
+    cost_and_placement(e,V,F,E,EMAP,EF,EI,cost,p);
+    C.row(e) = p;
+    Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
+  }
+
+  int prev_e = -1;
+  bool clean_finish = false;
+  while(true)
+  {
+    if(Q.empty())
+    {
+      break;
+    }
+    if(Q.begin()->first == std::numeric_limits<double>::infinity())
+    {
+      // min cost edge is infinite cost
+      break;
+    }
+    int e,e1,e2,f1,f2;
+    if(collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2))
+    {
+      if(stopping_condition(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2))
+      {
+        clean_finish = true;
+        break;
+      }
+    }else
+    {
+      if(prev_e == e)
+      {
+        assert(false && "Edge collapse no progress... bad stopping condition?");
+        break;
+      }
+      // Edge was not collapsed... must have been invalid. collapse_edge should
+      // have updated its cost to inf... continue
+    }
+    prev_e = e;
+  }
+  // remove all IGL_COLLAPSE_EDGE_NULL faces
+  MatrixXi F2(F.rows(),3);
+  int m = 0;
+  for(int f = 0;f<F.rows();f++)
+  {
+    if(
+      F(f,0) != IGL_COLLAPSE_EDGE_NULL || 
+      F(f,1) != IGL_COLLAPSE_EDGE_NULL || 
+      F(f,2) != IGL_COLLAPSE_EDGE_NULL)
+    {
+      F2.row(m++) = F.row(f);
+    }
+  }
+  F2.conservativeResize(m,F2.cols());
+  VectorXi _1;
+  remove_unreferenced(V,F2,U,G,_1);
+  return clean_finish;
+}

include/igl/decimate.h

@@ -0,0 +1,81 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_DECIMATE_H
+#define IGL_DECIMATE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+#include <set>
+namespace igl
+{
+  // Assumes (V,F) is a closed manifold mesh 
+  // Collapses edges until desired number of faces is achieved. This uses
+  // default edge cost and merged vertex placement functions {edge length, edge
+  // midpoint}.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3 list of face indices into V.
+  //   max_m  desired number of output faces
+  // Outputs:
+  //   U  #U by dim list of output vertex posistions (can be same ref as V)
+  //   G  #G by 3 list of output face indices into U (can be same ref as G)
+  // Returns true if m was reached (otherwise #G > m)
+  IGL_INLINE bool decimate(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const size_t max_m,
+    Eigen::MatrixXd & U,
+    Eigen::MatrixXi & G);
+  // Inputs:
+  //   cost_and_placement  function computing cost of collapsing an edge and 3d
+  //     position where it should be placed:
+  //     cost_and_placement(V,F,E,EMAP,EF,EI,cost,placement);
+  //   stopping_condition  function returning whether to stop collapsing edges
+  //     based on current state. Guaranteed to be called after _successfully_
+  //     collapsing edge e removing edges (e,e1,e2) and faces (f1,f2):
+  //     bool should_stop =
+  //       stopping_condition(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+  IGL_INLINE bool decimate(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const std::function<void(
+      const int,
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      double &,
+      Eigen::RowVectorXd &)> & cost_and_placement,
+    const std::function<bool(
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const std::set<std::pair<double,int> > &,
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,
+      const Eigen::MatrixXd &,
+      const int,
+      const int,
+      const int,
+      const int,
+      const int)> & stopping_condition,
+    Eigen::MatrixXd & U,
+    Eigen::MatrixXi & G);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "decimate.cpp"
+#endif
+#endif
+

include/igl/dqs.cpp

@@ -22,7 +22,6 @@ IGL_INLINE void igl::dqs(
   Eigen::PlainObjectBase<DerivedU> & U)
 {
   using namespace std;
-  using namespace igl;
   assert(V.rows() <= W.rows());
   assert(W.cols() == (int)vQ.size());
   assert(W.cols() == (int)vT.size());

include/igl/draw_mesh.h

@@ -20,7 +20,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V|#F by 3 eigen Matrix of 3D normals
   IGL_INLINE void draw_mesh(
     const Eigen::MatrixXd & V,
@@ -32,7 +32,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V|#F by 3 eigen Matrix of 3D normals
   //   C  #V|#F|1 by 3 eigen Matrix of RGB colors
   IGL_INLINE void draw_mesh(
@@ -42,7 +42,7 @@ namespace igl
     const Eigen::MatrixXd & C);
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V|#F by 3 eigen Matrix of 3D normals
   //   C  #V|#F|1 by 3 eigen Matrix of RGB colors
   //   TC  #V|#F|1 by 3 eigen Matrix of Texture Coordinates
@@ -58,7 +58,7 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
   //   C  #V by 3 eigen Matrix of mesh vertex RGB colors
   //   TC  #V by 3 eigen Matrix of mesh vertex UC coorindates between 0 and 1
@@ -84,14 +84,14 @@ namespace igl
   //
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
-  //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   F  #F by 3|4 eigen Matrix of face (triangle/quad) indices
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
-  //   NF  #F by 3 eigen Matrix of face (triangle) normal indices, <0 means no
-  //     normal
+  //   NF  #F by 3 eigen Matrix of face (triangle/quad) normal indices, <0
+  //     means no normal
   //   C  #V by 3 eigen Matrix of mesh vertex RGB colors
   //   TC  #V by 3 eigen Matrix of mesh vertex UC coorindates between 0 and 1
-  //   TF  #F by 3 eigen Matrix of face (triangle) texture indices, <0 means no
-  //     texture
+  //   TF  #F by 3 eigen Matrix of face (triangle/quad) texture indices, <0
+  //     means no texture
   //   W  #V by #H eigen Matrix of per mesh vertex, per handle weights
   //   W_index  Specifies the index of the "weight" vertex attribute: see
   //     glBindAttribLocation, if W_index is 0 then weights are ignored

include/igl/draw_rectangular_marquee.cpp

@@ -15,7 +15,6 @@ IGL_INLINE void igl::draw_rectangular_marquee(
   const int to_x,
   const int to_y)
 {
-  using namespace igl;
   using namespace std;
   int l;
   glGetIntegerv(GL_LIGHTING,&l);

include/igl/edge_collapse_is_valid.cpp

@@ -23,7 +23,6 @@ IGL_INLINE bool igl::edge_collapse_is_valid(
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   // For consistency with collapse_edge.cpp, let's determine edge flipness
   // (though not needed to check validity)
   const int eflip = E(e,0)>E(e,1);

include/igl/facet_components.cpp

@@ -8,7 +8,6 @@ IGL_INLINE void igl::facet_components(
   Eigen::PlainObjectBase<DerivedC> & C)
 {
   using namespace std;
-  using namespace igl;
   typedef typename DerivedF::Index Index;
   vector<vector<vector<Index > > > TT;
   vector<vector<vector<Index > > > TTi;
@@ -27,7 +26,6 @@ IGL_INLINE void igl::facet_components(
   Eigen::PlainObjectBase<Derivedcounts> & counts)
 {
   using namespace std;
-  using namespace igl;
   typedef TTIndex Index;
   const Index m = TT.size();
   C.resize(m,1);

include/igl/harmonic.cpp

@@ -26,7 +26,6 @@ IGL_INLINE bool igl::harmonic(
   const int k,
   Eigen::PlainObjectBase<DerivedW> & W)
 {
-  using namespace igl;
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
   typedef Matrix<Scalar,Dynamic,1> VectorXS;

include/igl/in_element.cpp

@@ -9,7 +9,6 @@ IGL_INLINE void igl::in_element(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   const int Qr = Q.rows();
   I.setConstant(Qr,1,-1);
 #pragma omp parallel for

include/igl/is_boundary_edge.cpp

@@ -0,0 +1,122 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "is_boundary_edge.h"
+#include "unique.h"
+#include "sort.h"
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DerivedB>
+void igl::is_boundary_edge(
+  const Eigen::PlainObjectBase<DerivedE> & E,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedB> & B)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Should be triangles
+  assert(F.cols() == 3);
+  // Should be edges
+  assert(E.cols() == 2);
+  // number of faces
+  const int m = F.rows();
+  // Collect all directed edges after E
+  MatrixXi EallE(E.rows()+3*m,2);
+  EallE.block(0,0,E.rows(),E.cols()) = E;
+  for(int e = 0;e<3;e++)
+  {
+    for(int f = 0;f<m;f++)
+    {
+      for(int c = 0;c<2;c++)
+      {
+        // 12 20 01
+        EallE(E.rows()+m*e+f,c) = F(f,(c+1+e)%3);
+      }
+    }
+  }
+  // sort directed edges into undirected edges
+  MatrixXi sorted_EallE,_;
+  sort(EallE,2,true,sorted_EallE,_);
+  // Determine unique undirected edges E and map to directed edges EMAP
+  MatrixXi uE;
+  VectorXi EMAP;
+  unique_rows(sorted_EallE,uE,_,EMAP);
+  // Counts of occurances
+  VectorXi N = VectorXi::Zero(uE.rows());
+  for(int e = 0;e<EMAP.rows();e++)
+  {
+    N(EMAP(e))++;
+  }
+  B.resize(E.rows());
+  // Look of occurances of 2: one for original and another for boundary
+  for(int e = 0;e<E.rows();e++)
+  {
+    B(e) = (N(EMAP(e)) == 2);
+  }
+}
+
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DerivedB,
+  typename DerivedEMAP>
+void igl::is_boundary_edge(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Should be triangles
+  assert(F.cols() == 3);
+  // number of faces
+  const int m = F.rows();
+  // Collect all directed edges after E
+  MatrixXi allE(3*m,2);
+  for(int e = 0;e<3;e++)
+  {
+    for(int f = 0;f<m;f++)
+    {
+      for(int c = 0;c<2;c++)
+      {
+        // 12 20 01
+        allE(m*e+f,c) = F(f,(c+1+e)%3);
+      }
+    }
+  }
+  // sort directed edges into undirected edges
+  MatrixXi sorted_allE,_;
+  sort(allE,2,true,sorted_allE,_);
+  // Determine unique undirected edges E and map to directed edges EMAP
+  unique_rows(sorted_allE,E,_,EMAP);
+  // Counts of occurances
+  VectorXi N = VectorXi::Zero(E.rows());
+  for(int e = 0;e<EMAP.rows();e++)
+  {
+    N(EMAP(e))++;
+  }
+  B.resize(E.rows());
+  // Look of occurances of 1
+  for(int e = 0;e<E.rows();e++)
+  {
+    B(e) = (N(EMAP(e)) == 1);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit instanciation:
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::is_boundary_edge<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/is_boundary_edge.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IS_BOUNDARY_EDGE_H
+#define IS_BOUNDARY_EDGE_H
+#include <Eigen/Dense>
+
+namespace igl
+{
+  //  IS_BOUNDARY_EDGE Determine for each edge E if it is a "boundary edge" in F.
+  //  Boundary edges are undirected edges which occur only once.
+  // 
+  //  Inputs:
+  //    E  #E by 2 list of edges
+  //    F  #F by 3 list of triangles
+  //  Outputs:
+  //    B  #E list bools. true iff unoriented edge occurs exactly once in F
+  //      (non-manifold and non-existant edges will be false)
+  // 
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DerivedB>
+  void is_boundary_edge(
+    const Eigen::PlainObjectBase<DerivedE> & E,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedB> & B);
+  // Wrapper where Edges should also be computed from F
+  //   E  #E by 2 list of edges
+  //   EMAP  #F*3 list of indices mapping allE to E
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DerivedB,
+    typename DerivedEMAP>
+  void is_boundary_edge(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "is_boundary_edge.cpp"
+#endif
+
+#endif

include/igl/is_symmetric.cpp

@@ -49,7 +49,6 @@ IGL_INLINE bool igl::is_symmetric(
   const epsilonT epsilon)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   if(A.rows() != A.cols())
   {

include/igl/is_vertex_manifold.cpp

@@ -74,7 +74,7 @@ IGL_INLINE bool igl::is_vertex_manifold(
         }
       }
     }
-    return one_ring_size == seen.size();
+    return one_ring_size == (FIndex) seen.size();
   };
 
   // Unreferenced vertices are considered non-manifold
@@ -87,3 +87,7 @@ IGL_INLINE bool igl::is_vertex_manifold(
   }
   return all;
 }
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::is_vertex_manifold<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/jet.cpp

@@ -12,7 +12,6 @@
 //void igl::jet(const int m, Eigen::MatrixXd & J)
 //{
 //  using namespace Eigen;
-//  using namespace igl;
 //  // Applications/MATLAB_R2012b.app/toolbox/matlab/graph3d/jet.m
 //  const int n = ceil(m/4);
 //  // resize output

include/igl/lens_flare.cpp

@@ -65,7 +65,6 @@ IGL_INLINE void igl::lens_flare_create(
   const float * C,
   std::vector<igl::Flare> & flares)
 {
-  using namespace igl;
   using namespace std;
   flares.resize(12);
   /* Shines */

include/igl/line_segment_in_rectangle.cpp

@@ -15,7 +15,6 @@ IGL_INLINE bool igl::line_segment_in_rectangle(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // http://stackoverflow.com/a/100165/148668
   const auto SegmentIntersectRectangle = [](double a_rectangleMinX,
                                  double a_rectangleMinY,

include/igl/linprog.cpp

@@ -25,7 +25,6 @@ IGL_INLINE bool igl::linprog(
   // http://www.mathworks.com/matlabcentral/fileexchange/2166-introduction-to-linear-algebra/content/strang/linprog.m
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   bool success = true;
   // number of constraints
   const int m = _A.rows();
@@ -259,7 +258,6 @@ IGL_INLINE bool igl::linprog(
   Eigen::VectorXd & x)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   const int m = A.rows();
   const int n = A.cols();

include/igl/lscm.cpp

@@ -9,13 +9,12 @@
 
 #include <igl/vector_area_matrix.h>
 #include <igl/cotmatrix.h>
-//#include <igl/kronecker_product.h>
 #include <igl/repdiag.h>
 #include <igl/min_quad_with_fixed.h>
 #include <igl/matlab_format.h>
 #include <iostream>
 
-IGL_INLINE void igl::lscm(
+IGL_INLINE bool igl::lscm(
   const Eigen::MatrixXd& V,
   const Eigen::MatrixXi& F,
   const Eigen::VectorXi& b,
@@ -48,11 +47,16 @@ IGL_INLINE void igl::lscm(
   SparseMatrix<double> Q = -L_flat + 2.*A;
   const VectorXd B_flat = VectorXd::Zero(V.rows()*2);
   igl::min_quad_with_fixed_data<double> data;
-  igl::min_quad_with_fixed_precompute(Q,b_flat,SparseMatrix<double>(),true,data);
+  if(!igl::min_quad_with_fixed_precompute(Q,b_flat,SparseMatrix<double>(),true,data))
+  {
+    return false;
+  }
 
   MatrixXd W_flat;
   if(!min_quad_with_fixed_solve(data,B_flat,bc_flat,VectorXd(),W_flat))
-    assert(false);
+  {
+    return false;
+  }
 
 
   assert(W_flat.rows() == V.rows()*2);
@@ -61,7 +65,7 @@ IGL_INLINE void igl::lscm(
   {
     V_uv.col(i) = W_flat.block(V_uv.rows()*i,0,V_uv.rows(),1);
   }
-
+  return true;
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/lscm.h

@@ -14,9 +14,12 @@
 
 namespace igl
 {
-  // Compute a Least-squares conformal map parametrization following the algorithm
-  // presented in: Spectral Conformal Parameterization,
-  //               Patrick Mullen, Yiying Tong, Pierre Alliez and Mathieu Desbrun
+  // Compute a Least-squares conformal map parametrization following the
+  // algorithm presented in: Spectral Conformal Parameterization, Patrick
+  // Mullen, Yiying Tong, Pierre Alliez and Mathieu Desbrun. Input should be a
+  // manifold mesh (also no unreferenced vertices) and "boundary" `b` should
+  // contain at least two vertices per connected component.
+  //
   //
   // Inputs:
   //   V  #V by 3 list of mesh vertex positions
@@ -25,15 +28,9 @@ namespace igl
   //   bc #b by 3 list of boundary values
   // Outputs:
   //   UV #V by 2 list of 2D mesh vertex positions in UV space
+  // Returns true only on solver success.
   //
-  // X  Note: if b and bc are empty, lscm automatically removes the null space
-  // X        by fixing two farthest points on the boundary
-  // Boundary conditions are part of the api. It would be strange to secretly
-  // use other boundary conditions. This is also weird if there's more than one
-  // boundary loop.
-  // X  Note: (V,F) must be a genus-0 mesh, with a single boundary
-  // No longer the case. Should probably just be manifold.
-  IGL_INLINE void lscm(
+  IGL_INLINE bool lscm(
     const Eigen::MatrixXd& V,
     const Eigen::MatrixXi& F,
     const Eigen::VectorXi& b,

include/igl/matlab_format.cpp

@@ -14,7 +14,6 @@ IGL_INLINE const Eigen::WithFormat< DerivedM > igl::matlab_format(
   const Eigen::PlainObjectBase<DerivedM> & M,
   const std::string name)
 {
-  using namespace igl;
   using namespace std;
   string prefix = "";
   if(!name.empty())
@@ -42,7 +41,6 @@ igl::matlab_format(
   const std::string name)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   Matrix<typename Eigen::SparseMatrix<DerivedS>::Scalar,Dynamic,1> I,J,V;
   Matrix<DerivedS,Dynamic,Dynamic> SIJV;

include/igl/median.cpp

@@ -14,7 +14,6 @@
 IGL_INLINE bool igl::median(const Eigen::VectorXd & V, double & m)
 {
   using namespace std;
-  using namespace igl;
   if(V.size() == 0)
   {
     return false;

include/igl/min_quad_with_fixed.cpp

@@ -38,7 +38,6 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
 //#define MIN_QUAD_WITH_FIXED_CPP_DEBUG
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   const Eigen::SparseMatrix<T> A = 0.5*A2;
 #ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
   cout<<"    pre"<<endl;
@@ -383,7 +382,6 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 {
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   typedef Matrix<T,Dynamic,1> VectorXT;
   typedef Matrix<T,Dynamic,Dynamic> MatrixXT;
   // number of known rows

include/igl/on_boundary.cpp

@@ -20,7 +20,6 @@ IGL_INLINE void igl::on_boundary(
   std::vector<std::vector<bool> > & C)
 {
   using namespace std;
-  using namespace igl;
 
   // Get a list of all faces
   vector<vector<IntegerT> > F(T.size()*4,vector<IntegerT>(3));
@@ -77,7 +76,6 @@ IGL_INLINE void igl::on_boundary(
   assert(T.cols() == 0 || T.cols() == 4);
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // Cop out: use vector of vectors version
   vector<vector<typename Eigen::PlainObjectBase<DerivedT>::Scalar> > vT;
   matrix_to_list(T,vT);

include/igl/outer_hull.cpp

@@ -42,6 +42,17 @@ IGL_INLINE void igl::outer_hull(
   typedef Matrix<typename DerivedN::Scalar,1,3> RowVector3N;
   const Index m = F.rows();
 
+  const auto & duplicate_simplex = [&F](const int f, const int g)->bool
+  {
+    return 
+      (F(f,0) == F(g,0) && F(f,1) == F(g,1) && F(f,2) == F(g,2)) ||
+      (F(f,1) == F(g,0) && F(f,2) == F(g,1) && F(f,0) == F(g,2)) ||
+      (F(f,2) == F(g,0) && F(f,0) == F(g,1) && F(f,1) == F(g,2)) ||
+      (F(f,0) == F(g,2) && F(f,1) == F(g,1) && F(f,2) == F(g,0)) ||
+      (F(f,1) == F(g,2) && F(f,2) == F(g,1) && F(f,0) == F(g,0)) ||
+      (F(f,2) == F(g,2) && F(f,0) == F(g,1) && F(f,1) == F(g,0));
+  };
+
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"outer hull..."<<endl;
 #endif
@@ -107,11 +118,32 @@ IGL_INLINE void igl::outer_hull(
       // duplicate faces the same way, regardless of their orientations
       di_I(fei,1) = (cons[fei]?1.:-1.)*f;
     }
-    VectorXi IM;
 
+    // Despite the effort to get stable normals the atan2 up doesn't
+    // compute (exactly) -θ for -n if it computes θ for n. So just
+    // explicitly check if there's a duplicate face
+    // Shitty O(val^2) implementation
+    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
+    {
+      const auto & fe = uE2E[ui][fei];
+      const auto f = fe % m;
+      for(size_t gei = fei+1;gei<uE2E[ui].size();gei++)
+      {
+        const auto & ge = uE2E[ui][gei];
+        const auto g = ge % m;
+        if(duplicate_simplex(f,g))
+        {
+#ifdef IGL_OUTER_HULL_DEBUG
+          cout<<"Forcing duplicate: "<<(f+1)<<","<<(g+1)<<endl;
+#endif
+          di_I(gei,0) = di_I(fei,0);
+        }
+      }
+    }
+    VectorXi IM;
     //igl::sort(di[ui],true,di[ui],IM);
-    // Sort, but break ties using index to ensure that duplicates always show
-    // up in same order.
+    // Sort, but break ties using "signed index" to ensure that duplicates
+    // always show up in same order.
     MatrixXd s_di_I;
     igl::sortrows(di_I,true,s_di_I,IM);
     di[ui].resize(uE2E[ui].size());
@@ -247,16 +279,19 @@ IGL_INLINE void igl::outer_hull(
         const int nfei_new = (diIM(e) + 2*val + e_cons*step*(flip(f)?-1:1))%val;
         const int nf = uE2E[EMAP(e)][nfei_new] % m;
         // Don't consider faces with identical dihedral angles
-        if(di[EMAP(e)][diIM(e)] != di[EMAP(e)][nfei_new])
+        if((di[EMAP(e)][diIM(e)] != di[EMAP(e)][nfei_new]))
 //#warning "THIS IS HACK, FIX ME"
 //        if( abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new]) < 1e-16 )
         {
-//#ifdef IGL_OUTER_HULL_DEBUG
-//        cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
-//          di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei_new]<<"| = "<<
-//            abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new])
-//            <<endl;
-//#endif
+#ifdef IGL_OUTER_HULL_DEBUG
+        cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
+          di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei_new]<<"| = "<<
+            abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new])
+            <<endl;
+#endif
+        
+
+
           // Only use this face if not already seen
           if(!FH[nf])
           {
@@ -517,7 +552,7 @@ IGL_INLINE void igl::outer_hull(
   Eigen::PlainObjectBase<Derivedflip> & flip)
 {
   Eigen::Matrix<typename DerivedV::Scalar,DerivedF::RowsAtCompileTime,3> N;
-  per_face_normals(V,F,N);
+  per_face_normals_stable(V,F,N);
   return outer_hull(V,F,N,G,J,flip);
 }
 

include/igl/per_corner_normals.cpp

@@ -18,7 +18,6 @@ IGL_INLINE void igl::per_corner_normals(
   const double corner_threshold,
   Eigen::PlainObjectBase<DerivedV> & CN)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   Eigen::PlainObjectBase<DerivedV> FN;
@@ -36,7 +35,6 @@ IGL_INLINE void igl::per_corner_normals(
   const double corner_threshold,
   Eigen::PlainObjectBase<DerivedCN> & CN)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
   vector<vector<int> > VF,VFi;
@@ -53,7 +51,6 @@ IGL_INLINE void igl::per_corner_normals(
   const double corner_threshold,
   Eigen::PlainObjectBase<DerivedV> & CN)
 {
-  using namespace igl;
   using namespace Eigen;
   using namespace std;
 

include/igl/per_face_normals.cpp

@@ -47,6 +47,65 @@ IGL_INLINE void igl::per_face_normals(
   return per_face_normals(V,F,Z,N);
 }
 
+template <typename DerivedV, typename DerivedF, typename DerivedN>
+IGL_INLINE void igl::per_face_normals_stable(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedN> & N)
+{
+  using namespace Eigen;
+  typedef Matrix<typename DerivedN::Scalar,DerivedN::RowsAtCompileTime,3> MatrixN3;
+  typedef Matrix<typename DerivedV::Scalar,DerivedF::RowsAtCompileTime,3> MatrixV3;
+  typedef Matrix<typename DerivedV::Scalar,3,3> MatrixV33;
+  typedef Matrix<typename DerivedV::Scalar,1,3> RowVectorV3;
+  typedef typename DerivedV::Scalar Scalar;
+
+  const size_t m = F.rows();
+
+  N.resize(F.rows(),3);
+  // Grad all points
+  for(size_t f = 0;f<m;f++)
+  {
+    const RowVectorV3 p0 = V.row(F(f,0));
+    const RowVectorV3 p1 = V.row(F(f,1));
+    const RowVectorV3 p2 = V.row(F(f,2));
+    const RowVectorV3 n0 = (p1 - p0).cross(p2 - p0);
+    const RowVectorV3 n1 = (p2 - p1).cross(p0 - p1);
+    const RowVectorV3 n2 = (p0 - p2).cross(p1 - p2);
+
+    // careful sum
+    for(int d = 0;d<3;d++)
+    {
+      // This is a little _silly_ in terms of complexity, but its recursive
+      // implementation is clean looking...
+      const std::function<Scalar(Scalar,Scalar,Scalar)> sum3 = 
+        [&sum3](Scalar a, Scalar b, Scalar c)->Scalar
+      {
+        if(fabs(c)>fabs(a))
+        {
+          return sum3(c,b,a);
+        }
+        // c < a
+        if(fabs(c)>fabs(b))
+        {
+          return sum3(a,c,b);
+        }
+        // c < a, c < b
+        if(fabs(b)>fabs(a))
+        {
+          return sum3(b,a,c);
+        }
+        return (a+b)+c;
+      };
+
+      N(f,d) = sum3(n0(d),n1(d),n2(d));
+    }
+    // sum better not be sure, or else NaN
+    N.row(f) /= N.row(f).norm();
+  }
+
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
@@ -54,4 +113,6 @@ template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Ei
 template void igl::per_face_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_face_normals<Eigen::Matrix<float, -1, -1, 1, -1, -1>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 1, -1, -1> >&);
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::per_face_normals_stable<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::per_face_normals_stable<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/per_face_normals.h

@@ -35,6 +35,13 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedN> & N);
+  // Special version where order of face indices is guaranteed not to effect
+  // output.
+  template <typename DerivedV, typename DerivedF, typename DerivedN>
+  IGL_INLINE void per_face_normals_stable(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedN> & N);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/quat_to_axis_angle.cpp

@@ -18,7 +18,6 @@ IGL_INLINE void igl::quat_to_axis_angle(
   Q_type *axis, 
   Q_type & angle)
 {
-  using namespace igl;
     if( fabs(q[3])>(1.0 + igl::EPS<Q_type>()) )
     {
         //axis[0] = axis[1] = axis[2] = 0; // no, keep the previous value

include/igl/random_dir.cpp

@@ -12,7 +12,6 @@
 IGL_INLINE Eigen::Vector3d igl::random_dir()
 {
   using namespace Eigen;
-  using namespace igl;
   double z =  (double)rand() / (double)RAND_MAX*2.0 - 1.0;
   double t =  (double)rand() / (double)RAND_MAX*2.0*PI;
   // http://www.altdevblogaday.com/2012/05/03/generating-uniformly-distributed-points-on-sphere/
@@ -25,7 +24,6 @@ IGL_INLINE Eigen::Vector3d igl::random_dir()
 IGL_INLINE Eigen::MatrixXd igl::random_dir_stratified(const int n)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   const double m = floor(sqrt(double(n)));
   MatrixXd N(n,3);

include/igl/random_quaternion.cpp

@@ -0,0 +1,23 @@
+#include "random_quaternion.h"
+
+template <typename Scalar>
+IGL_INLINE Eigen::Quaternion<Scalar> igl::random_quaternion()
+{
+  // http://mathproofs.blogspot.com/2005/05/uniformly-distributed-random-unit.html
+  const auto & unit_rand = []()->Scalar
+  {
+    return ((Scalar)rand() / (Scalar)RAND_MAX);
+  };
+  const Scalar t0 = 2.*M_PI*unit_rand();
+  const Scalar t1 = acos(1.-2.*unit_rand());
+  const Scalar t2 = 0.5*(M_PI*unit_rand() + acos(unit_rand()));
+  return Eigen::Quaternion<Scalar>(
+    1.*sin(t0)*sin(t1)*sin(t2),
+    1.*cos(t0)*sin(t1)*sin(t2),
+    1.*cos(t1)*sin(t2),
+    1.*cos(t2));
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/random_quaternion.h

@@ -0,0 +1,14 @@
+#ifndef IGL_RANDOM_QUATERNION_H
+#define IGL_RANDOM_QUATERNION_H
+#include "igl_inline.h"
+#include <Eigen/Geometry>
+namespace igl
+{
+  // Return a random quaternion via uniform sampling of the 4-sphere
+  template <typename Scalar>
+  IGL_INLINE Eigen::Quaternion<Scalar> random_quaternion();
+}
+#ifndef IGL_STATIC_LIBRARY
+#include "random_quaternion.cpp"
+#endif
+#endif

include/igl/readCSV.cpp

@@ -20,7 +20,6 @@ IGL_INLINE bool igl::readCSV(
   Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>& M)
 {
   using namespace std;
-  using namespace igl;
 
   std::vector<std::vector<Scalar> > Mt;
   

include/igl/readMESH.cpp

@@ -18,7 +18,6 @@ IGL_INLINE bool igl::readMESH(
   std::vector<std::vector<Index > > & F)
 {
   using namespace std;
-  using namespace igl;
   FILE * mesh_file = fopen(mesh_file_name.c_str(),"r");
   if(NULL==mesh_file)
   {
@@ -233,7 +232,6 @@ IGL_INLINE bool igl::readMESH(
   Eigen::PlainObjectBase<DerivedF>& F)
 {
   using namespace std;
-  using namespace igl;
   FILE * mesh_file = fopen(mesh_file_name.c_str(),"r");
   if(NULL==mesh_file)
   {

include/igl/readNODE.cpp

@@ -10,7 +10,7 @@
 #include <stdio.h>
 
 template <typename Scalar, typename Index>
-IGL_INLINE bool readNODE(
+IGL_INLINE bool igl::readNODE(
   const std::string node_file_name,
   std::vector<std::vector<Scalar > > & V,
   std::vector<std::vector<Index > > & I)
@@ -36,7 +36,6 @@ IGL_INLINE bool igl::readNODE(
   Eigen::PlainObjectBase<DerivedI>& I)
 {
   using namespace std;
-  using namespace igl;
   FILE * node_file = fopen(node_file_name.c_str(),"r");
   if(NULL==node_file)
   {

include/igl/readTGF.cpp

@@ -148,7 +148,6 @@ IGL_INLINE bool igl::readTGF(
   Eigen::MatrixXi & CE,
   Eigen::MatrixXi & PE)
 {
-  using namespace igl;
   std::vector<std::vector<double> > vC;
   std::vector<std::vector<int> > vE;
   std::vector<int> vP;

include/igl/read_triangle_mesh.cpp

@@ -30,7 +30,6 @@ IGL_INLINE bool igl::read_triangle_mesh(
   std::vector<std::vector<Index> > & F)
 {
   using namespace std;
-  using namespace igl;
   // dirname, basename, extension and filename
   string d,b,e,f;
   pathinfo(str,d,b,e,f);

include/igl/remove_duplicate_vertices.cpp

@@ -24,7 +24,6 @@ IGL_INLINE void igl::remove_duplicate_vertices(
   Eigen::PlainObjectBase<DerivedSVI>& SVI,
   Eigen::PlainObjectBase<DerivedSVJ>& SVJ)
 {
-  using namespace igl;
   if(epsilon > 0)
   {
     Eigen::PlainObjectBase<DerivedV> rV,rSV;

include/igl/reorder.cpp

@@ -37,4 +37,5 @@ template void igl::reorder<int>(std::vector<int, std::allocator<int> > const&, s
   template void igl::reorder<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&);
   template void igl::reorder<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > >(std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> >, std::allocator<igl::SortableRow<Eigen::Matrix<double, -1, 1, 0, -1, 1> > > >&);
 #  endif
+template void igl::reorder<long>(std::vector<long, std::allocator<long> > const&, std::vector<unsigned long, std::allocator<unsigned long> > const&, std::vector<long, std::allocator<long> >&);
 #endif

include/igl/sort.cpp

@@ -147,13 +147,13 @@ IGL_INLINE void igl::sort2(
 {
   using namespace Eigen;
   using namespace std;
-  Y = X;
+  typedef typename Eigen::PlainObjectBase<DerivedY>::Scalar YScalar;
+  Y = X.template cast<YScalar>();
   // get number of columns (or rows)
   int num_outer = (dim == 1 ? X.cols() : X.rows() );
   // get number of rows (or columns)
   int num_inner = (dim == 1 ? X.rows() : X.cols() );
   assert(num_inner == 2);(void)num_inner;
-  typedef typename Eigen::PlainObjectBase<DerivedX>::Scalar Scalar;
   typedef typename Eigen::PlainObjectBase<DerivedIX>::Scalar Index;
   IX.resize(X.rows(),X.cols());
   if(dim==1)
@@ -168,8 +168,8 @@ IGL_INLINE void igl::sort2(
   // loop over columns (or rows)
   for(int i = 0;i<num_outer;i++)
   {
-    Scalar & a = (dim==1 ? Y(0,i) : Y(i,0));
-    Scalar & b = (dim==1 ? Y(1,i) : Y(i,1));
+    YScalar & a = (dim==1 ? Y(0,i) : Y(i,0));
+    YScalar & b = (dim==1 ? Y(1,i) : Y(i,1));
     Index & ai = (dim==1 ? IX(0,i) : IX(i,0));
     Index & bi = (dim==1 ? IX(1,i) : IX(i,1));
     if((ascending && a>b) || (!ascending && a<b))
@@ -226,4 +226,6 @@ template void igl::sort_new<Eigen::Matrix<int, 1, 6, 1, 1, 6>, Eigen::Matrix<int
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::sort<Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, 4, 0, -1, 4> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sort<long>(std::vector<long, std::allocator<long> > const&, bool, std::vector<long, std::allocator<long> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
+template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/sort_triangles.cpp

@@ -32,7 +32,6 @@ IGL_INLINE void igl::sort_triangles(
   Eigen::PlainObjectBase<DerivedI> & I)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
 
@@ -101,7 +100,6 @@ void igl::sort_triangles(
   Eigen::PlainObjectBase<DerivedI> & I)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   // Put model, projection, and viewport matrices into double arrays
   Matrix4d MV;
@@ -136,7 +134,6 @@ void igl::sort_triangles_slow(
   Eigen::PlainObjectBase<DerivedI> & I)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   // Barycenter, centroid
   Eigen::Matrix<typename DerivedV::Scalar, DerivedF::RowsAtCompileTime,1> D,sD;
@@ -463,7 +460,6 @@ class Triangle
 //    "must be drawn before that one'");
 //  using namespace std;
 //  using namespace Eigen;
-//  using namespace igl;
 //  typedef Matrix<typename DerivedV::Scalar,3,1> Vec3;
 //  assert(V.cols() == 4);
 //  Matrix<typename DerivedV::Scalar, DerivedV::RowsAtCompileTime,3> VMVP =
@@ -542,7 +538,6 @@ class Triangle
 //  Eigen::PlainObjectBase<DerivedI> & I)
 //{
 //  using namespace Eigen;
-//  using namespace igl;
 //  using namespace std;
 //  // Put model, projection, and viewport matrices into double arrays
 //  Matrix4d MV;

include/igl/sortrows.cpp

@@ -61,7 +61,6 @@ IGL_INLINE void igl::sortrows(
   // improvement.
   using namespace std;
   using namespace Eigen;
-  using namespace igl;
   // Resize output
   Y.resize(X.rows(),X.cols());
   IX.resize(X.rows(),1);

include/igl/texture_from_tga.cpp

@@ -15,7 +15,6 @@
 IGL_INLINE bool igl::texture_from_tga(const std::string tga_file, GLuint & id)
 {
   using namespace std;
-  using namespace igl;
 
   // read pixels to tga file
   FILE * imgFile;

include/igl/trackball.cpp

@@ -22,7 +22,8 @@
 template <typename Q_type>
 static IGL_INLINE Q_type _QuatD(double w, double h)
 {
-  return (Q_type)(abs(w) < abs(h) ? abs(w) : abs(h)) - 4;
+  using namespace std;
+  return (Q_type)(std::abs(w) < std::abs(h) ? std::abs(w) : std::abs(h)) - 4;
 }
 template <typename Q_type>
 static IGL_INLINE Q_type _QuatIX(double x, double w, double h)

include/igl/triangle/triangulate.cpp

@@ -111,10 +111,6 @@ IGL_INLINE void igl::triangulate(
   free(in.segmentmarkerlist);
   free(in.holelist);
 
-  // Cleanup out
-  free(out.pointlist);
-  free(out.trianglelist);
-  free(out.segmentlist);
 
   // Return the mesh
   V2.resize(out.numberofpoints,2);
@@ -127,6 +123,11 @@ IGL_INLINE void igl::triangulate(
     for (unsigned j=0;j<3;++j)
       F2(i,j) = out.trianglelist[i*3+j];
 
+  // Cleanup out
+  free(out.pointlist);
+  free(out.trianglelist);
+  free(out.segmentlist);
+
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/triangle_triangle_adjacency.cpp

@@ -11,6 +11,7 @@
 #include "unique_simplices.h"
 #include "unique_edge_map.h"
 #include <algorithm>
+#include <iostream>
 
 template <typename Scalar, typename Index>
 IGL_INLINE void igl::triangle_triangle_adjacency_preprocess(const Eigen::PlainObjectBase<Scalar>& /*V*/,
@@ -211,4 +212,6 @@ template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&);
 #endif

include/igl/uniformly_sample_two_manifold.cpp

@@ -29,7 +29,6 @@ IGL_INLINE void igl::uniformly_sample_two_manifold(
   Eigen::MatrixXd & WS)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
   // Euclidean distance between two points on a mesh given as barycentric
@@ -351,7 +350,6 @@ IGL_INLINE void igl::uniformly_sample_two_manifold_at_vertices(
   Eigen::VectorXi & S)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
 
   // Copy weights and faces

include/igl/unique.cpp

@@ -208,7 +208,6 @@ IGL_INLINE void igl::unique_rows(
   Eigen::PlainObjectBase<DerivedIC>& IC)
 {
   using namespace std;
-  using namespace igl;
   using namespace Eigen;
   VectorXi IM;
   Eigen::PlainObjectBase<DerivedA> sortA;
@@ -260,4 +259,5 @@ template void igl::unique<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<i
 template void igl::unique<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::unique_rows<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::unique<int>(std::vector<int, std::allocator<int> > const&, std::vector<int, std::allocator<int> >&);
+template void igl::unique<long>(std::vector<long, std::allocator<long> > const&, std::vector<long, std::allocator<long> >&, std::vector<unsigned long, std::allocator<unsigned long> >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
 #endif

include/igl/unique_edge_map.cpp

@@ -42,4 +42,5 @@ template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::M
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+template void igl::unique_edge_map<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 #endif 

include/igl/unique_simplices.cpp

@@ -22,7 +22,6 @@ IGL_INLINE void igl::unique_simplices(
   Eigen::PlainObjectBase<DerivedIC>& IC)
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   // Sort each face
   MatrixXi sortF, unusedI;
@@ -62,4 +61,5 @@ template void igl::unique_simplices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen
 template void igl::unique_simplices<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::unique_simplices<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::unique_simplices<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::unique_simplices<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/upsample.cpp

@@ -25,7 +25,6 @@ IGL_INLINE void igl::upsample(
   // Use "in place" wrapper instead
   assert(&V != &NV);
   assert(&F != &NF);
-  using namespace igl;
   using namespace std;
   using namespace Eigen;
 

include/igl/vertex_triangle_adjacency.cpp

@@ -48,4 +48,6 @@ IGL_INLINE void igl::vertex_triangle_adjacency(
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, unsigned int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > >&);
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 #endif

include/igl/writeMESH.cpp

@@ -52,7 +52,6 @@ IGL_INLINE bool igl::writeMESH(
   const Eigen::PlainObjectBase<DerivedF> & F)
 {
   using namespace std;
-  using namespace igl;
   using namespace Eigen;
 
   //// This is (surprisingly) slower than the C-ish code below

include/igl/writeTGF.cpp

@@ -63,7 +63,6 @@ IGL_INLINE bool igl::writeTGF(
   const Eigen::MatrixXd & C,
   const Eigen::MatrixXi & E)
 {
-  using namespace igl;
   using namespace std;
   vector<vector<double> > vC;
   vector<vector<int> > vE;

tutorial/504_NRosyDesign/CMakeLists.txt

@@ -5,7 +5,8 @@ include("../CMakeLists.shared")
 
 find_package(LIBCOMISO REQUIRED)
 
-include_directories( ${LIBCOMISO_INCLUDE_DIR} )
+include_directories( ${LIBCOMISO_INCLUDE_DIR}
+                     ${LIBCOMISO_INCLUDE_DIRS} )
 
 set(SOURCES
 ${PROJECT_SOURCE_DIR}/main.cpp

tutorial/505_MIQ/CMakeLists.txt

@@ -5,7 +5,8 @@ include("../CMakeLists.shared")
 
 find_package(LIBCOMISO REQUIRED)
 
-include_directories( ${LIBCOMISO_INCLUDE_DIR} )
+include_directories( ${LIBCOMISO_INCLUDE_DIR}
+                     ${LIBCOMISO_INCLUDE_DIRS} )
 
 set(SOURCES
 ${PROJECT_SOURCE_DIR}/main.cpp

tutorial/505_MIQ/main.cpp

@@ -16,6 +16,7 @@ Eigen::MatrixXd B;
 
 // Scale for visualizing the fields
 double global_scale;
+bool extend_arrows = false;
 
 // Cross field
 Eigen::MatrixXd X1,X2;
@@ -70,6 +71,11 @@ void line_texture(Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> &te
 
 bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
 {
+  if (key == 'E')
+  {
+    extend_arrows = !extend_arrows;
+  }
+  
   if (key <'1' || key >'8')
     return false;
 
@@ -81,24 +87,24 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
   {
     // Cross field
     viewer.data.set_mesh(V, F);
-    viewer.data.add_edges(B, B + global_scale*X1 ,Eigen::RowVector3d(1,0,0));
-    viewer.data.add_edges(B, B + global_scale*X2 ,Eigen::RowVector3d(0,0,1));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*X1 : B, B + global_scale*X1 ,Eigen::RowVector3d(1,0,0));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*X2 : B, B + global_scale*X2 ,Eigen::RowVector3d(0,0,1));
   }
 
   if (key == '2')
   {
     // Bisector field
     viewer.data.set_mesh(V, F);
-    viewer.data.add_edges(B, B + global_scale*BIS1 ,Eigen::RowVector3d(1,0,0));
-    viewer.data.add_edges(B, B + global_scale*BIS2 ,Eigen::RowVector3d(0,0,1));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*BIS1 : B, B + global_scale*BIS1 ,Eigen::RowVector3d(1,0,0));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*BIS2 : B, B + global_scale*BIS2 ,Eigen::RowVector3d(0,0,1));
   }
 
   if (key == '3')
   {
     // Bisector field combed
     viewer.data.set_mesh(V, F);
-    viewer.data.add_edges(B, B + global_scale*BIS1_combed ,Eigen::RowVector3d(1,0,0));
-    viewer.data.add_edges(B, B + global_scale*BIS2_combed ,Eigen::RowVector3d(0,0,1));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*BIS1_combed : B, B + global_scale*BIS1_combed ,Eigen::RowVector3d(1,0,0));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*BIS2_combed : B, B + global_scale*BIS2_combed ,Eigen::RowVector3d(0,0,1));
   }
 
   if (key == '4')
@@ -142,8 +148,8 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     // Singularities and cuts, original field
     // Singularities and cuts
     viewer.data.set_mesh(V, F);
-    viewer.data.add_edges(B, B + global_scale*X1_combed ,Eigen::RowVector3d(1,0,0));
-    viewer.data.add_edges(B, B + global_scale*X2_combed ,Eigen::RowVector3d(0,0,1));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*X1_combed : B, B + global_scale*X1_combed ,Eigen::RowVector3d(1,0,0));
+    viewer.data.add_edges(extend_arrows ? B - global_scale*X2_combed : B, B + global_scale*X2_combed ,Eigen::RowVector3d(0,0,1));
 
     // Plot cuts
     int l_count = Seams.sum();
@@ -198,7 +204,7 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     viewer.data.set_uv(UV_seams,FUV_seams);
     viewer.core.show_texture = true;
   }
-
+  
   viewer.data.set_colors(Eigen::RowVector3d(1,1,1));
 
   // Replace the standard texture with an integer shift invariant texture

tutorial/506_FrameField/CMakeLists.txt

@@ -5,7 +5,8 @@ include("../CMakeLists.shared")
 
 find_package(LIBCOMISO REQUIRED)
 
-include_directories( ${LIBCOMISO_INCLUDE_DIR} )
+include_directories( ${LIBCOMISO_INCLUDE_DIR}
+                     ${LIBCOMISO_INCLUDE_DIRS} )
 
 set(SOURCES
 ${PROJECT_SOURCE_DIR}/main.cpp

tutorial/703_Decimation/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(703_Decimation)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES})

tutorial/703_Decimation/main.cpp

@@ -0,0 +1,133 @@
+#include <igl/collapse_edge.h>
+#include <igl/edge_flaps.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/circulation.h>
+#include <igl/Viewer/viewer.h>
+#include <Eigen/Core>
+#include <iostream>
+#include <set>
+
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  cout<<"Usage: ./progressive_hulls [filename.(off|obj|ply)]"<<endl;
+  cout<<"  [space]  toggle animation."<<endl;
+  cout<<"  'r'  reset."<<endl;
+  // Load a closed manifold mesh
+  string filename("../shared/fertility.off");
+  if(argc>=2)
+  {
+    filename = argv[1];
+  }
+  MatrixXd V,OV;
+  MatrixXi F,OF;
+  read_triangle_mesh(filename,OV,OF);
+
+  igl::Viewer viewer;
+
+  // Prepare array-based edge data structures and priority queue
+  VectorXi EMAP;
+  MatrixXi E,EF,EI;
+  typedef std::set<std::pair<double,int> > PriorityQueue;
+  PriorityQueue Q;
+  std::vector<PriorityQueue::iterator > Qit;
+  // If an edge were collapsed, we'd collapse it to these points:
+  MatrixXd C;
+  int num_collapsed;
+
+  // Function for computing cost of collapsing edge (lenght) and placement
+  // (midpoint)
+  const auto & shortest_edge_and_midpoint = [](
+    const int e,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & /*F*/,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & /*EMAP*/,
+    const Eigen::MatrixXi & /*EF*/,
+    const Eigen::MatrixXi & /*EI*/,
+    double & cost,
+    RowVectorXd & p)
+  {
+    cost = (V.row(E(e,0))-V.row(E(e,1))).norm();
+    p = 0.5*(V.row(E(e,0))+V.row(E(e,1)));
+  };
+
+
+  // Function to reset original mesh and data structures
+  const auto & reset = [&]()
+  {
+    F = OF;
+    V = OV;
+    edge_flaps(F,E,EMAP,EF,EI);
+    Qit.resize(E.rows());
+
+    C.resize(E.rows(),V.cols());
+    VectorXd costs(E.rows());
+    for(int e = 0;e<E.rows();e++)
+    {
+      double cost = e;
+      RowVectorXd p(1,3);
+      shortest_edge_and_midpoint(e,V,F,E,EMAP,EF,EI,cost,p);
+      C.row(e) = p;
+      Qit[e] = Q.insert(std::pair<double,int>(cost,e)).first;
+    }
+    num_collapsed = 0;
+    viewer.data.clear();
+    viewer.data.set_mesh(V,F);
+    viewer.data.set_face_based(true);
+  };
+
+  const auto &pre_draw = [&](igl::Viewer & viewer)->bool
+  {
+    // If animating then collapse 10% of edges
+    if(viewer.core.is_animating && !Q.empty())
+    {
+      bool something_collapsed = false;
+      // collapse edge 
+      const int max_iter = std::ceil(0.01*Q.size());
+      for(int j = 0;j<max_iter;j++)
+      {
+        if(!collapse_edge(shortest_edge_and_midpoint,V,F,E,EMAP,EF,EI,Q,Qit,C))
+        {
+          break;
+        }
+        something_collapsed = true;
+        num_collapsed++;
+      }
+
+      if(something_collapsed)
+      {
+        viewer.data.clear();
+        viewer.data.set_mesh(V,F);
+        viewer.data.set_face_based(true);
+      }
+    }
+    return false;
+  };
+
+  const auto &key_down = 
+    [&](igl::Viewer &viewer,unsigned char key,int mod)->bool
+  {
+    switch(key)
+    {
+      case ' ':
+        viewer.core.is_animating ^= 1;
+        break;
+      case 'R':
+      case 'r':
+        reset();
+        break;
+      default:
+        return false;
+    }
+    return true;
+  };
+
+  reset();
+  viewer.core.is_animating = true;
+  viewer.callback_key_down = key_down;
+  viewer.callback_pre_draw = pre_draw;
+  return viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -79,3 +79,4 @@ add_subdirectory("609_Boolean")
 endif()
 add_subdirectory("701_Statistics")
 add_subdirectory("702_WindingNumber")
+add_subdirectory("703_Decimation")

tutorial/cmake/FindEIGEN.cmake

@@ -17,10 +17,10 @@
 
 if(NOT Eigen_FIND_VERSION)
   if(NOT Eigen_FIND_VERSION_MAJOR)
-    set(Eigen_FIND_VERSION_MAJOR 2)
+    set(Eigen_FIND_VERSION_MAJOR 3)
   endif(NOT Eigen_FIND_VERSION_MAJOR)
   if(NOT Eigen_FIND_VERSION_MINOR)
-    set(Eigen_FIND_VERSION_MINOR 91)
+    set(Eigen_FIND_VERSION_MINOR 2)
   endif(NOT Eigen_FIND_VERSION_MINOR)
   if(NOT Eigen_FIND_VERSION_PATCH)
     set(Eigen_FIND_VERSION_PATCH 0)

tutorial/cmake/FindLIBCOMISO.cmake

@@ -38,12 +38,32 @@ FIND_LIBRARY(LIBCOMISO_LIBRARY NAMES CoMISo
     ${PROJECT_SOURCE_DIR}/../../../CoMISo/
     ${PROJECT_SOURCE_DIR}/../../../CoMISo/build/Build/lib/CoMISo/
     /Users/olkido/Documents/igl/MIQ/src/CoMISo/Build
+    /usr/local/lib
+    /usr/local/lib/CoMISo
+    /usr/lib
+    /usr/lib/CoMISo
 )
 #message(STATUS "${LIBCOMISO_LIBRARY}")
 
 if(LIBCOMISO_INCLUDE_DIR AND LIBCOMISO_LIBRARY)
 
-   set(LIBCOMISO_INCLUDE_DIR ${LIBCOMISO_INCLUDE_DIR} ${LIBCOMISO_INCLUDE_DIR}/CoMISo/gmm/include)
+   #message("${LIBCOMISO_INCLUDE_DIR}")
+
+   set(LIBCOMISO_INCLUDE_DIRS
+      ${LIBCOMISO_INCLUDE_DIR}
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/Solver
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/EigenSolver
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/NSolver
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/Config
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/Utils
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/QtWidgets
+      ${LIBCOMISO_INCLUDE_DIR}/CoMISo/gmm/include
+      )
+
+   #message("${LIBCOMISO_INCLUDE_DIRS}")
+
+   set(LIBCOMISO_INCLUDE_DIR ${LIBCOMISO_INCLUDE_DIR})
 
    add_definitions(-DINCLUDE_TEMPLATES)
    message(STATUS "Found LIBCOMISO: ${LIBCOMISO_INCLUDE_DIR} ${LIBCOMISO_LIBRARY}")

tutorial/images/edge-collapse.pdf.REMOVED.git-id

@@ -0,0 +1 @@
+8bf78140fdc48fa6597d2579892b6c88bc7701e0

tutorial/images/fertility-edge-collapse.gif.REMOVED.git-id

@@ -0,0 +1 @@
+9c5fd841eac14ffcb1da4d52c7f2d6075f36c9cd