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

Former-commit-id: b834fbdc9b0ca6b4671bbac50be19b823e139918

.gitignore

@@ -95,3 +95,4 @@ python/build3
 python/build4
 python/scripts/generated
 python/builddebug
+tutorial/.idea

.gitmodules

@@ -6,9 +6,6 @@ url=https://github.com/libigl/nanogui.git
 [submodule "external/embree"]
 	path = external/embree
 	url = https://github.com/embree/embree.git
-[submodule "external/pybind11"]
-	path = external/pybind11
-	url = https://github.com/wjakob/pybind11.git
 [submodule "external/cgal"]
 	path = external/cgal
 	url = https://github.com/CGAL/cgal.git

README.md

@@ -194,8 +194,10 @@ Libigl is used by many research groups around the world. In 2015, it won the
 Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:
 
+ - [Activision](http://www.activision.com)
  - [Adobe Research](http://www.adobe.com/technology/)  
  - [Electronic Arts, Inc](http://www.ea.com)
+ - [Google Research](https://research.google.com)
  - [Mesh](http://meshconsultants.ca/), consultants, Canada
  - [Pixar Research](http://graphics.pixar.com/research/)
  - [Spine by Esoteric Software](http://esotericsoftware.com/) is an animation tool dedicated to 2D characters.
@@ -219,6 +221,7 @@ few labs/companies/institutions using libigl:
  - [University of Cambridge](http://www.cam.ac.uk/), England
  - [University of Pennsylvania](http://cg.cis.upenn.edu/), USA
  - [University of Texas at Austin](http://www.cs.utexas.edu/users/evouga/), USA
+ - [University of Toronto](http://dgp.toronto.edu), Canada
  - [University of Victoria](https://www.csc.uvic.ca/Research/graphics/), Canada
  - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
  - [Université Toulouse III Paul Sabatier](http://www.univ-tlse3.fr/), France

include/igl/AABB.cpp

@@ -10,7 +10,6 @@
 #include "barycenter.h"
 #include "barycentric_coordinates.h"
 #include "colon.h"
-#include "colon.h"
 #include "doublearea.h"
 #include "matlab_format.h"
 #include "point_simplex_squared_distance.h"

include/igl/cat.cpp

@@ -237,14 +237,14 @@ IGL_INLINE void igl::cat(const std::vector<std::vector< Mat > > & A, Mat & C)
   using namespace std;
   // Start with empty matrix
   C.resize(0,0);
-  for(typename vector<vector< Mat > >::const_iterator rit = A.begin(); rit != A.end(); rit++)
+  for(const auto & row_vec : A)
   {
     // Concatenate each row horizontally
     // Start with empty matrix
     Mat row(0,0);
-    for(typename vector<vector< Mat > >::iterator cit = A.begin(); rit != A.end(); rit++)
+    for(const auto & element : row_vec)
     {
-      row = cat(2,row,*cit);
+      row = cat(2,row,element);
     }
     // Concatenate rows vertically
     C = cat(1,C,row);

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -109,7 +109,7 @@ namespace igl
           // Make a short name for the edge map
           typedef std::map<EMK,EMV> EdgeMap;
           // Maps edges of offending faces to all incident offending faces
-          EdgeMap edge2faces;
+          //EdgeMap edge2faces;
           std::vector<std::pair<const Box, const Box> > candidate_box_pairs;
 
         public:
@@ -330,7 +330,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   T(),
   lIF(),
   offending(),
-  edge2faces(),
+  //edge2faces(),
   params(params)
 {
   using namespace std;
@@ -427,7 +427,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   }
 
   remesh_intersections(
-    V,F,T,offending,edge2faces,params.stitch_all,VV,FF,J,IM);
+    V,F,T,offending,params.stitch_all,VV,FF,J,IM);
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("remesh_intersection");
@@ -477,7 +477,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
       // append face to edge's list
       Index i = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+1)%3) : F(f,(e+2)%3);
       Index j = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+2)%3) : F(f,(e+1)%3);
-      edge2faces[EMK(i,j)].push_back(f);
+      //edge2faces[EMK(i,j)].push_back(f);
     }
   }
 }

include/igl/copyleft/cgal/assign_scalar.cpp

@@ -29,9 +29,38 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   } while (d < interval.second);
 }
 
+IGL_INLINE void igl::copyleft::cgal::assign_scalar(
+  const typename CGAL::Epeck::FT & _cgal,
+  float& d)
+{
+  // FORCE evaluation of the exact type otherwise interval might be huge.
+  const typename CGAL::Epeck::FT cgal = _cgal.exact();
+  const auto interval = CGAL::to_interval(cgal);
+  d = interval.first;
+  do {
+      const float next = nextafter(d, float(interval.second));
+      if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
+      d = next;
+  } while (d < float(interval.second));
+}
+
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const double & c,
   double & d)
 {
   d = c;
 }
+
+IGL_INLINE void igl::copyleft::cgal::assign_scalar(
+  const float& c,
+  float& d)
+{
+  d = c;
+}
+
+IGL_INLINE void igl::copyleft::cgal::assign_scalar(
+  const float& c,
+  double& d)
+{
+  d = c;
+}

include/igl/copyleft/cgal/assign_scalar.h

@@ -25,9 +25,18 @@ namespace igl
       IGL_INLINE void assign_scalar(
         const typename CGAL::Epeck::FT & cgal,
         double & d);
+      IGL_INLINE void assign_scalar(
+        const typename CGAL::Epeck::FT & cgal,
+        float& d);
       IGL_INLINE void assign_scalar(
         const double & c,
         double & d);
+      IGL_INLINE void assign_scalar(
+        const float& c,
+        float & d);
+      IGL_INLINE void assign_scalar(
+        const float& c,
+        double& d);
     }
   }
 }

include/igl/copyleft/cgal/intersect_other.cpp

@@ -9,6 +9,8 @@
 #include "CGAL_includes.hpp"
 #include "mesh_to_cgal_triangle_list.h"
 #include "remesh_intersections.h"
+#include "../../slice_mask.h"
+#include "../../remove_unreferenced.h"
 
 #ifndef IGL_FIRST_HIT_EXCEPTION
 #define IGL_FIRST_HIT_EXCEPTION 10
@@ -30,10 +32,10 @@ namespace igl
         std::map<
           typename DerivedF::Index,
           std::vector<std::pair<typename DerivedF::Index, CGAL::Object> > > &
-          offending,
-        std::map<
-          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-          std::vector<typename DerivedF::Index> > & edge2faces)
+          offending)
+        //std::map<
+        //  std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+        //  std::vector<typename DerivedF::Index> > & edge2faces)
       {
         typedef typename DerivedF::Index Index;
         typedef std::pair<Index,Index> EMK;
@@ -46,7 +48,7 @@ namespace igl
             // append face to edge's list
             Index i = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+1)%3) : F(f,(e+2)%3);
             Index j = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+2)%3) : F(f,(e+1)%3);
-            edge2faces[EMK(i,j)].push_back(f);
+            //edge2faces[EMK(i,j)].push_back(f);
           }
         }
         offending[f].push_back({f_other,result});
@@ -58,14 +60,10 @@ namespace igl
         typename DerivedVB,
         typename DerivedFB,
         typename DerivedIF,
-        typename DerivedVVA,
-        typename DerivedFFA,
-        typename DerivedJA,
-        typename DerivedIMA,
-        typename DerivedVVB,
-        typename DerivedFFB,
-        typename DerivedJB,
-        typename DerivedIMB>
+        typename DerivedVVAB,
+        typename DerivedFFAB,
+        typename DerivedJAB,
+        typename DerivedIMAB>
       static IGL_INLINE bool intersect_other_helper(
         const Eigen::PlainObjectBase<DerivedVA> & VA,
         const Eigen::PlainObjectBase<DerivedFA> & FA,
@@ -73,14 +71,10 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedFB> & FB,
         const RemeshSelfIntersectionsParam & params,
         Eigen::PlainObjectBase<DerivedIF> & IF,
-        Eigen::PlainObjectBase<DerivedVVA> & VVA,
-        Eigen::PlainObjectBase<DerivedFFA> & FFA,
-        Eigen::PlainObjectBase<DerivedJA>  & JA,
-        Eigen::PlainObjectBase<DerivedIMA> & IMA,
-        Eigen::PlainObjectBase<DerivedVVB> & VVB,
-        Eigen::PlainObjectBase<DerivedFFB> & FFB,
-        Eigen::PlainObjectBase<DerivedJB>  & JB,
-        Eigen::PlainObjectBase<DerivedIMB> & IMB)
+        Eigen::PlainObjectBase<DerivedVVAB> & VVAB,
+        Eigen::PlainObjectBase<DerivedFFAB> & FFAB,
+        Eigen::PlainObjectBase<DerivedJAB>  & JAB,
+        Eigen::PlainObjectBase<DerivedIMAB> & IMAB)
       {
 
         using namespace std;
@@ -139,7 +133,7 @@ namespace igl
         box_up(TA,A_boxes);
         box_up(TB,B_boxes);
         OffendingMap offendingA,offendingB;
-        EdgeMap edge2facesA,edge2facesB;
+        //EdgeMap edge2facesA,edge2facesB;
 
         std::list<int> lIF;
         const auto cb = [&](const Box &a, const Box &b) -> void
@@ -163,8 +157,8 @@ namespace igl
             {
               CGAL::Object result = CGAL::intersection(A,B);
 
-              push_result(FA,fa,fb,result,offendingA,edge2facesA);
-              push_result(FB,fb,fa,result,offendingB,edge2facesB);
+              push_result(FA,fa,fb,result,offendingA);
+              push_result(FB,fb,fa,result,offendingB);
             }
           }
         };
@@ -202,8 +196,38 @@ namespace igl
         }
         if(!params.detect_only)
         {
-          remesh_intersections(VA,FA,TA,offendingA,edge2facesA,VVA,FFA,JA,IMA);
-          remesh_intersections(VB,FB,TB,offendingB,edge2facesB,VVB,FFB,JB,IMB);
+          // Obsolete, now remesh_intersections expects a single mesh
+          // remesh_intersections(VA,FA,TA,offendingA,VVA,FFA,JA,IMA);
+          // remesh_intersections(VB,FB,TB,offendingB,VVB,FFB,JB,IMB);
+          // Combine mesh and offending maps
+          DerivedVA VAB(VA.rows()+VB.rows(),3);
+          VAB<<VA,VB;
+          DerivedFA FAB(FA.rows()+FB.rows(),3);
+          FAB<<FA,(FB.array()+VA.rows());
+          Triangles TAB;
+          TAB.reserve(TA.size()+TB.size());
+          TAB.insert(TAB.end(),TA.begin(),TA.end());
+          TAB.insert(TAB.end(),TB.begin(),TB.end());
+          OffendingMap offending;
+          //offending.reserve(offendingA.size() + offendingB.size());
+          for (const auto itr : offendingA)
+          {
+            // Remap offenders in FB to FAB
+            auto offenders = itr.second;
+            for(auto & offender : offenders)
+            {
+              offender.first += FA.rows();
+            }
+            offending[itr.first] = offenders;
+          }
+          for (const auto itr : offendingB)
+          {
+            // Store offenders for FB according to place in FAB
+            offending[FA.rows() + itr.first] = itr.second;
+          }
+
+          remesh_intersections(
+            VAB,FAB,TAB,offending,params.stitch_all,VVAB,FFAB,JAB,IMAB);
         }
 
         return IF.rows() > 0;
@@ -218,38 +242,30 @@ template <
   typename DerivedVB,
   typename DerivedFB,
   typename DerivedIF,
-  typename DerivedVVA,
-  typename DerivedFFA,
-  typename DerivedJA,
-  typename DerivedIMA,
-  typename DerivedVVB,
-  typename DerivedFFB,
-  typename DerivedJB,
-  typename DerivedIMB>
+  typename DerivedVVAB,
+  typename DerivedFFAB,
+  typename DerivedJAB,
+  typename DerivedIMAB>
 IGL_INLINE bool igl::copyleft::cgal::intersect_other(
-  const Eigen::PlainObjectBase<DerivedVA> & VA,
-  const Eigen::PlainObjectBase<DerivedFA> & FA,
-  const Eigen::PlainObjectBase<DerivedVB> & VB,
-  const Eigen::PlainObjectBase<DerivedFB> & FB,
-  const RemeshSelfIntersectionsParam & params,
-  Eigen::PlainObjectBase<DerivedIF> & IF,
-  Eigen::PlainObjectBase<DerivedVVA> & VVA,
-  Eigen::PlainObjectBase<DerivedFFA> & FFA,
-  Eigen::PlainObjectBase<DerivedJA>  & JA,
-  Eigen::PlainObjectBase<DerivedIMA> & IMA,
-  Eigen::PlainObjectBase<DerivedVVB> & VVB,
-  Eigen::PlainObjectBase<DerivedFFB> & FFB,
-  Eigen::PlainObjectBase<DerivedJB>  & JB,
-  Eigen::PlainObjectBase<DerivedIMB> & IMB)
+    const Eigen::PlainObjectBase<DerivedVA> & VA,
+    const Eigen::PlainObjectBase<DerivedFA> & FA,
+    const Eigen::PlainObjectBase<DerivedVB> & VB,
+    const Eigen::PlainObjectBase<DerivedFB> & FB,
+    const RemeshSelfIntersectionsParam & params,
+    Eigen::PlainObjectBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedVVAB> & VVAB,
+    Eigen::PlainObjectBase<DerivedFFAB> & FFAB,
+    Eigen::PlainObjectBase<DerivedJAB>  & JAB,
+    Eigen::PlainObjectBase<DerivedIMAB> & IMAB)
 {
   if(params.detect_only)
   {
     return intersect_other_helper<CGAL::Epick>
-      (VA,FA,VB,FB,params,IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+      (VA,FA,VB,FB,params,IF,VVAB,FFAB,JAB,IMAB);
   }else
   {
     return intersect_other_helper<CGAL::Epeck>
-      (VA,FA,VB,FB,params,IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+      (VA,FA,VB,FB,params,IF,VVAB,FFAB,JAB,IMAB);
   }
 }
 
@@ -261,11 +277,11 @@ IGL_INLINE bool igl::copyleft::cgal::intersect_other(
   const bool first_only,
   Eigen::MatrixXi & IF)
 {
-  Eigen::MatrixXd VVA,VVB;
-  Eigen::MatrixXi FFA,FFB;
-  Eigen::VectorXi JA,JB,IMA,IMB;
+  Eigen::MatrixXd VVAB;
+  Eigen::MatrixXi FFAB;
+  Eigen::VectorXi JAB,IMAB;
   return intersect_other(
-    VA,FA,VB,FB,{true,first_only},IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+    VA,FA,VB,FB,{true,first_only},IF,VVAB,FFAB,JAB,IMAB);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/copyleft/cgal/intersect_other.h

@@ -39,28 +39,21 @@ namespace igl
       // Outputs:
       //   IF  #intersecting face pairs by 2 list of intersecting face pairs,
       //     indexing FA and FB
-      //   VVA  #VVA by 3 list of vertex positions
-      //   FFA  #FFA by 3 list of triangle indices into VVA
-      //   JA  #FFA list of indices into FA denoting birth triangle
-      //   IMA  #VVA list of indices into VVA of unique vertices.
-      //   VVB  #VVB by 3 list of vertex positions
-      //   FFB  #FFB by 3 list of triangle indices into VVB
-      //   JB  #FFB list of indices into FB denoting birth triangle
-      //   IMB  #VVB list of indices into VVB of unique vertices.
+      //   VVAB  #VVAB by 3 list of vertex positions
+      //   FFAB  #FFAB by 3 list of triangle indices into VVA
+      //   JAB  #FFAB list of indices into [FA;FB] denoting birth triangle
+      //   IMAB  #VVAB list of indices stitching duplicates (resulting from
+      //     mesh intersections) together
       template <
         typename DerivedVA,
         typename DerivedFA,
         typename DerivedVB,
         typename DerivedFB,
         typename DerivedIF,
-        typename DerivedVVA,
-        typename DerivedFFA,
-        typename DerivedJA,
-        typename DerivedIMA,
-        typename DerivedVVB,
-        typename DerivedFFB,
-        typename DerivedJB,
-        typename DerivedIMB>
+        typename DerivedVVAB,
+        typename DerivedFFAB,
+        typename DerivedJAB,
+        typename DerivedIMAB>
       IGL_INLINE bool intersect_other(
         const Eigen::PlainObjectBase<DerivedVA> & VA,
         const Eigen::PlainObjectBase<DerivedFA> & FA,
@@ -68,14 +61,10 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedFB> & FB,
         const RemeshSelfIntersectionsParam & params,
         Eigen::PlainObjectBase<DerivedIF> & IF,
-        Eigen::PlainObjectBase<DerivedVVA> & VVA,
-        Eigen::PlainObjectBase<DerivedFFA> & FFA,
-        Eigen::PlainObjectBase<DerivedJA>  & JA,
-        Eigen::PlainObjectBase<DerivedIMA> & IMA,
-        Eigen::PlainObjectBase<DerivedVVB> & VVB,
-        Eigen::PlainObjectBase<DerivedFFB> & FFB,
-        Eigen::PlainObjectBase<DerivedJB>  & JB,
-        Eigen::PlainObjectBase<DerivedIMB> & IMB);
+        Eigen::PlainObjectBase<DerivedVVAB> & VVAB,
+        Eigen::PlainObjectBase<DerivedFFAB> & FFAB,
+        Eigen::PlainObjectBase<DerivedJAB>  & JAB,
+        Eigen::PlainObjectBase<DerivedIMAB> & IMAB);
       // Legacy wrapper for detect only using common types.
       //
       // Inputs:

include/igl/copyleft/cgal/point_mesh_squared_distance.cpp

@@ -8,14 +8,21 @@
 #include "point_mesh_squared_distance.h"
 #include "mesh_to_cgal_triangle_list.h"
 
-template <typename Kernel>
+template <
+  typename Kernel,
+  typename DerivedP,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedC>
 IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
-  const Eigen::MatrixXd & P,
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::VectorXd & sqrD,
-  Eigen::VectorXi & I,
-  Eigen::MatrixXd & C)
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+        Eigen::PlainObjectBase<DerivedI> & I,
+        Eigen::PlainObjectBase<DerivedC> & C)
 {
   using namespace std;
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
@@ -29,10 +36,10 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
   return point_mesh_squared_distance(P,tree,T,sqrD,I,C);
 }
 
-template <typename Kernel>
+template <typename Kernel, typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance_precompute(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   CGAL::AABB_tree<
     CGAL::AABB_traits<Kernel, 
       CGAL::AABB_triangle_primitive<Kernel, 
@@ -77,9 +84,14 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance_precompute(
   tree.closest_point_and_primitive(Point_3(0,0,0));
 }
 
-template <typename Kernel>
+template <
+  typename Kernel,
+  typename DerivedP,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedC>
 IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
-  const Eigen::MatrixXd & P,
+  const Eigen::PlainObjectBase<DerivedP> & P,
   const CGAL::AABB_tree<
     CGAL::AABB_traits<Kernel, 
       CGAL::AABB_triangle_primitive<Kernel, 
@@ -88,9 +100,9 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
     >
   > & tree,
   const std::vector<CGAL::Triangle_3<Kernel> > & T,
-  Eigen::VectorXd & sqrD,
-  Eigen::VectorXi & I,
-  Eigen::MatrixXd & C)
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C)
 {
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
   typedef typename std::vector<Triangle_3>::iterator Iterator;
@@ -110,18 +122,14 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
     // Find closest point and primitive id
     Point_and_primitive_id pp = tree.closest_point_and_primitive(query);
     Point_3 closest_point = pp.first;
-    C(p,0) = CGAL::to_double(closest_point[0]);
-    C(p,1) = CGAL::to_double(closest_point[1]);
-    C(p,2) = CGAL::to_double(closest_point[2]);
-    sqrD(p) = CGAL::to_double((closest_point-query).squared_length());
+    assign_scalar(closest_point[0],C(p,0));
+    assign_scalar(closest_point[1],C(p,1));
+    assign_scalar(closest_point[2],C(p,2));
+    assign_scalar((closest_point-query).squared_length(),sqrD(p));
     I(p) = pp.second - T.begin();
   }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::copyleft::cgal::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
-template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/copyleft/cgal/point_mesh_squared_distance.h

@@ -33,24 +33,35 @@ namespace igl
       //
       // Known bugs: This only computes distances to triangles. So unreferenced
       // vertices and degenerate triangles (segments) are ignored.
-      template <typename Kernel>
+      template <
+        typename Kernel,
+        typename DerivedP,
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedsqrD,
+        typename DerivedI,
+        typename DerivedC>
       IGL_INLINE void point_mesh_squared_distance(
-        const Eigen::MatrixXd & P,
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & F,
-        Eigen::VectorXd & sqrD,
-        Eigen::VectorXi & I,
-        Eigen::MatrixXd & C);
+        const Eigen::PlainObjectBase<DerivedP> & P,
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+              Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+              Eigen::PlainObjectBase<DerivedI> & I,
+              Eigen::PlainObjectBase<DerivedC> & C);
       // Probably can do this in a way that we don't pass around `tree` and `T`
       //
       // Outputs:
       //   tree  CGAL's AABB tree
       //   T  list of CGAL triangles in order of F (for determining which was found
       //     in computation)
-      template <typename Kernel>
+      template <
+        typename Kernel,
+        typename DerivedV,
+        typename DerivedF
+        >
       IGL_INLINE void point_mesh_squared_distance_precompute(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & F,
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
         CGAL::AABB_tree<
           CGAL::AABB_traits<Kernel, 
             CGAL::AABB_triangle_primitive<Kernel, 
@@ -63,9 +74,14 @@ namespace igl
       //  see above
       // Outputs:
       //  see above
-      template <typename Kernel>
+      template <
+        typename Kernel,
+        typename DerivedP,
+        typename DerivedsqrD,
+        typename DerivedI,
+        typename DerivedC>
       IGL_INLINE void point_mesh_squared_distance(
-        const Eigen::MatrixXd & P,
+        const Eigen::PlainObjectBase<DerivedP> & P,
         const CGAL::AABB_tree<
           CGAL::AABB_traits<Kernel, 
             CGAL::AABB_triangle_primitive<Kernel, 
@@ -74,9 +90,9 @@ namespace igl
           >
         > & tree,
         const std::vector<CGAL::Triangle_3<Kernel> > & T,
-        Eigen::VectorXd & sqrD,
-        Eigen::VectorXi & I,
-        Eigen::MatrixXd & C);
+        Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+        Eigen::PlainObjectBase<DerivedI> & I,
+        Eigen::PlainObjectBase<DerivedC> & C);
     }
   }
 }

include/igl/copyleft/cgal/point_solid_signed_squared_distance.cpp

@@ -0,0 +1,50 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "point_solid_signed_squared_distance.h"
+#include "points_inside_component.h"
+#include "point_mesh_squared_distance.h"
+#include "../../list_to_matrix.h"
+#include "../../slice_mask.h"
+#include <vector>
+#include <Eigen/Core>
+
+template <
+  typename DerivedQ,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedD>
+IGL_INLINE void igl::copyleft::cgal::point_solid_signed_squared_distance(
+  const Eigen::PlainObjectBase<DerivedQ> & Q,
+  const Eigen::PlainObjectBase<DerivedVB> & VB,
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  // compute unsigned distances
+  Eigen::VectorXi I;
+  DerivedVB C;
+  point_mesh_squared_distance<CGAL::Epeck>(Q,VB,FB,D,I,C);
+  // Collect queries that have non-zero distance
+  Eigen::Array<bool,Eigen::Dynamic,1> NZ = D.array()!=0;
+  // Compute sign for non-zero distance queries
+  DerivedQ QNZ;
+  slice_mask(Q,NZ,1,QNZ);
+  Eigen::Array<bool,Eigen::Dynamic,1> DNZ;
+  igl::copyleft::cgal::points_inside_component(VB,FB,QNZ,DNZ);
+  // Apply sign to distances
+  DerivedD S = DerivedD::Zero(Q.rows(),1);
+  {
+    int k = 0;
+    for(int q = 0;q<Q.rows();q++)
+    {
+      if(NZ(q))
+      {
+        D(q) *= DNZ(k++) ? -1. : 1.;
+      }
+    }
+  }
+}

include/igl/copyleft/cgal/point_solid_signed_squared_distance.h

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_COPYLEFT_CGAL_POINT_SOLID_SIGNED_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_POINT_SOLID_SIGNED_SQUARED_DISTANCE_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // POINT_SOLID_SIGNED_SQUARED_DISTANCE Given a set of points (Q) and the
+      // boundary mesh (VB,FB) of a solid (as defined in [Zhou et al. 2016],
+      // determine the signed squared distance for each point q in Q so that d(q,B) is
+      // negative if inside and positive if outside.
+      //
+      // Inputs:
+      //   Q  #Q by 3 list of query point positions
+      //   VB  #VB by 3 list of mesh vertex positions of B
+      //   FB  #FB by 3 list of mesh triangle indices into VB
+      // Outputs:
+      //   D
+      template <
+        typename DerivedQ,
+        typename DerivedVB,
+        typename DerivedFB,
+        typename DerivedD>
+      IGL_INLINE void point_solid_signed_squared_distance(
+        const Eigen::PlainObjectBase<DerivedQ> & Q,
+        const Eigen::PlainObjectBase<DerivedVB> & VB,
+        const Eigen::PlainObjectBase<DerivedFB> & FB,
+        Eigen::PlainObjectBase<DerivedD> & D);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_solid_signed_squared_distance.cpp"
+#endif
+
+#endif

include/igl/copyleft/cgal/remesh_intersections.h

@@ -30,7 +30,6 @@ namespace igl
       //   offending #offending map taking face indices into F to pairs of order
       //     of first finding and list of intersection objects from all
       //     intersections
-      //   edge2faces  #edges <= #offending*3 to incident offending faces 
       //   stitch_all  if true, merge all vertices with thte same coordiante.
       // Outputs:
       //   VV  #VV by 3 list of vertex positions, if stitch_all = false then
@@ -58,9 +57,6 @@ namespace igl
           typename DerivedF::Index,
             std::vector<
             std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
-        const std::map<
-          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-          std::vector<typename DerivedF::Index> > & edge2faces,
         bool stitch_all,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
@@ -83,9 +79,6 @@ namespace igl
           typename DerivedF::Index,
             std::vector<
             std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
-        const std::map<
-          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-          std::vector<typename DerivedF::Index> > & edge2faces,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,

include/igl/copyleft/cgal/trim_with_solid.cpp

@@ -0,0 +1,106 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "trim_with_solid.h"
+#include "assign_scalar.h"
+#include "intersect_other.h"
+#include "point_solid_signed_squared_distance.h"
+
+#include "../../extract_manifold_patches.h"
+#include "../../list_to_matrix.h"
+#include "../../remove_unreferenced.h"
+#include "../../slice_mask.h"
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <vector>
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedD,
+  typename DerivedJ>
+IGL_INLINE void igl::copyleft::cgal::trim_with_solid(
+  const Eigen::PlainObjectBase<DerivedVA> & VA,
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<DerivedVB> & VB,
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  Eigen::PlainObjectBase<DerivedV> & Vd,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedD> & D,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  // resolve intersections using exact representation
+  typedef Eigen::Matrix<CGAL::Epeck::FT,Eigen::Dynamic,3> MatrixX3E;
+  typedef Eigen::Matrix<CGAL::Epeck::FT,Eigen::Dynamic,1> VectorXE;
+  typedef Eigen::Matrix<CGAL::Epeck::FT,1,3> RowVector3E;
+  MatrixX3E V;
+  Eigen::MatrixXi _1;
+  Eigen::VectorXi _2;
+  // Intersect A and B meshes and stitch together new faces
+  igl::copyleft::cgal::intersect_other(
+    VA,FA,VB,FB,{false,false,true},_1,V,F,J,_2);
+  // Partition result into manifold patches
+  Eigen::VectorXi P;
+  const size_t num_patches = igl::extract_manifold_patches(F,P);
+  // only keep faces from A
+  Eigen::Matrix<bool,Eigen::Dynamic,1> A = J.array()< FA.rows();
+  igl::slice_mask(Eigen::MatrixXi(F),A,1,F);
+  igl::slice_mask(Eigen::VectorXi(P),A,1,P);
+  igl::slice_mask(Eigen::VectorXi(J),A,1,J);
+  // Agregate representative query points for each patch
+  std::vector<bool> flag(num_patches);
+  std::vector<std::vector<CGAL::Epeck::FT> > vQ;
+  Eigen::VectorXi P2Q(num_patches);
+  for(int f = 0;f<P.rows();f++)
+  {
+    const auto p = P(f);
+    // if not yet processed this patch
+    if(!flag[p])
+    {
+      P2Q(p) = vQ.size();
+      std::vector<CGAL::Epeck::FT> q = {
+        (V(F(f,0),0)+ V(F(f,1),0)+ V(F(f,2),0))/3.,
+        (V(F(f,0),1)+ V(F(f,1),1)+ V(F(f,2),1))/3.,
+        (V(F(f,0),2)+ V(F(f,1),2)+ V(F(f,2),2))/3.};
+      vQ.emplace_back(q);
+      flag[p] = true;
+    }
+  }
+  MatrixX3E Q;
+  igl::list_to_matrix(vQ,Q);
+  VectorXE SP;
+  point_solid_signed_squared_distance(Q,VB,FB,SP);
+  Eigen::Matrix<bool,Eigen::Dynamic,1> DP = SP.array()>0;
+  // distribute flag to all faces
+  D.resize(F.rows());
+  for(int f = 0;f<F.rows();f++)
+  {
+    D(f) = DP(P2Q(P(f)));
+  }
+  Eigen::VectorXi _;
+  igl::remove_unreferenced(MatrixX3E(V),DerivedF(F),V,F,_);
+  const auto & assign = [](
+    const MatrixX3E & V, 
+    Eigen::PlainObjectBase<DerivedV> & Vd)
+  {
+    Vd.resize(V.rows(),3);
+    for(int v = 0;v<V.rows();v++)
+    {
+      for(int d = 0;d<3;d++) 
+      {
+        igl::copyleft::cgal::assign_scalar(V(v,d),Vd(v,d));
+      }
+    }
+  };
+  assign(V,Vd);
+}
+

include/igl/copyleft/cgal/trim_with_solid.h

@@ -0,0 +1,63 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_COPYLEFT_CGAL_TRIM_WITH_SOLID_H
+#define IGL_COPYLEFT_CGAL_TRIM_WITH_SOLID_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // TRIM_WITH_SOLID Given an arbitrary mesh (VA,FA) and the boundary mesh
+      // (VB,FB) of a solid (as defined in [Zhou et al. 2016]), Resolve intersections
+      // between A and B subdividing faces of A so that intersections with B exists
+      // only along edges and vertices (and coplanar faces). Then determine whether
+      // each of these faces is inside or outside of B. This can be used to extract
+      // the part of A inside or outside of B.
+      //
+      // Inputs:
+      //   VA  #VA by 3 list of mesh vertex positions of A
+      //   FA  #FA by 3 list of mesh triangle indices into VA
+      //   VB  #VB by 3 list of mesh vertex positions of B
+      //   FB  #FB by 3 list of mesh triangle indices into VB
+      // Outputs:
+      //   V  #V by 3 list of mesh vertex positions of output
+      //   F  #F by 3 list of mesh triangle indices into V
+      //   D  #F list of bools whether face is inside B
+      //   J  #F list of indices into FA revealing birth parent
+      //
+      template <
+        typename DerivedVA,
+        typename DerivedFA,
+        typename DerivedVB,
+        typename DerivedFB,
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedD,
+        typename DerivedJ>
+      IGL_INLINE void trim_with_solid(
+        const Eigen::PlainObjectBase<DerivedVA> & VA,
+        const Eigen::PlainObjectBase<DerivedFA> & FA,
+        const Eigen::PlainObjectBase<DerivedVB> & VB,
+        const Eigen::PlainObjectBase<DerivedFB> & FB,
+        Eigen::PlainObjectBase<DerivedV> & Vd,
+        Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<DerivedD> & D,
+        Eigen::PlainObjectBase<DerivedJ> & J);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "trim_with_solid.cpp"
+#endif
+#endif

include/igl/copyleft/tetgen/tetrahedralize.cpp

@@ -84,17 +84,17 @@ IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
     bool TV_rect = list_to_matrix(vTV,TV);
     if(!TV_rect)
     {
-      return false;
+      return 3;
     }
     bool TT_rect = list_to_matrix(vTT,TT);
     if(!TT_rect)
     {
-      return false;
+      return 3;
     }
     bool TF_rect = list_to_matrix(vTF,TF);
     if(!TF_rect)
     {
-      return false;
+      return 3;
     }
   }
   return e;

include/igl/edge_flaps.h

@@ -24,6 +24,8 @@ namespace igl
   //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
   //     e=(j->i)
   //   EI  #E by 2 list of edge flap corners (see above).
+  //
+  // TODO: This seems to be a duplicate of edge_topology.h
   IGL_INLINE void edge_flaps(
     const Eigen::MatrixXi & F,
     const Eigen::MatrixXi & E,

include/igl/edge_topology.h

@@ -16,11 +16,13 @@
 namespace igl 
 {
   // Initialize Edges and their topological relations
-  
+  //
   // Output:
   // EV  : #Ex2, Stores the edge description as pair of indices to vertices
   // FE : #Fx3, Stores the Triangle-Edge relation
   // EF : #Ex2: Stores the Edge-Triangle relation
+  //
+  // TODO: This seems to be a duplicate of edge_flaps.h
 template <typename DerivedV, typename DerivedF>
   IGL_INLINE void edge_topology(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/euler_characteristic.cpp

@@ -0,0 +1,32 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich <michaelrabinovich27@gmail.com@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 "euler_characteristic.h"
+
+#include <igl/edge_topology.h>
+
+template <typename Scalar, typename Index>
+IGL_INLINE int igl::euler_characteristic(
+  const Eigen::PlainObjectBase<Scalar> & V,
+  const Eigen::PlainObjectBase<Index> & F)
+{
+
+  int euler_v = V.rows();
+  Eigen::MatrixXi EV, FE, EF;
+  igl::edge_topology(V, F, EV, FE, EF);
+  int euler_e = EV.rows();
+  int euler_f = F.rows();
+
+  int euler_char = euler_v - euler_e + euler_f;
+  return euler_char;
+
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+#endif

include/igl/euler_characteristic.h

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich <michaelrabinovich27@gmail.com@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_EULER_CHARACTERISTIC_H
+#define IGL_EULER_CHARACTERISTIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <vector>
+namespace igl
+{
+  // Computes the Euler characteristic of a given mesh (V,F)
+  // Templates:
+  //   Scalar  should be a floating point number type
+  //   Index   should be an integer type
+  // Inputs:
+  //   V       #V by dim list of mesh vertex positions
+  //   F       #F by dim list of mesh faces (must be triangles)
+  // Outputs:
+  //   An int containing the Euler characteristic
+  template <typename Scalar, typename Index>
+  IGL_INLINE int euler_characteristic(
+    const Eigen::PlainObjectBase<Scalar> & V,
+    const Eigen::PlainObjectBase<Index> & F);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "euler_characteristic.cpp"
+#endif
+
+#endif

include/igl/extract_manifold_patches.cpp

@@ -1,18 +1,20 @@
 #include "extract_manifold_patches.h"
+#include "unique_edge_map.h"
 #include <cassert>
 #include <limits>
 #include <queue>
 
 template<
-typename DerivedF,
-typename DerivedEMAP,
-typename uE2EType,
-typename DerivedP>
+  typename DerivedF,
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename DerivedP>
 IGL_INLINE size_t igl::extract_manifold_patches(
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-        const std::vector<std::vector<uE2EType> >& uE2E,
-        Eigen::PlainObjectBase<DerivedP>& P) {
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  const std::vector<std::vector<uE2EType> >& uE2E,
+  Eigen::PlainObjectBase<DerivedP>& P) 
+{
     assert(F.cols() == 3);
     const size_t num_faces = F.rows();
 
@@ -67,6 +69,20 @@ IGL_INLINE size_t igl::extract_manifold_patches(
     return num_patches;
 }
 
+template<
+  typename DerivedF,
+  typename DerivedP>
+IGL_INLINE size_t igl::extract_manifold_patches(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedP>& P) 
+{
+  Eigen::MatrixXi E, uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<size_t> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  return igl::extract_manifold_patches(F, EMAP, uE2E, P);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 template unsigned long igl::extract_manifold_patches<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, 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> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/extract_manifold_patches.h

@@ -30,6 +30,12 @@ namespace igl {
       const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
       const std::vector<std::vector<uE2EType> >& uE2E,
       Eigen::PlainObjectBase<DerivedP>& P);
+    template <
+      typename DerivedF,
+      typename DerivedP>
+    IGL_INLINE size_t extract_manifold_patches(
+      const Eigen::PlainObjectBase<DerivedF>& F,
+      Eigen::PlainObjectBase<DerivedP>& P);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "extract_manifold_patches.cpp"

include/igl/flipped_triangles.cpp

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich <michaelrabinovich27@gmail.com@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 "flipped_triangles.h"
+
+#include "list_to_matrix.h"
+#include <vector>
+template <typename DerivedV, typename DerivedF, typename DerivedX>
+IGL_INLINE void igl::flipped_triangles(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedX> & X)
+{
+  assert(V.cols() == 2 && "V should contain 2D positions");
+  std::vector<typename DerivedX::Scalar> flip_idx;
+  for (int i = 0; i < F.rows(); i++) 
+  {
+    // https://www.cs.cmu.edu/~quake/robust.html
+    typedef Eigen::Matrix<typename DerivedV::Scalar,1,2> RowVector2S;
+    RowVector2S v1_n = V.row(F(i,0)); 
+    RowVector2S v2_n = V.row(F(i,1)); 
+    RowVector2S v3_n = V.row(F(i,2));
+    Eigen::Matrix<typename DerivedV::Scalar,3,3> T2_Homo;
+    T2_Homo.col(0) << v1_n(0),v1_n(1),1.;
+    T2_Homo.col(1) << v2_n(0),v2_n(1),1.;
+    T2_Homo.col(2) << v3_n(0),v3_n(1),1.;
+    double det = T2_Homo.determinant();
+    assert(det == det && "det should not be NaN");
+    if (det < 0) 
+    {
+      flip_idx.push_back(i);
+    }
+  }
+  igl::list_to_matrix(flip_idx,X);
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE Eigen::VectorXi igl::flipped_triangles(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F)
+{
+  Eigen::VectorXi X;
+  flipped_triangles(V,F,X);
+  return X;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template void igl::flipped_triangles<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::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> >&);
+#endif

include/igl/flipped_triangles.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich <michaelrabinovich27@gmail.com@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_FLIPPED_TRIANGLES_H
+#define IGL_FLIPPED_TRIANGLES_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+namespace igl
+{
+  // Finds the ids of the flipped triangles of the mesh V,F in the UV mapping uv
+  //
+  // Inputs:
+  //   V  #V by 2 list of mesh vertex positions
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  // Outputs:
+  //   X  #flipped list of containing the indices into F of the flipped triangles
+  // Wrapper with return type
+  template <typename DerivedV, typename DerivedF, typename DerivedX>
+  IGL_INLINE void flipped_triangles(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedX> & X);
+  template <typename Scalar, typename Index>
+  IGL_INLINE Eigen::VectorXi flipped_triangles(
+    const Eigen::PlainObjectBase<Scalar> & V,
+    const Eigen::PlainObjectBase<Index> & F);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flipped_triangles.cpp"
+#endif
+
+#endif

include/igl/get_seconds.cpp

@@ -6,24 +6,10 @@
 // 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 "get_seconds.h"
-// NULL for Linux
-#include <cstddef>
-
-#if _WIN32
-#  include <ctime>
+#include <chrono>
 IGL_INLINE double igl::get_seconds()
 {
-  // This does not work on mac os x with glut in the main loop
-  return double(clock())/CLOCKS_PER_SEC;
+  return 
+    std::chrono::duration<double>(
+      std::chrono::system_clock::now().time_since_epoch()).count();
 }
-#else
-#  include <sys/time.h>
-IGL_INLINE double igl::get_seconds()
-{
-  timeval time;
-  gettimeofday(&time, NULL);
-  return time.tv_sec + time.tv_usec / 1e6;
-  // This does not work on mac os x with glut in the main loop
-  //return double(clock())/CLOCKS_PER_SEC;
-}
-#endif

include/igl/harmonic.cpp

@@ -6,10 +6,15 @@
 // 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 "harmonic.h"
+#include "adjacency_matrix.h"
 #include "cotmatrix.h"
-#include "massmatrix.h"
+#include "diag.h"
 #include "invert_diag.h"
+#include "isdiag.h"
+#include "massmatrix.h"
 #include "min_quad_with_fixed.h"
+#include "speye.h"
+#include "sum.h"
 #include <Eigen/Sparse>
 
 template <
@@ -28,8 +33,7 @@ IGL_INLINE bool igl::harmonic(
 {
   using namespace Eigen;
   typedef typename DerivedV::Scalar Scalar;
-  typedef Matrix<Scalar,Dynamic,1> VectorXS;
-  SparseMatrix<Scalar> L,M,Mi;
+  SparseMatrix<Scalar> L,M;
   cotmatrix(V,F,L);
   switch(F.cols())
   {
@@ -39,18 +43,73 @@ IGL_INLINE bool igl::harmonic(
     case 4:
     default:
       massmatrix(V,F,MASSMATRIX_TYPE_BARYCENTRIC,M);
-      break;
+    break;
   }
+  return harmonic(L,M,b,bc,k,W);
+}
+
+template <
+  typename DerivedF,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedW>
+IGL_INLINE bool igl::harmonic(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<Derivedb> & b,
+  const Eigen::PlainObjectBase<Derivedbc> & bc,
+  const int k,
+  Eigen::PlainObjectBase<DerivedW> & W)
+{
+  using namespace Eigen;
+  typedef typename Derivedbc::Scalar Scalar;
+  SparseMatrix<Scalar> A;
+  adjacency_matrix(F,A);
+  // sum each row
+  SparseVector<Scalar> Asum;
+  sum(A,1,Asum);
+  // Convert row sums into diagonal of sparse matrix
+  SparseMatrix<Scalar> Adiag;
+  diag(Asum,Adiag);
+  SparseMatrix<Scalar> L = A-Adiag;
+  SparseMatrix<Scalar> M;
+  speye(L.rows(),M);
+  return harmonic(L,M,b,bc,k,W);
+}
+
+template <
+  typename DerivedL,
+  typename DerivedM,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedW>
+IGL_INLINE bool igl::harmonic(
+  const Eigen::SparseMatrix<DerivedL> & L,
+  const Eigen::SparseMatrix<DerivedM> & M,
+  const Eigen::PlainObjectBase<Derivedb> & b,
+  const Eigen::PlainObjectBase<Derivedbc> & bc,
+  const int k,
+  Eigen::PlainObjectBase<DerivedW> & W)
+{
+  const int n = L.rows();
+  assert(n == L.cols() && "L must be square");
+  assert(n == M.cols() && "M must be same size as L");
+  assert(n == M.rows() && "M must be square");
+  assert(igl::isdiag(M) && "Mass matrix should be diagonal");
+
+  Eigen::SparseMatrix<DerivedM> Mi;
   invert_diag(M,Mi);
-  SparseMatrix<Scalar> Q = -L;
+  Eigen::SparseMatrix<DerivedL> Q;
+  Q = -L;
   for(int p = 1;p<k;p++)
   {
     Q = (Q*Mi*-L).eval();
   }
-  const VectorXS B = VectorXS::Zero(V.rows(),1);
+  typedef DerivedL Scalar;
   min_quad_with_fixed_data<Scalar> data;
-  min_quad_with_fixed_precompute(Q,b,SparseMatrix<Scalar>(),true,data);
-  W.resize(V.rows(),bc.cols());
+  min_quad_with_fixed_precompute(Q,b,Eigen::SparseMatrix<Scalar>(),true,data);
+  W.resize(n,bc.cols());
+  typedef Eigen::Matrix<Scalar,Eigen::Dynamic,1> VectorXS;
+  const VectorXS B = VectorXS::Zero(n,1);
   for(int w = 0;w<bc.cols();w++)
   {
     const VectorXS bcw = bc.col(w);
@@ -64,9 +123,14 @@ IGL_INLINE bool igl::harmonic(
   return true;
 }
 
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template bool igl::harmonic<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<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template bool igl::harmonic<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::harmonic<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+
+template bool igl::harmonic<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<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> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template bool igl::harmonic<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<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> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::harmonic<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<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> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/harmonic.h

@@ -1,15 +1,16 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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_HARMONIC_H
 #define IGL_HARMONIC_H
 #include "igl_inline.h"
 #include <Eigen/Core>
-namespace igl 
+#include <Eigen/Sparse>
+namespace igl
 {
   // Compute k-harmonic weight functions "coordinates".
   //
@@ -36,7 +37,52 @@ namespace igl
     const Eigen::PlainObjectBase<Derivedbc> & bc,
     const int k,
     Eigen::PlainObjectBase<DerivedW> & W);
+  // Compute harmonic map using uniform laplacian operator
+  //
+  // Inputs:
+  //   F  #F by simplex-size list of element indices
+  //   b  #b boundary indices into V
+  //   bc #b by #W list of boundary values
+  //   k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
+  // Outputs:
+  //   W  #V by #W list of weights
+  //
+  template <
+    typename DerivedF,
+    typename Derivedb,
+    typename Derivedbc,
+    typename DerivedW>
+  IGL_INLINE bool harmonic(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<Derivedb> & b,
+    const Eigen::PlainObjectBase<Derivedbc> & bc,
+    const int k,
+    Eigen::PlainObjectBase<DerivedW> & W);
+  // Compute a harmonic map using a given Laplacian and mass matrix
+  //
+  // Inputs:
+  //   L  #V by #V discrete (integrated) Laplacian  
+  //   M  #V by #V mass matrix
+  //   b  #b boundary indices into V
+  //   bc  #b by #W list of boundary values
+  //   k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
+  // Outputs:
+  //   W  #V by #V list of weights
+  template <
+    typename DerivedL,
+    typename DerivedM,
+    typename Derivedb,
+    typename Derivedbc,
+    typename DerivedW>
+  IGL_INLINE bool harmonic(
+    const Eigen::SparseMatrix<DerivedL> & L,
+    const Eigen::SparseMatrix<DerivedM> & M,
+    const Eigen::PlainObjectBase<Derivedb> & b,
+    const Eigen::PlainObjectBase<Derivedbc> & bc,
+    const int k,
+    Eigen::PlainObjectBase<DerivedW> & W);
 };
+
 #ifndef IGL_STATIC_LIBRARY
 #include "harmonic.cpp"
 #endif

include/igl/hsv_to_rgb.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 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 
+//
+// 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 "hsv_to_rgb.h"
 #include <cmath>
@@ -12,21 +12,21 @@
 template <typename T>
 IGL_INLINE void igl::hsv_to_rgb(const T * hsv, T * rgb)
 {
-  igl::hsv_to_rgb( 
+  igl::hsv_to_rgb(
       hsv[0],hsv[1],hsv[2],
       rgb[0],rgb[1],rgb[2]);
 }
 
 template <typename T>
-IGL_INLINE void igl::hsv_to_rgb( 
-  const T & h, const T & s, const T & v, 
+IGL_INLINE void igl::hsv_to_rgb(
+  const T & h, const T & s, const T & v,
   T & r, T & g, T & b)
 {
   // From medit
   double f,p,q,t,hh;
   int    i;
   hh = ((int)h % 360) / 60.;
-  i = (int)floor(hh);    /* largest int <= h     */
+  i = (int)std::floor(hh);    /* largest int <= h     */
   f = hh - i;                    /* fractional part of h */
   p = v * (1.0 - s);
   q = v * (1.0 - (s * f));

include/igl/isdiag.cpp

@@ -0,0 +1,32 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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 "isdiag.h"
+
+template <typename Scalar>
+IGL_INLINE bool igl::isdiag(const Eigen::SparseMatrix<Scalar> & A)
+{
+  // Iterate over outside of A
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
+    {
+      if(it.row() != it.col() && it.value()!=0)
+      {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::isdiag<double>(class Eigen::SparseMatrix<double,0,int> const &);
+#endif

include/igl/isdiag.h

@@ -0,0 +1,26 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_ISDIAG_H
+#define IGL_ISDIAG_H
+#include "igl_inline.h"
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Determine if a given matrix is diagonal: all non-zeros lie on the
+  // main diagonal.
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  // Returns true iff and only if the matrix is diagonal.
+  template <typename Scalar>
+  IGL_INLINE bool isdiag(const Eigen::SparseMatrix<Scalar> & A);
+};
+#ifndef IGL_STATIC_LIBRARY
+#  include "isdiag.cpp"
+#endif
+#endif

include/igl/loop.cpp

@@ -0,0 +1,152 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Oded Stein <oded.stein@columbia.edu>
+//
+// 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 "loop.h"
+
+#include <igl/adjacency_list.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/unique.h>
+
+#include <vector>
+
+namespace igl
+{
+    
+    IGL_INLINE void loop(const int n_verts,
+                         const Eigen::MatrixXi& F,
+                         Eigen::SparseMatrix<double>& S,
+                         Eigen::MatrixXi& newF)
+    {
+        
+        typedef Eigen::SparseMatrix<double> SparseMat;
+        typedef Eigen::Triplet<double> Triplet_t;
+        
+        //Ref. https://graphics.stanford.edu/~mdfisher/subdivision.html
+        //Heavily borrowing from igl::upsample
+        
+        Eigen::MatrixXi FF, FFi;
+        triangle_triangle_adjacency(F, FF, FFi);
+        std::vector<std::vector<int>> adjacencyList;
+        adjacency_list(F, adjacencyList, true);
+        
+        //Compute the number and positions of the vertices to insert (on edges)
+        Eigen::MatrixXi NI = Eigen::MatrixXi::Constant(FF.rows(), FF.cols(), -1);
+        Eigen::MatrixXi NIdoubles = Eigen::MatrixXi::Zero(FF.rows(), FF.cols());
+        Eigen::VectorXi vertIsOnBdry = Eigen::VectorXi::Zero(n_verts);
+        int counter = 0;
+        for(int i=0; i<FF.rows(); ++i)
+        {
+            for(int j=0; j<3; ++j)
+            {
+                if(NI(i,j) == -1)
+                {
+                    NI(i,j) = counter;
+                    NIdoubles(i,j) = 0;
+                    if (FF(i,j) != -1) {
+                        //If it is not a boundary
+                        NI(FF(i,j), FFi(i,j)) = counter;
+                        NIdoubles(i,j) = 1;
+                    } else {
+                        //Mark boundary vertices for later
+                        vertIsOnBdry(F(i,j)) = 1;
+                        vertIsOnBdry(F(i,(j+1)%3)) = 1;
+                    }
+                    ++counter;
+                }
+            }
+        }
+        
+        const int& n_odd = n_verts;
+        const int& n_even = counter;
+        const int n_newverts = n_odd + n_even;
+        
+        //Construct vertex positions
+        std::vector<Triplet_t> tripletList;
+        for(int i=0; i<n_odd; ++i) {
+            //Old vertices
+            const std::vector<int>& localAdjList = adjacencyList[i];
+            if(vertIsOnBdry(i)==1) {
+                //Boundary vertex
+                tripletList.emplace_back(i, localAdjList.front(), 1./8.);
+                tripletList.emplace_back(i, localAdjList.back(), 1./8.);
+                tripletList.emplace_back(i, i, 3./4.);
+            } else {
+                const int n = localAdjList.size();
+                const double dn = n;
+                double beta;
+                if(n==3)
+                    beta = 3./16.;
+                else
+                    beta = 3./8./dn;
+                for(int j=0; j<n; ++j)
+                    tripletList.emplace_back(i, localAdjList[j], beta);
+                tripletList.emplace_back(i, i, 1.-dn*beta);
+            }
+        }
+        for(int i=0; i<FF.rows(); ++i) {
+            //New vertices
+            for(int j=0; j<3; ++j) {
+                if(NIdoubles(i,j)==0) {
+                    if(FF(i,j)==-1) {
+                        //Boundary vertex
+                        tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 1./2.);
+                        tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+1)%3), 1./2.);
+                    } else {
+                        tripletList.emplace_back(NI(i,j) + n_odd, F(i,j), 3./8.);
+                        tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+1)%3), 3./8.);
+                        tripletList.emplace_back(NI(i,j) + n_odd, F(i, (j+2)%3), 1./8.);
+                        tripletList.emplace_back(NI(i,j) + n_odd, F(FF(i,j), (FFi(i,j)+2)%3), 1./8.);
+                    }
+                }
+            }
+        }
+        S.resize(n_newverts, n_verts);
+        S.setFromTriplets(tripletList.begin(), tripletList.end());
+        
+        // Build the new topology (Every face is replaced by four)
+        newF.resize(F.rows()*4, 3);
+        for(int i=0; i<F.rows();++i)
+        {
+            Eigen::VectorXi VI(6);
+            VI << F(i,0), F(i,1), F(i,2), NI(i,0) + n_odd, NI(i,1) + n_odd, NI(i,2) + n_odd;
+            
+            Eigen::VectorXi f0(3), f1(3), f2(3), f3(3);
+            f0 << VI(0), VI(3), VI(5);
+            f1 << VI(1), VI(4), VI(3);
+            f2 << VI(3), VI(4), VI(5);
+            f3 << VI(4), VI(2), VI(5);
+            
+            newF.row((i*4)+0) = f0;
+            newF.row((i*4)+1) = f1;
+            newF.row((i*4)+2) = f2;
+            newF.row((i*4)+3) = f3;
+        }
+        
+    }
+    
+    
+    IGL_INLINE void loop(const Eigen::MatrixXd& V,
+                         const Eigen::MatrixXi& F,
+                         Eigen::MatrixXd& newV,
+                         Eigen::MatrixXi& newF,
+                         const int number_of_subdivs)
+    {
+        typedef Eigen::SparseMatrix<double> SparseMat;
+        typedef Eigen::Triplet<double> Triplet_t;
+        
+        newV = V;
+        newF = F;
+        for(int i=0; i<number_of_subdivs; ++i) {
+            Eigen::MatrixXi tempF = newF;
+            SparseMat S;
+            loop(newV.rows(), tempF, S, newF);
+            newV = S*newV;
+        }
+    }
+    
+}

include/igl/loop.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Oded Stein <oded.stein@columbia.edu>
+//
+// 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_LOOP_H
+#define IGL_LOOP_H
+
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+    // LOOP Given the triangle mesh [V, F], where n_verts = V.rows(), computes newV and a sparse matrix S s.t. [newV, newF] is the subdivided mesh where newV = S*V.
+    //
+    // Inputs:
+    //  n_verts an integer (number of mesh vertices)
+    //  F an m by 3 matrix of integers of triangle faces
+    // Outputs:
+    //  S a sparse matrix (will become the subdivision matrix)
+    //  newF a matrix containing the new faces
+    IGL_INLINE void loop(const int n_verts,
+                         const Eigen::MatrixXi& F,
+                         Eigen::SparseMatrix<double>& S,
+                         Eigen::MatrixXi& newF);
+    
+    // LOOP Given the triangle mesh [V, F], computes number_of_subdivs steps of loop subdivision and outputs the new mesh [newV, newF]
+    //
+    // Inputs:
+    //  V an n by 3 matrix of vertices
+    //  F an m by 3 matrix of integers of triangle faces
+    //  number_of_subdivs an integer that specifies how many subdivision steps to do
+    // Outputs:
+    //  newV a matrix containing the new vertices
+    //  newF a matrix containing the new faces
+    IGL_INLINE void loop(const Eigen::MatrixXd& V,
+                         const Eigen::MatrixXi& F,
+                         Eigen::MatrixXd& newV,
+                         Eigen::MatrixXi& newF,
+                         const int number_of_subdivs = 1);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "loop.cpp"
+#endif
+
+#endif

include/igl/lscm.cpp

@@ -63,7 +63,7 @@ IGL_INLINE bool igl::lscm(
   V_uv.resize(V.rows(),2);
   for (unsigned i=0;i<V_uv.cols();++i)
   {
-    V_uv.col(i) = W_flat.block(V_uv.rows()*i,0,V_uv.rows(),1);
+    V_uv.col(V_uv.cols()-i-1) = W_flat.block(V_uv.rows()*i,0,V_uv.rows(),1);
   }
   return true;
 }

include/igl/min_quad_with_fixed.cpp

@@ -94,7 +94,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
 
   SparseMatrix<T> Auu;
   slice(A,data.unknown,data.unknown,Auu);
-  assert(Auu.size() > 0 && "There should be at least one unknown.");
+  assert(Auu.size() != 0 && Auu.rows() > 0 && "There should be at least one unknown.");
 
   // Positive definiteness is *not* determined, rather it is given as a
   // parameter

include/igl/parallel_for.h

@@ -123,7 +123,12 @@ inline bool igl::parallel_for(
   // Estimate number of threads in the pool
   // http://ideone.com/Z7zldb
   const static size_t sthc = std::thread::hardware_concurrency();
-  const size_t nthreads = loop_size<min_parallel?0:(sthc==0?8:sthc);
+  const size_t nthreads = 
+#ifdef IGL_PARALLEL_FOR_FORCE_SERIAL
+    0;
+#else
+    loop_size<min_parallel?0:(sthc==0?8:sthc);
+#endif
   if(nthreads==0)
   {
     // serial

include/igl/parula.cpp

@@ -69,4 +69,5 @@ template void igl::parula<double>(double, double*);
 template void igl::parula<double>(double, double&, double&, double&);
 template void igl::parula<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::parula<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::parula<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, double, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/segment_segment_intersect.cpp

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Francisca Gil Ureta <gilureta@cs.nyu.edu>
+//
+// 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 "segment_segment_intersect.h"
+
+#include <Eigen/Geometry>
+
+template<typename DerivedSource, typename DerivedDir>
+IGL_INLINE bool igl::segments_intersect(
+        const Eigen::PlainObjectBase <DerivedSource> &p,
+        const Eigen::PlainObjectBase <DerivedDir> &r,
+        const Eigen::PlainObjectBase <DerivedSource> &q,
+        const Eigen::PlainObjectBase <DerivedDir> &s,
+        double &a_t,
+        double &a_u,
+        double eps
+)
+{
+    // http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
+    // Search intersection between two segments
+    // p + t*r :  t \in [0,1]
+    // q + u*s :  u \in [0,1]
+
+    // p + t * r = q + u * s  // x s
+    // t(r x s) = (q - p) x s
+    // t = (q - p) x s / (r x s)
+
+    // (r x s) ~ 0 --> directions are parallel, they will never cross
+    Eigen::RowVector3d rxs = r.cross(s);
+    if (rxs.norm() <= eps)
+        return false;
+
+    int sign;
+
+    double u;
+    // u = (q − p) × r / (r × s)
+    Eigen::RowVector3d u1 = (q - p).cross(r);
+    sign = ((u1.dot(rxs)) > 0) ? 1 : -1;
+    u = u1.norm() / rxs.norm();
+    u = u * sign;
+
+    if ((u - 1.) > eps || u < -eps)
+        return false;
+
+    double t;
+    // t = (q - p) x s / (r x s)
+    Eigen::RowVector3d t1 = (q - p).cross(s);
+    sign = ((t1.dot(rxs)) > 0) ? 1 : -1;
+    t = t1.norm() / rxs.norm();
+    t = t * sign;
+
+    if (t < -eps || fabs(t) < eps)
+        return false;
+
+    a_t = t;
+    a_u = u;
+
+    return true;
+};
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::segments_intersect<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, double&, double&, double);
+#endif

include/igl/segment_segment_intersect.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Francisca Gil Ureta <gilureta@cs.nyu.edu>
+//
+// 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_SEGMENT_SEGMENT_INTERSECT_H
+#define IGL_SEGMENT_SEGMENT_INTERSECT_H
+
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+
+    // Determine whether two line segments A,B intersect
+    // A: p + t*r :  t \in [0,1]
+    // B: q + u*s :  u \in [0,1]
+    // Inputs:
+    //   p  3-vector origin of segment A
+    //   r  3-vector direction of segment A
+    //   q  3-vector origin of segment B
+    //   s  3-vector direction of segment B
+    //  eps precision
+    // Outputs:
+    //   t  scalar point of intersection along segment A, t \in [0,1]
+    //   u  scalar point of intersection along segment B, u \in [0,1]
+    // Returns true if intersection
+    template<typename DerivedSource, typename DerivedDir>
+    IGL_INLINE bool segments_intersect(
+            const Eigen::PlainObjectBase <DerivedSource> &p,
+            const Eigen::PlainObjectBase <DerivedDir> &r,
+            const Eigen::PlainObjectBase <DerivedSource> &q,
+            const Eigen::PlainObjectBase <DerivedDir> &s,
+            double &t,
+            double &u,
+            double eps = 1e-6
+    );
+
+}
+#ifndef IGL_STATIC_LIBRARY
+#   include "segment_segment_intersect.cpp"
+#endif
+#endif //IGL_SEGMENT_SEGMENT_INTERSECT_H

include/igl/sortrows.cpp

@@ -62,24 +62,46 @@ IGL_INLINE void igl::sortrows(
   using namespace std;
   using namespace Eigen;
   // Resize output
-  Y.resize(X.rows(),X.cols());
-  IX.resize(X.rows(),1);
-  for(int i = 0;i<X.rows();i++)
+  const size_t num_rows = X.rows();
+  const size_t num_cols = X.cols();
+  Y.resize(num_rows,num_cols);
+  IX.resize(num_rows,1);
+  for(int i = 0;i<num_rows;i++)
   {
     IX(i) = i;
   }
-  std::sort(
-    IX.data(),
-    IX.data()+IX.size(),
-    igl::IndexRowLessThan<const Eigen::PlainObjectBase<DerivedX> & >(X));
-  // if not ascending then reverse
-  if(!ascending)
-  {
-    std::reverse(IX.data(),IX.data()+IX.size());
+  if (ascending) {
+    auto index_less_than = [&X, num_cols](size_t i, size_t j) {
+      for (size_t c=0; c<num_cols; c++) {
+        if (X.coeff(i, c) < X.coeff(j, c)) return true;
+        else if (X.coeff(j,c) < X.coeff(i,c)) return false;
+      }
+      return false;
+    };
+      std::sort(
+        IX.data(),
+        IX.data()+IX.size(),
+        index_less_than
+        );
+  } else {
+    auto index_greater_than = [&X, num_cols](size_t i, size_t j) {
+      for (size_t c=0; c<num_cols; c++) {
+        if (X.coeff(i, c) > X.coeff(j, c)) return true;
+        else if (X.coeff(j,c) > X.coeff(i,c)) return false;
+      }
+      return false;
+    };
+      std::sort(
+        IX.data(),
+        IX.data()+IX.size(),
+        index_greater_than
+        );
   }
-  for(int i = 0;i<X.rows();i++)
-  {
-    Y.row(i) = X.row(IX(i));
+  for (size_t j=0; j<num_cols; j++) {
+      for(int i = 0;i<num_rows;i++)
+      {
+          Y(i,j) = X(IX(i), j);
+      }
   }
 }
 

include/igl/streamlines.cpp

@@ -0,0 +1,166 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Francisca Gil Ureta <gilureta@cs.nyu.edu>
+//
+// 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 "edge_topology.h"
+#include "sort_vectors_ccw.h"
+#include "streamlines.h"
+#include "per_face_normals.h"
+#include "polyvector_field_matchings.h"
+#include "segment_segment_intersect.h"
+#include "triangle_triangle_adjacency.h"
+#include "barycenter.h"
+#include "slice.h"
+
+#include <Eigen/Geometry>
+
+
+
+IGL_INLINE void igl::streamlines_init(
+        const Eigen::MatrixXd V,
+        const Eigen::MatrixXi F,
+        const Eigen::MatrixXd &temp_field,
+        const bool treat_as_symmetric,
+        StreamlineData &data,
+        StreamlineState &state,
+        double percentage
+){
+    using namespace Eigen;
+    using namespace std;
+
+    igl::edge_topology(V, F, data.E, data.F2E, data.E2F);
+    igl::triangle_triangle_adjacency(F, data.TT);
+
+    // prepare vector field
+    // --------------------------
+    int half_degree = temp_field.cols() / 3;
+    int degree = treat_as_symmetric ? half_degree * 2 : half_degree;
+    data.degree = degree;
+
+    Eigen::MatrixXd FN;
+    Eigen::VectorXi order, inv_order_unused;
+    Eigen::RowVectorXd sorted;
+
+    igl::per_face_normals(V, F, FN);
+    data.field.setZero(F.rows(), degree * 3);
+    for (unsigned i = 0; i < F.rows(); ++i){
+        const Eigen::RowVectorXd &n = FN.row(i);
+        Eigen::RowVectorXd temp(1, degree * 3);
+        if (treat_as_symmetric)
+            temp << temp_field.row(i), -temp_field.row(i);
+        else
+            temp = temp_field.row(i);
+        igl::sort_vectors_ccw(temp, n, order, true, sorted, false, inv_order_unused);
+
+        // project vectors to tangent plane
+        for (int j = 0; j < degree; ++j)
+        {
+            Eigen::RowVector3d pd = sorted.segment(j * 3, 3);
+            pd = (pd - (n.dot(pd)) * n).normalized();
+            data.field.block(i, j * 3, 1, 3) = pd;
+        }
+    }
+    Eigen::VectorXd curl;
+    igl::polyvector_field_matchings(data.field, V, F, false, data.match_ab, data.match_ba, curl);
+
+    // create seeds for tracing
+    // --------------------------
+    Eigen::VectorXi samples;
+    int nsamples;
+
+    nsamples = percentage * F.rows();
+    Eigen::VectorXd r;
+    r.setRandom(nsamples, 1);
+    r = (1 + r.array()) / 2.;
+    samples = (r.array() * F.rows()).cast<int>();
+    data.nsample = nsamples;
+
+    Eigen::MatrixXd BC, BC_sample;
+    igl::barycenter(V, F, BC);
+    igl::slice(BC, samples, 1, BC_sample);
+
+    // initialize state for tracing vector field
+
+    state.start_point = BC_sample.replicate(degree,1);
+    state.end_point = state.start_point;
+
+    state.current_face = samples.replicate(1, degree);
+
+    state.current_direction.setZero(nsamples, degree);
+    for (int i = 0; i < nsamples; ++i)
+        for (int j = 0; j < degree; ++j)
+            state.current_direction(i, j) = j;
+
+}
+
+IGL_INLINE void igl::streamlines_next(
+        const Eigen::MatrixXd V,
+        const Eigen::MatrixXi F,
+        const StreamlineData & data,
+        StreamlineState & state
+){
+    using namespace Eigen;
+    using namespace std;
+
+    int degree = data.degree;
+    int nsample = data.nsample;
+
+    state.start_point = state.end_point;
+
+    for (int i = 0; i < degree; ++i)
+    {
+        for (int j = 0; j < nsample; ++j)
+        {
+            int f0 = state.current_face(j,i);
+            if (f0 == -1) // reach boundary
+                continue;
+            int m0 = state.current_direction(j, i);
+
+            // the starting point of the vector
+            const Eigen::RowVector3d &p = state.start_point.row(j + nsample * i);
+            // the direction where we are trying to go
+            const Eigen::RowVector3d &r = data.field.block(f0, 3 * m0, 1, 3);
+
+
+            // new state,
+            int f1, m1;
+
+            for (int k = 0; k < 3; ++k)
+            {
+                f1 = data.TT(f0, k);
+
+                // edge vertices
+                const Eigen::RowVector3d &q = V.row(F(f0, k));
+                const Eigen::RowVector3d &qs = V.row(F(f0, (k + 1) % 3));
+                // edge direction
+                Eigen::RowVector3d s = qs - q;
+
+                double u;
+                double t;
+                if (igl::segments_intersect(p, r, q, s, t, u))
+                {
+                    // point on next face
+                    state.end_point.row(j + nsample * i) = p + t * r;
+                    state.current_face(j,i) = f1;
+
+                    // matching direction on next face
+                    int e1 = data.F2E(f0, k);
+                    if (data.E2F(e1, 0) == f0)
+                        m1 = data.match_ab(e1, m0);
+                    else
+                        m1 = data.match_ba(e1, m0);
+
+                    state.current_direction(j, i) = m1;
+                    break;
+                }
+
+            }
+
+
+        }
+    }
+}

include/igl/streamlines.h

@@ -0,0 +1,88 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Francisca Gil Ureta <gilureta@cs.nyu.edu>
+//
+// 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_STREAMLINES_H
+#define IGL_STREAMLINES_H
+
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+    struct StreamlineData
+    {
+        Eigen::MatrixXi TT;         //  #F by #3 adjacent matrix
+        Eigen::MatrixXi E;          //  #E by #3
+        Eigen::MatrixXi F2E;        //  #Fx3, Stores the Triangle-Edge relation
+        Eigen::MatrixXi E2F;        //  #Ex2, Stores the Edge-Triangle relation
+        Eigen::MatrixXd field;      //  #F by 3N list of the 3D coordinates of the per-face vectors
+                                    //      (N degrees stacked horizontally for each triangle)
+        Eigen::MatrixXi match_ab;   //  #E by N matrix, describing for each edge the matching a->b, where a
+                                    //      and b are the faces adjacent to the edge (i.e. vector #i of
+                                    //      the vector set in a is matched to vector #mab[i] in b)
+        Eigen::MatrixXi match_ba;   //  #E by N matrix, describing the inverse relation to match_ab
+        int nsample;                //  #S, number of sample points
+        int degree;                 //  #N, degrees of the vector field
+    };
+
+    struct StreamlineState
+    {
+        Eigen::MatrixXd start_point;        //  #N*S by 3 starting points of segment (stacked vertically for each degree)
+        Eigen::MatrixXd end_point;          //  #N*S by 3 endpoints points of segment (stacked vertically for each degree)
+        Eigen::MatrixXi current_face;       //  #S by N face indices (stacked horizontally for each degree)
+        Eigen::MatrixXi current_direction;  //  #S by N field direction indices (stacked horizontally for each degree)
+
+    };
+
+
+    // Given a mesh and a field the function computes the /data/ necessary for tracing the field'
+    // streamlines, and creates the initial /state/ for the tracing.
+    // Inputs:
+    //   V             #V by 3 list of mesh vertex coordinates
+    //   F             #F by 3 list of mesh faces
+    //   temp_field    #F by 3n list of the 3D coordinates of the per-face vectors
+    //                    (n-degrees stacked horizontally for each triangle)
+    //   treat_as_symmetric
+    //              if true, adds n symmetry directions to the field (N = 2n). Else N = n
+    //   percentage    [0-1] percentage of faces sampled
+    // Outputs:
+    //   data          struct containing topology information of the mesh and field
+    //   state         struct containing the state of the tracing
+    IGL_INLINE void streamlines_init(
+            const Eigen::MatrixXd V,
+            const Eigen::MatrixXi F,
+            const Eigen::MatrixXd &temp_field,
+            const bool treat_as_symmetric,
+            StreamlineData &data,
+            StreamlineState &state,
+            double percentage = 0.3
+
+    );
+
+    // The function computes the next state for each point in the sample
+    //   V             #V by 3 list of mesh vertex coordinates
+    //   F             #F by 3 list of mesh faces
+    //   data          struct containing topology information
+    //   state         struct containing the state of the tracing
+    IGL_INLINE void streamlines_next(
+            const Eigen::MatrixXd V,
+            const Eigen::MatrixXi F,
+            const StreamlineData & data,
+            StreamlineState & state
+
+    );
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "streamlines.cpp"
+#endif
+
+#endif

include/igl/unique.cpp

@@ -12,7 +12,6 @@
 #include "sortrows.h"
 #include "list_to_matrix.h"
 #include "matrix_to_list.h"
-#include "get_seconds.h"
 
 #include <algorithm>
 #include <iostream>
@@ -213,21 +212,32 @@ IGL_INLINE void igl::unique_rows(
   sortrows(A,true,sortA,IM);
 
 
-  vector<int> vIA(sortA.rows());
-  for(int i=0;i<(int)sortA.rows();i++)
+  const int num_rows = sortA.rows();
+  const int num_cols = sortA.cols();
+  vector<int> vIA(num_rows);
+  for(int i=0;i<num_rows;i++)
   {
     vIA[i] = i;
   }
+
+  auto index_equal = [&sortA, &num_cols](const size_t i, const size_t j) {
+    for (size_t c=0; c<num_cols; c++) {
+      if (sortA.coeff(i,c) != sortA.coeff(j,c))
+        return false;
+    }
+    return true;
+  };
   vIA.erase(
     std::unique(
     vIA.begin(),
     vIA.end(),
-    igl::IndexRowEquals<const Eigen::PlainObjectBase<DerivedA> &>(sortA)),vIA.end());
+    index_equal
+    ),vIA.end());
 
   IC.resize(A.rows(),1);
   {
     int j = 0;
-    for(int i = 0;i<(int)sortA.rows();i++)
+    for(int i = 0;i<num_rows;i++)
     {
       if(sortA.row(vIA[j]) != sortA.row(i))
       {
@@ -236,10 +246,11 @@ IGL_INLINE void igl::unique_rows(
       IC(IM(i,0),0) = j;
     }
   }
-  C.resize(vIA.size(),A.cols());
-  IA.resize(vIA.size(),1);
+  const int unique_rows = vIA.size();
+  C.resize(unique_rows,A.cols());
+  IA.resize(unique_rows,1);
   // Reindex IA according to IM
-  for(int i = 0;i<(int)vIA.size();i++)
+  for(int i = 0;i<unique_rows;i++)
   {
     IA(i,0) = IM(vIA[i],0);
     C.row(i) = A.row(IA(i,0));

include/igl/viewer/Viewer.cpp

@@ -84,7 +84,7 @@ static int global_KMod = 0;
 
 static void glfw_mouse_press(GLFWwindow* window, int button, int action, int modifier)
 {
-  bool tw_used = 
+  bool tw_used =
 #ifdef IGL_VIEWER_WITH_NANOGUI
     __viewer->screen->mouseButtonCallbackEvent(button,action,modifier);
 #else
@@ -261,6 +261,7 @@ namespace viewer
     ngui->addVariable("Show vertex labels", core.show_vertid);
     ngui->addVariable("Show faces labels", core.show_faceid);
 
+    screen->setVisible(true);
     screen->performLayout();
 #endif
 
@@ -275,33 +276,6 @@ namespace viewer
 
   IGL_INLINE Viewer::Viewer()
   {
-    // This mess is to help debug problems arising when compiling
-    // libiglviewer.a with(without) IGL_STATIC_LIBRARY defined and including
-    // Viewer.h without(with) IGL_STATIC_LIBRARY defined.
-#ifdef IGL_STATIC_LIBRARY
-    std::cout<<"igl_with_nanogui_defined_consistently: "<<igl_with_nanogui_defined_consistently()<<std::endl;
-    std::cout<<"igl_with_nanogui_defined_at_compile: "<<  igl_with_nanogui_defined_at_compile()  <<std::endl;
-    std::cout<<"igl_with_nanogui_defined_at_include: "<<  igl_with_nanogui_defined_at_include()  <<std::endl;
-    // First try to first assert
-    assert(igl_with_nanogui_defined_consistently() && 
-      "Must compile and include with IGL_VIEWER_WITH_NANOGUI defined consistently");
-#ifdef NDEBUG
-    // Second print warning since it's hopeless that this will run if wrong.
-    if(!igl_with_nanogui_defined_consistently())
-    {
-      std::cerr<<
-        "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl<<std::endl<<
-        "WARNING: Seems that IGL_WITH_NANOGUI " <<
-        (igl_with_nanogui_defined_at_compile() ? "was" : "was not") <<
-        " defined"<<std::endl<<"during compilation of Viewer.cpp and "<<
-        (igl_with_nanogui_defined_at_include() ? "was" : "was not") <<
-        " defined"<<std::endl<<"during inclusion of Viewer.h"<<std::endl <<
-        "You're about to get some nasty memory errors."<<std::endl<<std::endl<<
-        "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<std::endl;
-    }
-#endif
-#endif
-
 #ifdef IGL_VIEWER_WITH_NANOGUI
     ngui = nullptr;
     screen = nullptr;
@@ -629,11 +603,11 @@ namespace viewer
       center = data.V.colwise().sum()/data.V.rows();
     }
 
-    Eigen::Vector3f coord = 
+    Eigen::Vector3f coord =
       igl::project(
-        Eigen::Vector3f(center(0),center(1),center(2)), 
-        (core.view * core.model).eval(), 
-        core.proj, 
+        Eigen::Vector3f(center(0),center(1),center(2)),
+        (core.view * core.model).eval(),
+        core.proj,
         core.viewport);
     down_mouse_z = coord[2];
     down_rotation = core.trackball_angle;
@@ -1044,20 +1018,5 @@ namespace viewer
     launch_shut();
     return EXIT_SUCCESS;
   }
-  IGL_INLINE bool Viewer::igl_with_nanogui_defined_at_compile()
-  {
-#ifdef IGL_VIEWER_WITH_NANOGUI
-    return true;
-#else
-    return false;
-#endif
-  }
-  IGL_INLINE bool Viewer::igl_with_nanogui_defined_consistently()
-  {
-    return 
-      igl_with_nanogui_defined_at_compile() == 
-      igl_with_nanogui_defined_at_include();
-  }
 } // end namespace
 }
-

include/igl/viewer/Viewer.h

@@ -153,28 +153,10 @@ namespace viewer
     void* callback_key_down_data;
     void* callback_key_up_data;
 
-  public:
-
-    static IGL_INLINE bool igl_with_nanogui_defined_at_compile();
-    static IGL_INLINE bool igl_with_nanogui_defined_consistently();
-
   public:
       EIGEN_MAKE_ALIGNED_OPERATOR_NEW
   };
 
-  bool igl_with_nanogui_defined_at_include();
-#ifndef IGL_VIEWER_VIEWER_CPP
-  bool igl_with_nanogui_defined_at_include()
-  {
-    // this must be inlined here.
-#ifdef IGL_VIEWER_WITH_NANOGUI
-    return true;
-#else
-    return false;
-#endif
-  }
-#endif
-
 } // end namespace
 } // end namespace
 

index.html

@@ -39,8 +39,8 @@ like MATLAB.</p>
 just include igl headers (e.g. <code>#include &lt;igl/cotmatrix.h&gt;</code>) and run. Each
 header file contains a single function (e.g. <code>igl/cotmatrix.h</code> contains
 <code>igl::cotmatrix()</code>). Most are tailored to operate on a generic triangle mesh
-stored in an n-by&#8211;3 matrix of vertex positions V and an m-by&#8211;3 matrix of
-triangle indices F.</p>
+stored in an n-by&#8211;3 matrix of vertex positions <code>V</code> and an m-by&#8211;3 matrix of
+triangle indices <code>F</code>.</p>
 
 <p><em>Optionally</em> the library may also be <a href="optional/">pre-compiled</a> into a statically
 linked library, for faster compile times with your projects. This only effects

python/CMakeLists.txt

@@ -48,7 +48,7 @@ endif()
 include_directories(${PYTHON_INCLUDE_DIR} include)
 
 ## include pybind
-include_directories(${PROJECT_SOURCE_DIR}/../external/pybind11/include)
+include_directories(${PROJECT_SOURCE_DIR}/../external/nanogui/ext/pybind11/include)
 
 ## include libigl
 option(LIBIGL_USE_STATIC_LIBRARY "Use LibIGL as static library" OFF)
@@ -66,6 +66,9 @@ option(LIBIGL_WITH_PNG              "Use PNG"            ON)
 option(LIBIGL_WITH_TETGEN           "Use Tetgen"         ON)
 option(LIBIGL_WITH_TRIANGLE         "Use Triangle"       ON)
 option(LIBIGL_WITH_XML              "Use XML"            ON)
+option(LIBIGL_WITH_PYTHON           "Use Python"         ON)
+option(LIBIGL_WITH_COPYLEFT         "Use Copyleft"       ON)
+
 
 if(LIBIGL_WITH_CGAL) # Do not remove or move this block, cgal strange build system fails without it
   find_package(CGAL REQUIRED)
@@ -112,12 +115,17 @@ if (LIBIGL_WITH_CGAL)
   list(APPEND SHARED_SOURCES "modules/copyleft/py_igl_cgal.cpp")
 endif ()
 
+if (LIBIGL_WITH_COPYLEFT)
+  add_definitions(-DPY_COPYLEFT)
+  list(APPEND SHARED_SOURCES "modules/copyleft/py_igl_copyleft.cpp")
+endif ()
+
 if (LIBIGL_WITH_PNG)
   add_definitions(-DPY_PNG)
   list(APPEND SHARED_SOURCES "modules/py_igl_png.cpp")
 endif ()
 
-
+#----
 ## Prepare the python library
 add_library(pyigl SHARED
   python_shared.cpp
@@ -132,6 +140,10 @@ set_target_properties(pyigl PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE
 
 target_link_libraries(pyigl ${LIBIGL_LIBRARIES} ${LIBIGL_EXTRA_LIBRARIES})
 
+# Copy the nanogui bindings
+#get_target_property(NANOGUI_LIB nanogui_python LOCATION)
+#file(COPY ${NANOGUI_LIB} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../)
+
 if (WIN32)
   if (MSVC)
     # Enforce size-based optimization and link time code generation on MSVC (~30% smaller binaries in experiments)
@@ -147,7 +159,6 @@ if (WIN32)
   # target_link_libraries(igl ${PYTHON_LIBRARY})
   target_link_libraries(pyigl ${PYTHON_LIBRARIES})
 
-
 elseif (UNIX)
   # It's quite common to have multiple copies of the same Python version
   # installed on one's system. E.g.: one copy from the OS and another copy
@@ -190,3 +201,15 @@ elseif (UNIX)
     endif()
   endif()
 endif()
+
+if (LIBIGL_WITH_PYTHON AND LIBIGL_WITH_NANOGUI)
+  # Copy the nanogui python lib after compilation
+  get_target_property(NANOGUI_LIB nanogui_python LOCATION)
+  add_custom_target(copy ALL)
+  add_custom_command(
+        TARGET copy
+        COMMAND ${CMAKE_COMMAND} -E copy ${NANOGUI_LIB} ${CMAKE_CURRENT_BINARY_DIR}/../
+        )
+  add_dependencies(copy pyigl)
+
+endif()

python/iglhelpers.py

@@ -4,6 +4,8 @@ import pyigl as igl
 
 def p2e(m):
     if isinstance(m, np.ndarray):
+        if not m.flags['C_CONTIGUOUS']:
+            raise TypeError("p2e only support C-contiguous order")
         if m.dtype.type == np.int32:
             return igl.eigen.MatrixXi(m)
         elif m.dtype.type == np.float64:
@@ -33,11 +35,11 @@ def p2e(m):
 
 def e2p(m):
     if isinstance(m, igl.eigen.MatrixXd):
-        return np.array(m, dtype='float64')
+        return np.array(m, dtype='float64', order='C')
     elif isinstance(m, igl.eigen.MatrixXi):
-        return np.array(m, dtype='int32')
+        return np.array(m, dtype='int32', order='C')
     elif isinstance(m, igl.eigen.MatrixXb):
-        return np.array(m, dtype='bool')
+        return np.array(m, dtype='bool', order='C')
     elif isinstance(m, igl.eigen.SparseMatrixd):
         coo = np.array(m.toCOO())
         I = coo[:, 0]

python/modules/copyleft/py_igl_copyleft.cpp

@@ -0,0 +1,20 @@
+//#include <Eigen/Geometry>
+//#include <Eigen/Dense>
+//#include <Eigen/Sparse>
+
+
+#include "../../python_shared.h"
+
+#include <igl/copyleft/marching_cubes.h>
+#include <igl/copyleft/swept_volume.h>
+
+
+void python_export_igl_copyleft(py::module &me) {
+
+  py::module m = me.def_submodule(
+    "copyleft", "Wrappers for libigl functions that are copyleft");
+
+  #include "../../py_igl/copyleft/py_marching_cubes.cpp"
+  #include "../../py_igl/copyleft/py_swept_volume.cpp"
+
+}

python/modules/py_igl_embree.cpp

@@ -6,6 +6,7 @@
 #include "../python_shared.h"
 
 #include <igl/embree/ambient_occlusion.h>
+#include <igl/embree/reorient_facets_raycast.h>
 
 
 void python_export_igl_embree(py::module &me) {
@@ -14,5 +15,6 @@ void python_export_igl_embree(py::module &me) {
     "embree", "Wrappers for libigl functions that use embree");
 
   #include "../py_igl/embree/py_ambient_occlusion.cpp"
+  #include "../py_igl/embree/py_reorient_facets_raycast.cpp"
 
 }

python/modules/py_igl_viewer.cpp

@@ -8,6 +8,10 @@
 #include <igl/viewer/ViewerData.h>
 #include <igl/viewer/OpenGL_state.h>
 #include <igl/serialize.h>
+#ifdef IGL_VIEWER_WITH_NANOGUI
+#include "../../../external/nanogui/include/nanogui/formhelper.h"
+#include "../../../external/nanogui/include/nanogui/screen.h"
+#endif
 
 void python_export_igl_viewer(py::module &m)
 {
@@ -312,6 +316,14 @@ py::class_<igl::viewer::ViewerCore> viewercore_class(me, "ViewerCore");
 // UI Enumerations
     py::class_<igl::viewer::Viewer> viewer_class(me, "Viewer");
 
+    #ifdef IGL_VIEWER_WITH_NANOGUI
+    py::object fh = (py::object) py::module::import("nanogui").attr("FormHelper");
+    py::class_<nanogui::FormHelper> form_helper_class(me, "FormHelper", fh);
+
+    py::object screen = (py::object) py::module::import("nanogui").attr("Screen");
+    py::class_<nanogui::Screen> screen_class(me, "Screen", screen);
+    #endif
+
     py::enum_<igl::viewer::Viewer::MouseButton>(viewer_class, "MouseButton")
         .value("Left", igl::viewer::Viewer::MouseButton::Left)
         .value("Middle", igl::viewer::Viewer::MouseButton::Middle)
@@ -323,6 +335,12 @@ py::class_<igl::viewer::ViewerCore> viewercore_class(me, "ViewerCore");
     .def_readwrite("data", &igl::viewer::Viewer::data)
     .def_readwrite("core", &igl::viewer::Viewer::core)
     .def_readwrite("opengl", &igl::viewer::Viewer::opengl)
+
+    #ifdef IGL_VIEWER_WITH_NANOGUI
+    .def_readwrite("ngui", &igl::viewer::Viewer::ngui)
+    .def_readwrite("screen", &igl::viewer::Viewer::screen)
+    #endif
+
     .def("launch", &igl::viewer::Viewer::launch, py::arg("resizable") = true, py::arg("fullscreen") = false)
     .def("launch_init", &igl::viewer::Viewer::launch_init, py::arg("resizable") = true, py::arg("fullscreen") = false)
     .def("launch_rendering", &igl::viewer::Viewer::launch_rendering, py::arg("loop") = true)

python/modules/py_vector.cpp

@@ -153,6 +153,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("trace", [](const Type &m) {return m.trace();})
         .def("norm", [](const Type &m) {return m.norm();})
         .def("squaredNorm", [](const Type &m) {return m.squaredNorm();})
+        .def("squaredMean", [](const Type &m) {return m.array().square().mean();})
 
         .def("minCoeff", [](const Type &m) {return m.minCoeff();} )
         .def("maxCoeff", [](const Type &m) {return m.maxCoeff();} )
@@ -179,6 +180,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("rowwiseMean", [](const Type &m) {return Type(m.rowwise().mean());} )
         .def("rowwiseNorm", [](const Type &m) {return Type(m.rowwise().norm());} )
         .def("rowwiseNormalized", [](const Type &m) {return Type(m.rowwise().normalized());} )
+        .def("rowwiseReverse", [](const Type &m) {return Type(m.rowwise().reverse());} )
         .def("rowwiseMinCoeff", [](const Type &m) {return Type(m.rowwise().minCoeff());} )
         .def("rowwiseMaxCoeff", [](const Type &m) {return Type(m.rowwise().maxCoeff());} )
 
@@ -187,6 +189,8 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("colwiseProd", [](const Type &m) {return Type(m.colwise().prod());} )
         .def("colwiseMean", [](const Type &m) {return Type(m.colwise().mean());} )
         .def("colwiseNorm", [](const Type &m) {return Type(m.colwise().norm());} )
+        .def("colwiseNormalized", [](const Type &m) {return Type(m.colwise().normalized());} )
+        .def("colwiseReverse", [](const Type &m) {return Type(m.colwise().reverse());} )
         .def("colwiseMinCoeff", [](const Type &m) {return Type(m.colwise().minCoeff());} )
         .def("colwiseMaxCoeff", [](const Type &m) {return Type(m.colwise().maxCoeff());} )
 
@@ -688,6 +692,25 @@ void python_export_vector(py::module &m) {
     .def("solve",[](const Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s, const Eigen::MatrixXd& rhs) { return Eigen::MatrixXd(s.solve(rhs)); })
     ;
 
+    py::class_<Eigen::Affine3d > affine3d(me, "Affine3d");
+
+    affine3d
+    .def(py::init<>())
+    .def("setIdentity",[](Eigen::Affine3d& a){
+        return a.setIdentity();
+    })
+    .def("rotate",[](Eigen::Affine3d& a, double angle, Eigen::MatrixXd axis) {
+        assert_is_Vector3("axis", axis);
+        return a.rotate(Eigen::AngleAxisd(angle, Eigen::Vector3d(axis)));
+    })
+    .def("translate",[](Eigen::Affine3d& a, Eigen::MatrixXd offset) {
+        assert_is_Vector3("offset", offset);
+        return a.translate(Eigen::Vector3d(offset));
+    })
+    .def("matrix", [](Eigen::Affine3d& a) -> Eigen::MatrixXd {
+        return Eigen::MatrixXd(a.matrix());
+    })
+    ;
     /* Bindings for Quaterniond*/
     //py::class_<Eigen::Quaterniond > quaterniond(me, "Quaterniond");
     //

python/py_doc.cpp

@@ -112,6 +112,7 @@ const char *__doc_igl_cat = R"igl_Qu8mg5v7(// Perform concatenation of a two mat
   // Outputs:
   //   C  output matrix
   //   )igl_Qu8mg5v7";
+const char *__doc_igl_collapse_edge = R"igl_Qu8mg5v7(See collapse_edge for the documentation.)igl_Qu8mg5v7";
 const char *__doc_igl_colon = R"igl_Qu8mg5v7(// Colon operator like matlab's colon operator. Enumerats values between low
   // and hi with step step.
   // Templates:
@@ -219,6 +220,45 @@ const char *__doc_igl_copyleft_comiso_nrosy = R"igl_Qu8mg5v7(// Generate a N-RoS
     // Outputs:
     //   R       #F by 3 the representative vectors of the interpolated field
     //   S       #V by 1 the singularity index for each vertex (0 = regular))igl_Qu8mg5v7";
+const char *__doc_igl_copyleft_marching_cubes = R"igl_Qu8mg5v7(// marching_cubes( values, points, x_res, y_res, z_res, vertices, faces )
+    //
+    // performs marching cubes reconstruction on the grid defined by values, and
+    // points, and generates vertices and faces
+    //
+    // Input:
+    //  values  #number_of_grid_points x 1 array -- the scalar values of an
+    //    implicit function defined on the grid points (<0 in the inside of the
+    //    surface, 0 on the border, >0 outside)
+    //  points  #number_of_grid_points x 3 array -- 3-D positions of the grid
+    //    points, ordered in x,y,z order:
+    //      points[index] = the point at (x,y,z) where :
+    //      x = (index % (xres -1),
+    //      y = (index / (xres-1)) %(yres-1),
+    //      z = index / (xres -1) / (yres -1) ).
+    //      where x,y,z index x, y, z dimensions
+    //      i.e. index = x + y*xres + z*xres*yres
+    //  xres  resolutions of the grid in x dimension
+    //  yres  resolutions of the grid in y dimension
+    //  zres  resolutions of the grid in z dimension
+    // Output:
+    //   vertices  #V by 3 list of mesh vertex positions
+    //   faces  #F by 3 list of mesh triangle indices
+    //)igl_Qu8mg5v7";
+const char *__doc_igl_copyleft_swept_volume = R"igl_Qu8mg5v7(// Compute the surface of the swept volume of a solid object with surface
+    // (V,F) mesh under going rigid motion.
+    // 
+    // Inputs:
+    //   V  #V by 3 list of mesh positions in reference pose
+    //   F  #F by 3 list of mesh indices into V
+    //   transform  function handle so that transform(t) returns the rigid
+    //     transformation at time t∈[0,1]
+    //   steps  number of time steps: steps=3 --> t∈{0,0.5,1}
+    //   grid_res  number of grid cells on the longest side containing the
+    //     motion (isolevel+1 cells will also be added on each side as padding)
+    //   isolevel  distance level to be contoured as swept volume
+    // Outputs:
+    //   SV  #SV by 3 list of mesh positions of the swept surface
+    //   SF  #SF by 3 list of mesh faces into SV)igl_Qu8mg5v7";
 const char *__doc_igl_copyleft_tetgen_tetrahedralize = R"igl_Qu8mg5v7(// Mesh the interior of a surface mesh (V,F) using tetgen
       //
       // Inputs:
@@ -338,12 +378,19 @@ const char *__doc_igl_edge_lengths = R"igl_Qu8mg5v7(// Constructs a list of leng
   //    or
   //   T  #T by 4 list of mesh elements (must be tets)
   // Outputs:
-  //   L  #F by {1|3|6} list of edge lengths 
+  //   L  #F by {1|3|6} list of edge lengths
   //     for edges, column of lengths
   //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
   //     for tets, columns correspond to edges
   //     [3 0],[3 1],[3 2],[1 2],[2 0],[0 1]
   //)igl_Qu8mg5v7";
+const char *__doc_igl_edge_topology = R"igl_Qu8mg5v7(// Initialize Edges and their topological relations
+  //
+  // Output:
+  // EV  : #Ex2, Stores the edge description as pair of indices to vertices
+  // FE : #Fx3, Stores the Triangle-Edge relation
+  // EF : #Ex2: Stores the Edge-Triangle relation
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_eigs = R"igl_Qu8mg5v7(See eigs for the documentation.)igl_Qu8mg5v7";
 const char *__doc_igl_embree_ambient_occlusion = R"igl_Qu8mg5v7(// Compute ambient occlusion per given point
     //
@@ -355,6 +402,23 @@ const char *__doc_igl_embree_ambient_occlusion = R"igl_Qu8mg5v7(// Compute ambie
     //    S  #P list of ambient occlusion values between 1 (fully occluded) and
     //      0 (not occluded)
     //)igl_Qu8mg5v7";
+const char *__doc_igl_embree_reorient_facets_raycast = R"igl_Qu8mg5v7(// Orient each component (identified by C) of a mesh (V,F) using ambient
+    // occlusion such that the front side is less occluded than back side, as
+    // described in "A Simple Method for Correcting Facet Orientations in
+    // Polygon Meshes Based on Ray Casting" [Takayama et al. 2014].
+    //
+    // Inputs:
+    //   V  #V by 3 list of vertex positions
+    //   F  #F by 3 list of triangle indices
+    //   rays_total  Total number of rays that will be shot
+    //   rays_minimum  Minimum number of rays that each patch should receive
+    //   facet_wise  Decision made for each face independently, no use of patches
+    //     (i.e., each face is treated as a patch)
+    //   use_parity  Use parity mode
+    //   is_verbose  Verbose output to cout
+    // Outputs:
+    //   I  #F list of whether face has been flipped
+    //   C  #F list of patch ID (output of bfs_orient > manifold patches))igl_Qu8mg5v7";
 const char *__doc_igl_find_cross_field_singularities = R"igl_Qu8mg5v7(// Inputs:
   //   V                #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F                #F by 3 eigen Matrix of face (quad) indices
@@ -402,12 +466,28 @@ const char *__doc_igl_gaussian_curvature = R"igl_Qu8mg5v7(// Compute discrete lo
   // Output:
   //   K  #V by 1 eigen Matrix of discrete gaussian curvature values
   //)igl_Qu8mg5v7";
-const char *__doc_igl_grad = R"igl_Qu8mg5v7(// Gradient of a scalar function defined on piecewise linear elements (mesh)
-  // is constant on each triangle i,j,k:
-  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
-  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
-  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
-  // 90 degrees
+const char *__doc_igl_get_seconds = R"igl_Qu8mg5v7(// Return the current time in seconds since program start
+  // 
+  // Example:
+  //    const auto & tictoc = []()
+  //    {
+  //      static double t_start = igl::get_seconds();
+  //      double diff = igl::get_seconds()-t_start;
+  //      t_start += diff;
+  //      return diff;
+  //    };
+  //    tictoc();
+  //    ... // part 1
+  //    cout<<"part 1: "<<tictoc()<<endl;
+  //    ... // part 2
+  //    cout<<"part 2: "<<tictoc()<<endl;
+  //    ... // etc)igl_Qu8mg5v7";
+const char *__doc_igl_grad = R"igl_Qu8mg5v7(// Gradient of a scalar function defined on piecewise linear elements (mesh)
+  // is constant on each triangle i,j,k:
+  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
+  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
+  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
+  // 90 degrees
   //)igl_Qu8mg5v7";
 const char *__doc_igl_harmonic = R"igl_Qu8mg5v7(// Compute k-harmonic weight functions "coordinates".
   //
@@ -421,12 +501,41 @@ const char *__doc_igl_harmonic = R"igl_Qu8mg5v7(// Compute k-harmonic weight fun
   // Outputs:
   //   W  #V by #W list of weights
   //)igl_Qu8mg5v7";
+const char *__doc_igl_hsv_to_rgb = R"igl_Qu8mg5v7(// Convert RGB to HSV
+  //
+  // Inputs:
+  //   h  hue value (degrees: [0,360])
+  //   s  saturation value ([0,1])
+  //   v  value value ([0,1])
+  // Outputs:
+  //   r  red value ([0,1]) 
+  //   g  green value ([0,1])
+  //   b  blue value ([0,1]))igl_Qu8mg5v7";
+const char *__doc_igl_internal_angles = R"igl_Qu8mg5v7(// Compute internal angles for a triangle mesh
+  //
+  // Inputs:
+  //   V  #V by dim eigen Matrix of mesh vertex nD positions
+  //   F  #F by poly-size eigen Matrix of face (triangle) indices
+  // Output:
+  //   K  #F by poly-size eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //
+  // Known Issues:
+  //   if poly-size ≠ 3 then dim must equal 3.)igl_Qu8mg5v7";
 const char *__doc_igl_invert_diag = R"igl_Qu8mg5v7(// Templates:
   //   T  should be a eigen sparse matrix primitive type like int or double
   // Inputs:
   //   X  an m by n sparse matrix
   // Outputs:
   //   Y  an m by n sparse matrix)igl_Qu8mg5v7";
+const char *__doc_igl_is_irregular_vertex = R"igl_Qu8mg5v7(// Determine if a vertex is irregular, i.e. it has more than 6 (triangles)
+  // or 4 (quads) incident edges. Vertices on the boundary are ignored.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3[4] list of triangle[quads] indices
+  // Returns #V vector of bools revealing whether vertices are singular
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_jet = R"igl_Qu8mg5v7(// JET like MATLAB's jet
   //
   // Inputs:
@@ -694,6 +803,12 @@ const char *__doc_igl_quad_planarity = R"igl_Qu8mg5v7(// Compute planarity of th
   // Output:
   //   P  #F by 1 eigen Matrix of mesh face (quad) planarities
   //)igl_Qu8mg5v7";
+const char *__doc_igl_randperm = R"igl_Qu8mg5v7(// Like matlab's randperm(n) but minus 1
+  //
+  // Inputs:
+  //   n  number of elements
+  // Outputs:
+  //   I  n list of rand permutation of 0:n-1)igl_Qu8mg5v7";
 const char *__doc_igl_readDMAT = R"igl_Qu8mg5v7(See readDMAT for the documentation.)igl_Qu8mg5v7";
 const char *__doc_igl_readMESH = R"igl_Qu8mg5v7(// load a tetrahedral volume mesh from a .mesh file
   //
@@ -873,6 +988,40 @@ const char *__doc_igl_sortrows = R"igl_Qu8mg5v7(// Act like matlab's [Y,I] = sor
   //     reference as X)
   //   I  m list of indices so that
   //     Y = X(I,:);)igl_Qu8mg5v7";
+const char *__doc_igl_streamlines_init = R"igl_Qu8mg5v7(    // Given a mesh and a field the function computes the /data/ necessary for tracing the field'
+    // streamlines, and creates the initial /state/ for the tracing.
+    // Inputs:
+    //   V             #V by 3 list of mesh vertex coordinates
+    //   F             #F by 3 list of mesh faces
+    //   temp_field    #F by 3n list of the 3D coordinates of the per-face vectors
+    //                    (n-degrees stacked horizontally for each triangle)
+    //   treat_as_symmetric
+    //              if true, adds n symmetry directions to the field (N = 2n). Else N = n
+    //   percentage    [0-1] percentage of faces sampled
+    // Outputs:
+    //   data          struct containing topology information of the mesh and field
+    //   state         struct containing the state of the tracing )igl_Qu8mg5v7";
+const char *__doc_igl_streamlines_next = R"igl_Qu8mg5v7( // The function computes the next state for each point in the sample
+    //   V             #V by 3 list of mesh vertex coordinates
+    //   F             #F by 3 list of mesh faces
+    //   data          struct containing topology information
+    //   state         struct containing the state of the tracing )igl_Qu8mg5v7";
+const char *__doc_igl_triangle_triangle_adjacency = R"igl_Qu8mg5v7(// Constructs the triangle-triangle adjacency matrix for a given
+  // mesh (V,F).
+  //
+  // Templates:
+  //   Scalar derived type of eigen matrix for V (e.g. derived from
+  //     MatrixXd)
+  //   Index  derived type of eigen matrix for F (e.g. derived from
+  //     MatrixXi)
+  // Inputs:
+  //   F  #F by simplex_size list of mesh faces (must be triangles)
+  // Outputs:
+  //   TT   #F by #3 adjacent matrix, the element i,j is the id of the triangle adjacent to the j edge of triangle i
+  //   TTi  #F by #3 adjacent matrix, the element i,j is the id of edge of the triangle TT(i,j) that is adjacent with triangle i
+  // NOTE: the first edge of a triangle is [0,1] the second [1,2] and the third [2,3].
+  //       this convention is DIFFERENT from cotmatrix_entries.h
+  // Known bug: this should not need to take V as input.)igl_Qu8mg5v7";
 const char *__doc_igl_triangle_triangulate = R"igl_Qu8mg5v7(// Triangulate the interior of a polygon using the triangle library.
     //
     // Inputs:
@@ -941,6 +1090,41 @@ const char *__doc_igl_upsample = R"igl_Qu8mg5v7(// Subdivide a mesh without movi
   //
   // Known issues:
   //   - assumes (V,F) is edge-manifold.)igl_Qu8mg5v7";
+const char *__doc_igl_winding_number = R"igl_Qu8mg5v7(// WINDING_NUMBER Compute the sum of solid angles of a triangle/tetrahedron
+  // described by points (vectors) V
+  //
+  // Templates:
+  //   dim  dimension of input
+  // Inputs:
+  //  V  n by 3 list of vertex positions
+  //  F  #F by 3 list of triangle indices, minimum index is 0
+  //  O  no by 3 list of origin positions
+  // Outputs:
+  //  S  no by 1 list of winding numbers
+  //
+  // 3d)igl_Qu8mg5v7";
+const char *__doc_igl_winding_number_3 = R"igl_Qu8mg5v7(// Inputs:
+  //   V  pointer to array containing #V by 3 vertex positions along rows,
+  //     given in column major order
+  //   n  number of mesh vertices
+  //   F  pointer to array containing #F by 3 face indices along rows,
+  //     given in column major order
+  //   m  number of faces
+  //   O  pointer to array containing #O by 3 query positions along rows,
+  //     given in column major order
+  //   no  number of origins
+  // Outputs:
+  //   S  no by 1 list of winding numbers)igl_Qu8mg5v7";
+const char *__doc_igl_winding_number_2 = R"igl_Qu8mg5v7(//// Only one evaluation origin
+  //template <typename DerivedF>
+  //IGL_INLINE void winding_number_3(
+  //  const double * V,
+  //  const int n,
+  //  const DerivedF * F,
+  //  const int m,
+  //  const double * O,
+  //  double * S);
+  // 2d)igl_Qu8mg5v7";
 const char *__doc_igl_writeMESH = R"igl_Qu8mg5v7(// save a tetrahedral volume mesh to a .mesh file
   //
   // Templates:

python/py_doc.h

@@ -7,6 +7,7 @@ extern const char *__doc_igl_barycentric_coordinates;
 extern const char *__doc_igl_boundary_facets;
 extern const char *__doc_igl_boundary_loop;
 extern const char *__doc_igl_cat;
+extern const char *__doc_igl_collapse_edge;
 extern const char *__doc_igl_colon;
 extern const char *__doc_igl_comb_cross_field;
 extern const char *__doc_igl_comb_frame_field;
@@ -14,6 +15,8 @@ extern const char *__doc_igl_compute_frame_field_bisectors;
 extern const char *__doc_igl_copyleft_cgal_mesh_boolean;
 extern const char *__doc_igl_copyleft_comiso_miq;
 extern const char *__doc_igl_copyleft_comiso_nrosy;
+extern const char *__doc_igl_copyleft_marching_cubes;
+extern const char *__doc_igl_copyleft_swept_volume;
 extern const char *__doc_igl_copyleft_tetgen_tetrahedralize;
 extern const char *__doc_igl_cotmatrix;
 extern const char *__doc_igl_covariance_scatter_matrix;
@@ -23,17 +26,23 @@ extern const char *__doc_igl_doublearea;
 extern const char *__doc_igl_doublearea_single;
 extern const char *__doc_igl_doublearea_quad;
 extern const char *__doc_igl_edge_lengths;
+extern const char *__doc_igl_edge_topology;
 extern const char *__doc_igl_eigs;
 extern const char *__doc_igl_embree_ambient_occlusion;
+extern const char *__doc_igl_embree_reorient_facets_raycast;
 extern const char *__doc_igl_find_cross_field_singularities;
 extern const char *__doc_igl_fit_rotations;
 extern const char *__doc_igl_fit_rotations_planar;
 extern const char *__doc_igl_fit_rotations_SSE;
 extern const char *__doc_igl_floor;
 extern const char *__doc_igl_gaussian_curvature;
+extern const char *__doc_igl_get_seconds;
 extern const char *__doc_igl_grad;
 extern const char *__doc_igl_harmonic;
+extern const char *__doc_igl_hsv_to_rgb;
+extern const char *__doc_igl_internal_angles;
 extern const char *__doc_igl_invert_diag;
+extern const char *__doc_igl_is_irregular_vertex;
 extern const char *__doc_igl_jet;
 extern const char *__doc_igl_local_basis;
 extern const char *__doc_igl_lscm;
@@ -56,6 +65,7 @@ extern const char *__doc_igl_point_mesh_squared_distance;
 extern const char *__doc_igl_polar_svd;
 extern const char *__doc_igl_principal_curvature;
 extern const char *__doc_igl_quad_planarity;
+extern const char *__doc_igl_randperm;
 extern const char *__doc_igl_readDMAT;
 extern const char *__doc_igl_readMESH;
 extern const char *__doc_igl_readOBJ;
@@ -71,10 +81,16 @@ extern const char *__doc_igl_slice_into;
 extern const char *__doc_igl_slice_mask;
 extern const char *__doc_igl_slice_tets;
 extern const char *__doc_igl_sortrows;
+extern const char *__doc_igl_streamlines_init;
+extern const char *__doc_igl_streamlines_next;
+extern const char *__doc_igl_triangle_triangle_adjacency;
 extern const char *__doc_igl_triangle_triangulate;
 extern const char *__doc_igl_unique;
 extern const char *__doc_igl_unique_rows;
 extern const char *__doc_igl_unproject_onto_mesh;
 extern const char *__doc_igl_upsample;
+extern const char *__doc_igl_winding_number;
+extern const char *__doc_igl_winding_number_3;
+extern const char *__doc_igl_winding_number_2;
 extern const char *__doc_igl_writeMESH;
 extern const char *__doc_igl_writeOBJ;

python/py_igl.cpp

@@ -14,6 +14,7 @@
 #include <igl/boundary_facets.h>
 #include <igl/boundary_loop.h>
 #include <igl/cat.h>
+#include <igl/collapse_edge.h>
 #include <igl/colon.h>
 #include <igl/comb_cross_field.h>
 #include <igl/comb_frame_field.h>
@@ -24,14 +25,19 @@
 #include <igl/cut_mesh_from_singularities.h>
 #include <igl/doublearea.h>
 #include <igl/edge_lengths.h>
+#include <igl/edge_topology.h>
 #include <igl/eigs.h>
 #include <igl/find_cross_field_singularities.h>
 #include <igl/fit_rotations.h>
 #include <igl/floor.h>
 #include <igl/gaussian_curvature.h>
+#include <igl/get_seconds.h>
 #include <igl/grad.h>
 #include <igl/harmonic.h>
+#include <igl/hsv_to_rgb.h>
+#include <igl/internal_angles.h>
 #include <igl/invert_diag.h>
+#include <igl/is_irregular_vertex.h>
 #include <igl/jet.h>
 #include <igl/local_basis.h>
 #include <igl/lscm.h>
@@ -49,6 +55,7 @@
 #include <igl/polar_svd.h>
 #include <igl/principal_curvature.h>
 #include <igl/quad_planarity.h>
+#include <igl/randperm.h>
 #include <igl/readDMAT.h>
 #include <igl/readMESH.h>
 #include <igl/readOBJ.h>
@@ -62,9 +69,12 @@
 #include <igl/slice_mask.h>
 #include <igl/slice_tets.h>
 #include <igl/sortrows.h>
+#include <igl/streamlines.h>
+#include <igl/triangle_triangle_adjacency.h>
 #include <igl/unique.h>
 #include <igl/unproject_onto_mesh.h>
 #include <igl/upsample.h>
+#include <igl/winding_number.h>
 #include <igl/writeMESH.h>
 #include <igl/writeOBJ.h>
 
@@ -83,6 +93,7 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_boundary_facets.cpp"
 #include "py_igl/py_boundary_loop.cpp"
 #include "py_igl/py_cat.cpp"
+#include "py_igl/py_collapse_edge.cpp"
 #include "py_igl/py_colon.cpp"
 #include "py_igl/py_comb_cross_field.cpp"
 #include "py_igl/py_comb_frame_field.cpp"
@@ -93,14 +104,19 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_cut_mesh_from_singularities.cpp"
 #include "py_igl/py_doublearea.cpp"
 #include "py_igl/py_edge_lengths.cpp"
+#include "py_igl/py_edge_topology.cpp"
 #include "py_igl/py_eigs.cpp"
 #include "py_igl/py_find_cross_field_singularities.cpp"
 #include "py_igl/py_fit_rotations.cpp"
 #include "py_igl/py_floor.cpp"
 #include "py_igl/py_gaussian_curvature.cpp"
+#include "py_igl/py_get_seconds.cpp"
 #include "py_igl/py_grad.cpp"
 #include "py_igl/py_harmonic.cpp"
+#include "py_igl/py_hsv_to_rgb.cpp"
+#include "py_igl/py_internal_angles.cpp"
 #include "py_igl/py_invert_diag.cpp"
+#include "py_igl/py_is_irregular_vertex.cpp"
 #include "py_igl/py_jet.cpp"
 #include "py_igl/py_local_basis.cpp"
 #include "py_igl/py_lscm.cpp"
@@ -118,6 +134,7 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_polar_svd.cpp"
 #include "py_igl/py_principal_curvature.cpp"
 #include "py_igl/py_quad_planarity.cpp"
+#include "py_igl/py_randperm.cpp"
 #include "py_igl/py_readDMAT.cpp"
 #include "py_igl/py_readMESH.cpp"
 #include "py_igl/py_readOBJ.cpp"
@@ -131,9 +148,12 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_slice_mask.cpp"
 #include "py_igl/py_slice_tets.cpp"
 #include "py_igl/py_sortrows.cpp"
+#include "py_igl/py_streamlines.cpp"
+#include "py_igl/py_triangle_triangle_adjacency.cpp"
 #include "py_igl/py_unique.cpp"
 #include "py_igl/py_unproject_onto_mesh.cpp"
 #include "py_igl/py_upsample.cpp"
+#include "py_igl/py_winding_number.cpp"
 #include "py_igl/py_writeMESH.cpp"
 #include "py_igl/py_writeOBJ.cpp"
 

python/py_igl/copyleft/py_marching_cubes.cpp

@@ -0,0 +1,21 @@
+
+
+m.def("marching_cubes", []
+(
+  const Eigen::MatrixXd& values,
+  const Eigen::MatrixXd& points,
+  const unsigned int x_res,
+  const unsigned int y_res,
+  const unsigned int z_res,
+  Eigen::MatrixXd& vertices,
+  Eigen::MatrixXi& faces
+)
+{
+  assert_is_VectorX("values", values);
+  Eigen::VectorXd valuesv;
+  if (values.size() != 0)
+    valuesv = values;
+  return igl::copyleft::marching_cubes(valuesv, points, x_res, y_res, z_res, vertices, faces);
+}, __doc_igl_copyleft_marching_cubes,
+py::arg("values"), py::arg("points"), py::arg("x_res"), py::arg("y_res"), py::arg("z_res"), py::arg("vertices"), py::arg("faces"));
+

python/py_igl/copyleft/py_swept_volume.cpp

@@ -0,0 +1,19 @@
+
+m.def("swept_volume", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const std::function<Eigen::Affine3d (const double)> & transform,
+  const size_t steps,
+  const size_t grid_res,
+  const size_t isolevel,
+  Eigen::MatrixXd& SV,
+  Eigen::MatrixXi& SF
+)
+{
+  return igl::copyleft::swept_volume(V, F, transform, steps, grid_res, isolevel, SV, SF);
+}, __doc_igl_copyleft_swept_volume,
+py::arg("V"), py::arg("F"), py::arg("transform"), py::arg("steps"), py::arg("grid_res"), py::arg("isolevel"), py::arg("SV"), py::arg("SF"));
+
+
+

python/py_igl/embree/py_reorient_facets_raycast.cpp

@@ -0,0 +1,37 @@
+
+
+m.def("reorient_facets_raycast", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  int rays_total,
+  int rays_minimum,
+  bool facet_wise,
+  bool use_parity,
+  bool is_verbose,
+  Eigen::MatrixXi& I,
+  Eigen::MatrixXi& C
+)
+{
+  Eigen::VectorXi Iv;
+  Eigen::VectorXi Cv;
+  igl::embree::reorient_facets_raycast(V, F, rays_total, rays_minimum, facet_wise, use_parity, is_verbose, Iv, Cv);
+  I = Iv;
+  C = Cv;
+}, __doc_igl_embree_reorient_facets_raycast,
+py::arg("V"), py::arg("F"), py::arg("rays_total"), py::arg("rays_minimum"), py::arg("facet_wise"), py::arg("use_parity"), py::arg("is_verbose"), py::arg("I"), py::arg("C"));
+
+m.def("reorient_facets_raycast", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXi& FF,
+  Eigen::MatrixXi& I
+)
+{
+  Eigen::VectorXi Iv;
+  igl::embree::reorient_facets_raycast(V, F, FF, Iv);
+  I = Iv;
+}, __doc_igl_embree_reorient_facets_raycast,
+py::arg("V"), py::arg("F"), py::arg("FF"), py::arg("I"));
+

python/py_igl/py_collapse_edge.cpp

@@ -0,0 +1,87 @@
+// COMPLETE BINDINGS ========================
+
+
+
+
+
+
+// INCOMPLETE BINDINGS ========================
+
+
+//m.def("collapse_edge", []
+//(
+//  int e,
+//  Eigen::RowVectorXd & p,
+//  Eigen::MatrixXd& V,
+//  Eigen::MatrixXi& F,
+//  Eigen::MatrixXi& E,
+//  Eigen::MatrixXi& EMAP,
+//  Eigen::MatrixXi& EF,
+//  Eigen::MatrixXi& EI,
+//  int & e1,
+//  int & e2,
+//  int & f1,
+//  int & f2
+//)
+//{
+//  return igl::collapse_edge(e, p, V, F, E, EMAP, EF, EI, e1, e2, f1, f2);
+//}, __doc_igl_collapse_edge,
+//py::arg("e"), py::arg("p"), py::arg("V"), py::arg("F"), py::arg("E"), py::arg("EMAP"), py::arg("EF"), py::arg("EI"), py::arg("e1"), py::arg("e2"), py::arg("f1"), py::arg("f2"));
+
+//m.def("collapse_edge", []
+//(
+//  int e,
+//  Eigen::RowVectorXd & p,
+//  Eigen::MatrixXd& V,
+//  Eigen::MatrixXi& F,
+//  Eigen::MatrixXi& E,
+//  Eigen::MatrixXi& EMAP,
+//  Eigen::MatrixXi& EF,
+//  Eigen::MatrixXi& EI
+//)
+//{
+//  return igl::collapse_edge(e, p, V, F, E, EMAP, EF, EI);
+//}, __doc_igl_collapse_edge,
+//py::arg("e"), py::arg("p"), py::arg("V"), py::arg("F"), py::arg("E"), py::arg("EMAP"), py::arg("EF"), py::arg("EI"));
+
+//m.def("collapse_edge", []
+//(
+//  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::MatrixXi& 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
+//)
+//{
+//  return igl::collapse_edge(cost_and_placement, V, F, E, EMAP, EF, EI, Q, Qit, C);
+//}, __doc_igl_collapse_edge,
+//py::arg("cost_and_placement"), py::arg("V"), py::arg("F"), py::arg("E"), py::arg("EMAP"), py::arg("EF"), py::arg("EI"), py::arg("Q"), py::arg("Qit"), py::arg("C"));
+
+//m.def("collapse_edge", []
+//(
+//  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::MatrixXi& 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
+//)
+//{
+//  return igl::collapse_edge(cost_and_placement, V, F, E, EMAP, EF, EI, Q, Qit, C, e, e1, e2, f1, f2);
+//}, __doc_igl_collapse_edge,
+//py::arg("cost_and_placement"), py::arg("V"), py::arg("F"), py::arg("E"), py::arg("EMAP"), py::arg("EF"), py::arg("EI"), py::arg("Q"), py::arg("Qit"), py::arg("C"), py::arg("e"), py::arg("e1"), py::arg("e2"), py::arg("f1"), py::arg("f2"));
+

python/py_igl/py_edge_topology.cpp

@@ -0,0 +1,13 @@
+m.def("edge_topology", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXi& EV,
+  Eigen::MatrixXi& FE,
+  Eigen::MatrixXi& EF
+)
+{
+  return igl::edge_topology(V, F, EV, FE, EF);
+}, __doc_igl_edge_lengths,
+py::arg("V"), py::arg("F"), py::arg("EV"), py::arg("FE"), py::arg("EF"));
+

python/py_igl/py_get_seconds.cpp

@@ -0,0 +1,6 @@
+m.def("get_seconds", []
+()
+{
+  return igl::get_seconds();
+}, __doc_igl_get_seconds);
+

python/py_igl/py_hsv_to_rgb.cpp

@@ -0,0 +1,44 @@
+// COMPLETE BINDINGS ========================
+
+
+m.def("hsv_to_rgb", []
+(
+  const Eigen::MatrixXd& H,
+  Eigen::MatrixXd& R
+)
+{
+  return igl::hsv_to_rgb(H, R);
+}, __doc_igl_hsv_to_rgb,
+py::arg("H"), py::arg("R"));
+
+
+
+
+
+// INCOMPLETE BINDINGS ========================
+
+
+//m.def("hsv_to_rgb", []
+//(
+//  T * hsv,
+//  T * rgb
+//)
+//{
+//  return igl::hsv_to_rgb(hsv, rgb);
+//}, __doc_igl_hsv_to_rgb,
+//py::arg("hsv"), py::arg("rgb"));
+
+//m.def("hsv_to_rgb", []
+//(
+//  T & h,
+//  T & s,
+//  T & v,
+//  T & r,
+//  T & g,
+//  T & b
+//)
+//{
+//  return igl::hsv_to_rgb(h, s, v, r, g, b);
+//}, __doc_igl_hsv_to_rgb,
+//py::arg("h"), py::arg("s"), py::arg("v"), py::arg("r"), py::arg("g"), py::arg("b"));
+

python/py_igl/py_internal_angles.cpp

@@ -0,0 +1,23 @@
+
+
+m.def("internal_angles", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXd& K
+)
+{
+  return igl::internal_angles(V, F, K);
+}, __doc_igl_internal_angles,
+py::arg("V"), py::arg("F"), py::arg("K"));
+
+m.def("internal_angles", []
+(
+  const Eigen::MatrixXd& L,
+  Eigen::MatrixXd& K
+)
+{
+  return igl::internal_angles(L, K);
+}, __doc_igl_internal_angles,
+py::arg("L"), py::arg("K"));
+

python/py_igl/py_is_irregular_vertex.cpp

@@ -0,0 +1,12 @@
+
+
+m.def("is_irregular_vertex", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F
+)
+{
+  return igl::is_irregular_vertex(V, F);
+}, __doc_igl_is_irregular_vertex,
+py::arg("V"), py::arg("F"));
+

python/py_igl/py_parula.cpp

@@ -8,6 +8,17 @@
 //}, __doc_igl_parula,
 //py::arg("f"), py::arg("rgb"));
 
+m.def("parula", []
+(
+const double f
+)
+{
+  double r, g, b;
+  igl::parula(f, r, g, b);
+  return std::make_tuple(r,g,b);
+}, __doc_igl_parula,
+py::arg("f"));
+
 m.def("parula", []
 (
   const double f,

python/py_igl/py_randperm.cpp

@@ -0,0 +1,10 @@
+m.def("randperm", []
+(
+  int n,
+  Eigen::MatrixXi& I
+)
+{
+  return igl::randperm(n, I);
+}, __doc_igl_randperm,
+py::arg("n"), py::arg("I"));
+

python/py_igl/py_slice_tets.cpp

@@ -6,15 +6,15 @@ m.def("slice_tets", []
   Eigen::MatrixXd& U,
   Eigen::MatrixXi& G,
   Eigen::MatrixXi& J,
-  Eigen::SparseMatrix<double> & BC
+  Eigen::SparseMatrix<double>& BC
 )
 {
   assert_is_VectorX("plane", plane);
-  Eigen::VectorXd pl;
+  Eigen::VectorXd planev;
   if (plane.size() != 0)
-    pl = plane;
+    planev = plane;
   Eigen::VectorXi Jv;
-  igl::slice_tets(V, T, pl, U, G, Jv, BC);
+  igl::slice_tets(V, T, planev, U, G, Jv, BC);
   J = Jv;
 }, __doc_igl_slice_tets,
 py::arg("V"), py::arg("T"), py::arg("plane"), py::arg("U"), py::arg("G"), py::arg("J"), py::arg("BC"));

python/py_igl/py_streamlines.cpp

@@ -0,0 +1,53 @@
+py::class_<igl::StreamlineData> StreamlineData(m, "StreamlineData");
+StreamlineData
+.def(py::init<>())
+.def_readwrite("TT", &igl::StreamlineData::TT)
+.def_readwrite("E", &igl::StreamlineData::E)
+.def_readwrite("F2E", &igl::StreamlineData::F2E)
+.def_readwrite("E2F", &igl::StreamlineData::E2F)
+.def_readwrite("field", &igl::StreamlineData::field)
+.def_readwrite("match_ab", &igl::StreamlineData::match_ab)
+.def_readwrite("match_ba", &igl::StreamlineData::match_ba)
+.def_readwrite("nsample", &igl::StreamlineData::nsample)
+.def_readwrite("degree", &igl::StreamlineData::degree)
+;
+
+py::class_<igl::StreamlineState> StreamlineState(m, "StreamlineState");
+StreamlineState
+.def(py::init<>())
+.def_readwrite("start_point", &igl::StreamlineState::start_point)
+.def_readwrite("end_point", &igl::StreamlineState::end_point)
+.def_readwrite("current_face", &igl::StreamlineState::current_face)
+.def_readwrite("current_direction", &igl::StreamlineState::current_direction)
+.def("copy", [](const igl::StreamlineState &m) { return igl::StreamlineState(m); })
+;
+
+m.def("streamlines_init", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& temp_field,
+  const bool treat_as_symmetric,
+  igl::StreamlineData &data,
+  igl::StreamlineState &state,
+  double percentage
+)
+{
+  return igl::streamlines_init(V, F, temp_field, treat_as_symmetric, data, state, percentage);
+
+},__doc_igl_streamlines_init,
+py::arg("V"), py::arg("F"), py::arg("temp_field"), py::arg("treat_as_symmetric"),
+py::arg("data"), py::arg("state"), py::arg("percentage")=0.3);
+
+m.def("streamlines_next", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const igl::StreamlineData &data,
+  igl::StreamlineState &state
+)
+{
+  return igl::streamlines_next(V, F, data, state);
+
+},__doc_igl_streamlines_next,
+py::arg("V"), py::arg("F"), py::arg("data"), py::arg("state"));

python/py_igl/py_triangle_triangle_adjacency.cpp

@@ -0,0 +1,21 @@
+
+m.def("triangle_triangle_adjacency", []
+(
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXi& TT,
+  Eigen::MatrixXi& TTi
+)
+{
+  return igl::triangle_triangle_adjacency(F, TT, TTi);
+}, __doc_igl_triangle_triangle_adjacency,
+py::arg("F"), py::arg("TT"), py::arg("TTi"));
+
+m.def("triangle_triangle_adjacency", []
+(
+  const Eigen::MatrixXi& F,
+  Eigen::MatrixXi& TT
+)
+{
+  return igl::triangle_triangle_adjacency(F, TT);
+}, __doc_igl_triangle_triangle_adjacency,
+py::arg("F"), py::arg("TT"));

python/py_igl/py_winding_number.cpp

@@ -0,0 +1,18 @@
+// COMPLETE BINDINGS ========================
+
+
+m.def("winding_number", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& O,
+  Eigen::MatrixXd& W
+)
+{
+  Eigen::VectorXd Wv;
+  igl::winding_number(V, F, O, Wv);
+  W = Wv;
+}, __doc_igl_winding_number,
+py::arg("V"), py::arg("F"), py::arg("O"), py::arg("W"));
+
+

python/python_shared.cpp

@@ -34,6 +34,10 @@ extern void python_export_igl_cgal(py::module &);
 extern void python_export_igl_png(py::module &);
 #endif
 
+#ifdef PY_COPYLEFT
+extern void python_export_igl_copyleft(py::module &);
+#endif
+
 PYBIND11_PLUGIN(pyigl) {
     py::module m("pyigl", R"pyigldoc(
         Python wrappers for libigl
@@ -56,6 +60,7 @@ PYBIND11_PLUGIN(pyigl) {
            boundary_facets
            boundary_loop
            cat
+           collapse_edge
            colon
            comb_cross_field
            comb_frame_field
@@ -63,6 +68,8 @@ PYBIND11_PLUGIN(pyigl) {
            copyleft_cgal_mesh_boolean
            copyleft_comiso_miq
            copyleft_comiso_nrosy
+           copyleft_marching_cubes
+           copyleft_swept_volume
            copyleft_tetgen_tetrahedralize
            cotmatrix
            covariance_scatter_matrix
@@ -70,15 +77,21 @@ PYBIND11_PLUGIN(pyigl) {
            cut_mesh_from_singularities
            doublearea
            edge_lengths
+           edge_topology
            eigs
            embree_ambient_occlusion
+           embree_reorient_facets_raycast
            find_cross_field_singularities
            fit_rotations
            floor
            gaussian_curvature
+           get_seconds
            grad
            harmonic
+           hsv_to_rgb
+           internal_angles
            invert_diag
+           is_irregular_vertex
            jet
            local_basis
            lscm
@@ -98,6 +111,7 @@ PYBIND11_PLUGIN(pyigl) {
            polar_svd
            principal_curvature
            quad_planarity
+           randperm
            readDMAT
            readMESH
            readOBJ
@@ -111,10 +125,14 @@ PYBIND11_PLUGIN(pyigl) {
            slice_mask
            slice_tets
            sortrows
+           streamlines_init
+           streamlines_next
+           triangle_triangle_adjacency
            triangle_triangulate
            unique
            unproject_onto_mesh
            upsample
+           winding_number
            writeMESH
            writeOBJ
 
@@ -152,5 +170,9 @@ PYBIND11_PLUGIN(pyigl) {
     python_export_igl_png(m);
     #endif
 
+    #ifdef PY_COPYLEFT
+    python_export_igl_copyleft(m);
+    #endif
+
     return m.ptr();
 }

python/scripts/generate_bindings.py

@@ -79,7 +79,7 @@ def map_parameter_types(name, cpp_type, parsed_types, errors, enum_types):
 
     if len(parsed_types) == 0:
         errors.append("Empty typechain: %s" % cpp_type)
-        if cpp_type == "int" or cpp_type == "bool":
+        if cpp_type == "int" or cpp_type == "bool" or cpp_type == "unsigned int":
             return cpp_type, True
         else:
             return cpp_type, False

python/scripts/python_shared.mako

@@ -34,6 +34,10 @@ extern void python_export_igl_cgal(py::module &);
 extern void python_export_igl_png(py::module &);
 #endif
 
+#ifdef PY_COPYLEFT
+extern void python_export_igl_copyleft(py::module &);
+#endif
+
 PYBIND11_PLUGIN(pyigl) {
     py::module m("pyigl", R"pyigldoc(
         Python wrappers for libigl
@@ -82,5 +86,9 @@ PYBIND11_PLUGIN(pyigl) {
     python_export_igl_png(m);
     #endif
 
+    #ifdef PY_COPYLEFT
+    python_export_igl_copyleft(m);
+    #endif
+
     return m.ptr();
 }

python/tutorial/106_ViewerMenu.py

@@ -0,0 +1,83 @@
+# Add the igl library to the modules search path
+import sys, os
+sys.path.insert(0, os.getcwd() + "/../")
+
+import pyigl as igl
+import nanogui
+
+V1 = igl.eigen.MatrixXd()
+F1 = igl.eigen.MatrixXi()
+
+V2 = igl.eigen.MatrixXd()
+F2 = igl.eigen.MatrixXi()
+
+float_variable = 0.1
+bool_variable = True
+dir = 0
+
+
+def make_accessors(name):
+    def setter(value):
+        globals()[name] = value
+
+    def getter():
+        return globals()[name]
+    return setter, getter
+
+
+def viewer_init(viewer):
+    # add new group
+    viewer.ngui.addGroup("New Group")
+
+    # Expose the using general callback
+    viewer.ngui.addDoubleVariable("double", *make_accessors("float_variable"))
+
+    def setter(val):
+        global bool_variable
+        bool_variable = val
+
+    def getter():
+        global bool_variable
+        return bool_variable
+
+    # ... or using a custom callback
+    viewer.ngui.addBoolVariable("bool", setter, getter)
+
+    viewer.ngui.addEnumVariable("Direction", *make_accessors("dir")) \
+        .setItems(["Up", "Down", "Left", "Right"])
+
+    # Add a button
+    def cb():
+        print("Hello")
+    viewer.ngui.addButton("Print Hello", cb)
+
+    #Add an additional menu window
+    viewer.ngui.addWindow(nanogui.Vector2i(220, 10), "New Window")
+
+    # add accessor
+    viewer.ngui.addDoubleVariable("double", *make_accessors("float_variable"))
+
+    #Generate menu
+    viewer.screen.performLayout()
+
+    return False
+
+
+
+def main():
+    # Load a mesh in OFF format
+    igl.readOFF("../../tutorial/shared/bunny.off", V1, F1)
+
+    # Init the viewer
+    viewer = igl.viewer.Viewer()
+
+    # Extend viewer menu
+    viewer.callback_init = viewer_init
+
+    # Plot the mesh
+    viewer.data.set_mesh(V1, F1)
+    viewer.launch()
+
+
+if __name__ == "__main__":
+    main()

python/tutorial/701_Statistics.py

@@ -0,0 +1,56 @@
+import sys, os
+
+# Add the igl library to the modules search path
+import math
+
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+
+from shared import TUTORIAL_SHARED_PATH, check_dependencies
+
+dependencies = []
+check_dependencies(dependencies)
+
+if __name__ == "__main__":
+    V = igl.eigen.MatrixXd()
+    F = igl.eigen.MatrixXi()
+
+    # Load meshes in OFF format
+    igl.readOBJ(TUTORIAL_SHARED_PATH + "horse_quad.obj", V, F)
+
+    # Count the number of irregular vertices, the border is ignored
+    irregular = igl.is_irregular_vertex(V, F)
+    vertex_count = V.rows()
+    irregular_vertex_count = sum(irregular)
+    irregular_ratio = irregular_vertex_count / vertex_count
+
+    print("Irregular vertices: \n%d/%d (%.2f%%)\n" % (
+        irregular_vertex_count, vertex_count, irregular_ratio * 100))
+
+    # Compute areas, min, max and standard deviation
+    area = igl.eigen.MatrixXd()
+    igl.doublearea(V, F, area)
+    area /= 2.0
+
+    area_avg = area.mean()
+    area_min = area.minCoeff() / area_avg
+    area_max = area.maxCoeff() / area_avg
+    area_ns = (area - area_avg) / area_avg
+    area_sigma = math.sqrt(area_ns.squaredMean())
+
+    print("Areas (Min/Max)/Avg_Area Sigma: \n%.2f/%.2f (%.2f)\n" % (
+        area_min, area_max, area_sigma))
+
+    # Compute per face angles, min, max and standard deviation
+    angles = igl.eigen.MatrixXd()
+    igl.internal_angles(V, F, angles)
+    angles = 360.0 * (angles / (2 * math.pi))
+
+    angle_avg = angles.mean()
+    angle_min = angles.minCoeff()
+    angle_max = angles.maxCoeff()
+    angle_ns = angles - angle_avg
+    angle_sigma = math.sqrt(angle_ns.squaredMean())
+
+    print("Angles in degrees (Min/Max) Sigma: \n%.2f/%.2f (%.2f)\n" % (
+        angle_min, angle_max, angle_sigma))

python/tutorial/702_WindingNumber.py

@@ -0,0 +1,112 @@
+import sys, os
+
+# Add the igl library to the modules search path
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+
+from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
+
+dependencies = ["viewer"]
+check_dependencies(dependencies)
+
+
+def append_mesh(C_vis, F_vis, V_vis, V, F, color):
+    F_vis.conservativeResize(F_vis.rows() + F.rows(), 3)
+    F_vis.setBottomRows(F.rows(), F + V_vis.rows())
+    V_vis.conservativeResize(V_vis.rows() + V.rows(), 3)
+    V_vis.setBottomRows(V.rows(), V)
+    C_vis.conservativeResize(C_vis.rows() + F.rows(), 3)
+    colorM = igl.eigen.MatrixXd(F.rows(), C_vis.cols())
+    colorM.rowwiseSet(color)
+    C_vis.setBottomRows(F.rows(), colorM)
+
+
+def update(viewer):
+    global V, F, T, W, slice_z, overlay
+    plane = igl.eigen.MatrixXd([0, 0, 1, -((1 - slice_z) * V.col(2).minCoeff() + slice_z * V.col(2).maxCoeff())])
+    V_vis = igl.eigen.MatrixXd()
+    F_vis = igl.eigen.MatrixXi()
+    J = igl.eigen.MatrixXi()
+    bary = igl.eigen.SparseMatrixd()
+    igl.slice_tets(V, T, plane, V_vis, F_vis, J, bary)
+    W_vis = igl.eigen.MatrixXd()
+    igl.slice(W, J, W_vis)
+    C_vis = igl.eigen.MatrixXd()
+    igl.parula(W_vis, False, C_vis)
+
+    if overlay == 1:  # OVERLAY_INPUT
+        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[1., 0.894, 0.227]]))
+    elif overlay == 2:  # OVERLAY_OUTPUT
+        append_mesh(C_vis, F_vis, V_vis, V, F, igl.eigen.MatrixXd([[0.8, 0.8, 0.8]]))
+
+    viewer.data.clear()
+    viewer.data.set_mesh(V_vis, F_vis)
+    viewer.data.set_colors(C_vis)
+    viewer.data.set_face_based(True)
+
+
+def key_down(viewer, key, modifier):
+    global overlay, slice_z
+
+    if key == ord(' '):
+        overlay = (overlay + 1) % 3
+    elif key == ord('.'):
+        slice_z = min(slice_z + 0.01, 0.99)
+    elif key == ord(','):
+        slice_z = max(slice_z - 0.01, 0.01)
+
+    update(viewer)
+
+    return False
+
+
+if __name__ == "__main__":
+    keys = {"space": "toggle showing input mesh, output mesh or slice through tet-mesh of convex hull",
+            ". / ,": "push back/pull forward slicing plane"}
+
+    print_usage(keys)
+
+    V = igl.eigen.MatrixXd()
+    BC = igl.eigen.MatrixXd()
+    W = igl.eigen.MatrixXd()
+    T = igl.eigen.MatrixXi()
+    F = igl.eigen.MatrixXi()
+    G = igl.eigen.MatrixXi()
+
+    slice_z = 0.5
+    overlay = 0
+
+    # Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
+    # surface mesh _after_ self-intersection resolution
+    igl.readMESH(TUTORIAL_SHARED_PATH + "big-sigcat.mesh", V, T, F)
+
+    # Compute barycenters of all tets
+    igl.barycenter(V, T, BC)
+
+    # Compute generalized winding number at all barycenters
+    print("Computing winding number over all %i tets..." % T.rows())
+    igl.winding_number(V, F, BC, W)
+
+    # Extract interior tets
+    Wt = sum(W > 0.5)
+    CT = igl.eigen.MatrixXi(Wt, 4)
+    k = 0
+    for t in range(T.rows()):
+        if W[t] > 0.5:
+            CT.setRow(k, T.row(t))
+            k += 1
+
+    # find bounary facets of interior tets
+    igl.boundary_facets(CT, G)
+
+    # boundary_facets seem to be reversed...
+    G = G.rowwiseReverse()
+
+    # normalize
+    W = (W - W.minCoeff()) / (W.maxCoeff() - W.minCoeff())
+
+    # Plot the generated mesh
+    viewer = igl.viewer.Viewer()
+    update(viewer)
+    viewer.callback_key_down = key_down
+    viewer.launch()

python/tutorial/705_MarchingCubes.py

@@ -0,0 +1,94 @@
+import sys, os
+
+# Add the igl library to the modules search path
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+
+from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
+
+dependencies = ["copyleft", "viewer"]
+check_dependencies(dependencies)
+
+
+def key_down(viewer, key, modifier):
+    if key == ord('1'):
+        viewer.data.clear()
+        viewer.data.set_mesh(V, F)
+    elif key == ord('2'):
+        viewer.data.clear()
+        viewer.data.set_mesh(SV, SF)
+    elif key == ord('3'):
+        viewer.data.clear()
+        viewer.data.set_mesh(BV, BF)
+
+    return True
+
+
+if __name__ == "__main__":
+    keys = {"1": "show original mesh",
+            "2": "show marching cubes contour of signed distance",
+            "3": "show marching cubes contour of indicator function"}
+
+    print_usage(keys)
+
+    V = igl.eigen.MatrixXd()
+    F = igl.eigen.MatrixXi()
+
+    # Read in inputs as double precision floating point meshes
+    igl.read_triangle_mesh(TUTORIAL_SHARED_PATH + "armadillo.obj", V, F)
+
+    # number of vertices on the largest side
+    s = 50
+    Vmin = V.colwiseMinCoeff()
+    Vmax = V.colwiseMaxCoeff()
+    h = (Vmax - Vmin).maxCoeff() / s
+    res = (s * ((Vmax - Vmin) / (Vmax - Vmin).maxCoeff())).castint()
+
+    def lerp(res, Vmin, Vmax, di, d):
+        return Vmin[d] + float(di) / (res[d] - 1) * (Vmax[d] - Vmin[d])
+
+    # create grid
+    print("Creating grid...")
+    GV = igl.eigen.MatrixXd(res[0] * res[1] * res[2], 3)
+    for zi in range(res[2]):
+        z = lerp(res, Vmin, Vmax, zi, 2)
+        for yi in range(res[1]):
+            y = lerp(res, Vmin, Vmax, yi, 1)
+            for xi in range(res[0]):
+                x = lerp(res, Vmin, Vmax, xi, 0)
+                GV.setRow(xi + res[0] * (yi + res[1] * zi), igl.eigen.MatrixXd([[x, y, z]]))
+
+    # compute values
+    print("Computing distances...")
+    S = igl.eigen.MatrixXd()
+    B = igl.eigen.MatrixXd()
+    I = igl.eigen.MatrixXi()
+    C = igl.eigen.MatrixXd()
+    N = igl.eigen.MatrixXd()
+
+    igl.signed_distance(GV, V, F, igl.SIGNED_DISTANCE_TYPE_PSEUDONORMAL, S, I, C, N)
+    # Convert distances to binary inside-outside data --> aliasing artifacts
+    B = S.copy()
+    for e in range(B.rows()):
+        if B[e] > 0:
+            B[e] = 1
+        else:
+            if B[e] < 0:
+                B[e] = -1
+            else:
+                B[e] = 0
+
+    print("Marching cubes...")
+    SV = igl.eigen.MatrixXd()
+    BV = igl.eigen.MatrixXd()
+    SF = igl.eigen.MatrixXi()
+    BF = igl.eigen.MatrixXi()
+
+    igl.copyleft.marching_cubes(S, GV, res[0], res[1], res[2], SV, SF)
+    igl.copyleft.marching_cubes(B, GV, res[0], res[1], res[2], BV, BF)
+
+    # Plot the generated mesh
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(SV, SF)
+    viewer.callback_key_down = key_down
+    viewer.launch()

python/tutorial/706_FacetOrientation.py

@@ -0,0 +1,78 @@
+import sys, os
+
+# Add the igl library to the modules search path
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+
+from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
+
+dependencies = ["embree", "viewer"]
+check_dependencies(dependencies)
+
+
+def key_down(viewer, key, modifier):
+    global facetwise, is_showing_reoriented, FF
+    if key == ord('F') or key == ord('f'):
+        facetwise = (facetwise + 1) % 2
+    elif key == ord('S') or key == ord('s'):
+        scramble_colors()
+    elif key == ord(' '):
+        is_showing_reoriented = ~is_showing_reoriented
+
+    viewer.data.clear()
+    viewer.data.set_mesh(V, FF[facetwise] if is_showing_reoriented else F)
+    viewer.data.set_colors(RGBcolors[facetwise])
+
+    return True
+
+def scramble_colors():
+    global C, viewer, RGBcolors
+    for p in range(2):
+        R = igl.eigen.MatrixXi()
+        igl.randperm(C[p].maxCoeff() + 1, R)
+        C[p] = igl.slice(R, igl.eigen.MatrixXi(C[p]))
+        HSV = igl.eigen.MatrixXd(C[p].rows(), 3)
+        HSV.setCol(0, 360.0 * C[p].castdouble() / C[p].maxCoeff())
+        HSVright = igl.eigen.MatrixXd(HSV.rows(), 2)
+        HSVright.setConstant(1.0)
+        HSV.setRightCols(2, HSVright)
+        igl.hsv_to_rgb(HSV, RGBcolors[p])
+    viewer.data.set_colors(RGBcolors[facetwise])
+
+
+
+if __name__ == "__main__":
+    keys = {"space": "toggle between original and reoriented faces",
+            "F,f": "toggle between patchwise and facetwise reorientation",
+            "S,s": "scramble colors"}
+    print_usage(keys)
+
+    V = igl.eigen.MatrixXd()
+    F = igl.eigen.MatrixXi()
+    C = [igl.eigen.MatrixXi(), igl.eigen.MatrixXi()]
+    RGBcolors = [igl.eigen.MatrixXd(), igl.eigen.MatrixXd()]
+    FF = [igl.eigen.MatrixXi(), igl.eigen.MatrixXi()]
+    is_showing_reoriented = False
+    facetwise = 0
+
+    igl.read_triangle_mesh(TUTORIAL_SHARED_PATH + "truck.obj", V, F)
+
+    # Compute patches
+    for p in range(2):
+        I = igl.eigen.MatrixXi()
+        igl.embree.reorient_facets_raycast(V, F, F.rows() * 100, 10, p == 1, False, False, I, C[p])
+        # apply reorientation
+        FF[p].conservativeResize(F.rows(), F.cols())
+        for i in range(I.rows()):
+            if I[i]:
+                FF[p].setRow(i, F.row(i).rowwiseReverse())
+            else:
+                FF[p].setRow(i, F.row(i))
+
+    # Plot the generated mesh
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(V, FF[facetwise] if is_showing_reoriented else F)
+    viewer.data.set_face_based(True)
+    scramble_colors()
+    viewer.callback_key_down = key_down
+    viewer.launch()

python/tutorial/707_SweptVolume.py

@@ -0,0 +1,83 @@
+import sys, os
+
+# Add the igl library to the modules search path
+from math import pi, cos
+
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+
+from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
+
+dependencies = ["copyleft", "viewer"]
+check_dependencies(dependencies)
+
+
+def key_down(viewer, key, modifier):
+    global show_swept_volume, SV, SF, V, F
+    if key == ord(' '):
+        show_swept_volume = not show_swept_volume
+        viewer.data.clear()
+
+        if show_swept_volume:
+            viewer.data.set_mesh(SV, SF)
+            viewer.data.uniform_colors(igl.eigen.MatrixXd([0.2, 0.2, 0.2]), igl.eigen.MatrixXd([1.0, 1.0, 1.0]), igl.eigen.MatrixXd([1.0, 1.0, 1.0])) # TODO replace with constants from cpp
+        else:
+            viewer.data.set_mesh(V, F)
+
+        viewer.core.is_animating = not show_swept_volume
+        viewer.data.set_face_based(True)
+
+    return True
+
+
+def pre_draw(viewer):
+    global show_swept_volume, V
+    if not show_swept_volume:
+        T = transform(0.25 * igl.get_seconds())
+        VT = V * T.matrix().block(0, 0, 3, 3).transpose()
+        trans = T.matrix().block(0, 3, 3, 1).transpose()
+        Vtrans = igl.eigen.MatrixXd(VT.rows(), VT.cols())
+        Vtrans.rowwiseSet(trans)
+        VT += Vtrans
+        viewer.data.set_vertices(VT)
+        viewer.data.compute_normals()
+    return False
+
+
+# Define a rigid motion
+def transform(t):
+    T = igl.eigen.Affine3d()
+    T.setIdentity()
+    T.rotate(t * 2 * pi, igl.eigen.MatrixXd([0, 1, 0]))
+    T.translate(igl.eigen.MatrixXd([0, 0.125 * cos(2 * pi * t), 0]))
+    return T
+
+
+if __name__ == "__main__":
+    keys = {"space": "toggle between transforming original mesh and swept volume"}
+    print_usage(keys)
+
+    V = igl.eigen.MatrixXd()
+    SV = igl.eigen.MatrixXd()
+    VT = igl.eigen.MatrixXd()
+    F = igl.eigen.MatrixXi()
+    SF = igl.eigen.MatrixXi()
+    show_swept_volume = False
+    grid_size = 50
+    time_steps = 200
+    isolevel = 1
+
+    igl.read_triangle_mesh(TUTORIAL_SHARED_PATH + "bunny.off", V, F)
+
+    print("Computing swept volume...")
+    igl.copyleft.swept_volume(V, F, transform, time_steps, grid_size, isolevel, SV, SF)
+    print("...finished.")
+
+    # Plot the generated mesh
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(V, F)
+    viewer.data.set_face_based(True)
+    viewer.core.is_animating = not show_swept_volume
+    viewer.callback_pre_draw = pre_draw
+    viewer.callback_key_down = key_down
+    viewer.launch()

python/tutorial/607_Picking.py → python/tutorial/708_Picking.py

@@ -5,19 +5,12 @@ sys.path.insert(0, os.getcwd() + "/../")
 import pyigl as igl
 
 
-from shared import TUTORIAL_SHARED_PATH, check_dependencies
+from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
 
 dependencies = ["viewer"]
 check_dependencies(dependencies)
 
 
-# Mesh with per-face color
-V = igl.eigen.MatrixXd()
-F = igl.eigen.MatrixXi()
-C = igl.eigen.MatrixXd()
-
-viewer = igl.viewer.Viewer()
-
 def mouse_down(viewer, a, b):
     bc = igl.eigen.MatrixXd()
 
@@ -27,6 +20,7 @@ def mouse_down(viewer, a, b):
     hit = igl.unproject_onto_mesh(coord, viewer.core.view * viewer.core.model,
       viewer.core.proj, viewer.core.viewport, V, F, fid, bc)
     if hit:
+        # paint hit red
         C.setRow(fid[0, 0], igl.eigen.MatrixXd([[1, 0, 0]]))
         viewer.data.set_colors(C)
         return True
@@ -34,16 +28,25 @@ def mouse_down(viewer, a, b):
     return False
 
 
-print("Usage: [LeftMouseClick] to select a face")
+if __name__ == "__main__":
+    keys = {"click": "Pick face on shape"}
+    print_usage(keys)
+
+    # Mesh with per-face color
+    V = igl.eigen.MatrixXd()
+    F = igl.eigen.MatrixXi()
+    C = igl.eigen.MatrixXd()
 
-# Load a mesh in OFF format
-igl.readOFF(TUTORIAL_SHARED_PATH + "fertility.off", V, F)
+    # Load a mesh in OFF format
+    igl.readOFF(TUTORIAL_SHARED_PATH + "fertility.off", V, F)
 
-# Initialize white
-C.setConstant(F.rows(), 3, 1.0)
+    # Initialize white
+    C.setConstant(F.rows(), 3, 1.0)
 
-viewer.data.set_mesh(V, F)
-viewer.data.set_colors(C)
-viewer.core.show_lines = False
-viewer.callback_mouse_down = mouse_down
-viewer.launch()
+    # Show mesh
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(V, F)
+    viewer.data.set_colors(C)
+    viewer.core.show_lines = False
+    viewer.callback_mouse_down = mouse_down
+    viewer.launch()

python/tutorial/709_VectorFieldVisualizer.py

@@ -0,0 +1,126 @@
+import sys, os
+import numpy as np
+
+# Add the igl library to the modules search path
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+from shared import TUTORIAL_SHARED_PATH, check_dependencies
+
+dependencies = ["viewer"]
+check_dependencies(dependencies)
+
+# Input mesh
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+
+data = igl.StreamlineData()
+state = igl.StreamlineState()
+
+treat_as_symmetric = True
+
+# animation params
+anim_t = 0
+anim_t_dir = 1
+
+
+def representative_to_nrosy(V, F, R, N, Y):
+    B1 = igl.eigen.MatrixXd()
+    B2 = igl.eigen.MatrixXd()
+    B3 = igl.eigen.MatrixXd()
+
+    igl.local_basis(V, F, B1, B2, B3)
+
+    Y.resize(F.rows(), 3 * N)
+    for i in range(0, F.rows()):
+        x = R.row(i) * B1.row(i).transpose()
+        y = R.row(i) * B2.row(i).transpose()
+        angle = np.arctan2(y, x)
+
+        for j in range(0, N):
+            anglej = angle + np.pi * float(j) / float(N)
+            xj = float(np.cos(anglej))
+            yj = float(np.sin(anglej))
+            Y.setBlock(i, j * 3, 1, 3, xj * B1.row(i) + yj * B2.row(i))
+
+
+def pre_draw(viewer):
+    if not viewer.core.is_animating:
+        return False
+
+    global anim_t
+    global start_point
+    global end_point
+
+    igl.streamlines_next(V, F, data, state)
+
+    value = (anim_t % 100) / 100.0
+
+    if value > 0.5:
+        value = 1 - value
+    value /= 0.5
+    r, g, b = igl.parula(value)
+    viewer.data.add_edges(state.start_point, state.end_point, igl.eigen.MatrixXd([[r, g, b]]))
+
+    anim_t += anim_t_dir
+
+    return False
+
+
+def key_down(viewer, key, modifier):
+    if key == ord(' '):
+        viewer.core.is_animating = not viewer.core.is_animating
+        return True
+
+    return False
+
+
+def main():
+    # Load a mesh in OFF format
+    igl.readOFF(TUTORIAL_SHARED_PATH + "bumpy.off", V, F)
+
+    # Create a Vector Field
+    temp_field = igl.eigen.MatrixXd()
+    b = igl.eigen.MatrixXi([[0]])
+    bc = igl.eigen.MatrixXd([[1, 1, 1]])
+    S = igl.eigen.MatrixXd()  # unused
+
+    degree = 3
+    igl.comiso.nrosy(V, F, b, bc, igl.eigen.MatrixXi(), igl.eigen.MatrixXd(), igl.eigen.MatrixXd(), 1, 0.5, temp_field, S)
+    temp_field2 = igl.eigen.MatrixXd()
+    representative_to_nrosy(V, F, temp_field, degree, temp_field2)
+
+    # Initialize tracer
+    igl.streamlines_init(V, F, temp_field2, treat_as_symmetric, data, state)
+
+    # Setup viewer
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(V, F)
+    viewer.callback_pre_draw = pre_draw
+    viewer.callback_key_down = key_down
+
+    viewer.core.show_lines = False
+
+    viewer.core.is_animating = False
+    viewer.core.animation_max_fps = 30.0
+
+    # Paint mesh grayish
+    C = igl.eigen.MatrixXd()
+    C.setConstant(viewer.data.V.rows(), 3, .9)
+    viewer.data.set_colors(C)
+
+    # Draw vector field on sample points
+    state0 = state.copy()
+
+    igl.streamlines_next(V, F, data, state0)
+    v = state0.end_point - state0.start_point
+    v = v.rowwiseNormalized()
+
+    viewer.data.add_edges(state0.start_point,
+                          state0.start_point + 0.059 * v,
+                          igl.eigen.MatrixXd([[1.0, 1.0, 1.0]]))
+
+    print("Press [space] to toggle animation")
+    viewer.launch()
+
+if __name__ == "__main__":
+    main()

python/tutorial/shared.py

@@ -14,3 +14,9 @@ def check_dependencies(deps):
 
     if not all_available:
         sys.exit(-1)
+
+
+def print_usage(key_dict):
+    print("Usage:")
+    for k in key_dict.keys():
+        print("%s : %s" %(k, key_dict[k]))

shared/cmake/CMakeLists.txt

@@ -19,7 +19,7 @@ option(LIBIGL_WITH_TETGEN           "Use Tetgen"         OFF)
 option(LIBIGL_WITH_TRIANGLE         "Use Triangle"       OFF)
 option(LIBIGL_WITH_VIEWER           "Use OpenGL viewer"  OFF)
 option(LIBIGL_WITH_XML              "Use XML"            OFF)
-
+option(LIBIGL_WITH_PYTHON           "Use Python"         OFF)
 
 ### Compilation configuration ###
 if(MSVC)
@@ -102,14 +102,14 @@ endif()
 if(LIBIGL_WITH_ANTTWEAKBAR)
   set(ANTTWEAKBAR_DIR "${LIBIGL_EXTERNAL}/anttweakbar")
   set(ANTTWEAKBAR_INCLUDE_DIR "${ANTTWEAKBAR_DIR}/include")
-  set(ANTTWEAKBAR_C_SRC_FILES 
+  set(ANTTWEAKBAR_C_SRC_FILES
     "${ANTTWEAKBAR_DIR}/src/TwEventGLFW.c"
     "${ANTTWEAKBAR_DIR}/src/TwEventGLUT.c"
     "${ANTTWEAKBAR_DIR}/src/TwEventSDL.c"
     "${ANTTWEAKBAR_DIR}/src/TwEventSDL12.c"
     "${ANTTWEAKBAR_DIR}/src/TwEventSDL13.c"
     )
-  set(ANTTWEAKBAR_CPP_SRC_FILES 
+  set(ANTTWEAKBAR_CPP_SRC_FILES
     "${ANTTWEAKBAR_DIR}/src/LoadOGL.cpp"
     "${ANTTWEAKBAR_DIR}/src/LoadOGLCore.cpp"
     "${ANTTWEAKBAR_DIR}/src/TwBar.cpp"
@@ -129,15 +129,15 @@ if(LIBIGL_WITH_ANTTWEAKBAR)
     #"${ANTTWEAKBAR_DIR}/src/TwDirect3D11.cpp"
     #"${ANTTWEAKBAR_DIR}/src/TwDirect3D9.cpp"
   list(
-    APPEND 
-    ANTTWEAKBAR_SRC_FILES 
+    APPEND
+    ANTTWEAKBAR_SRC_FILES
     "${ANTTWEAKBAR_C_SRC_FILES}"
     "${ANTTWEAKBAR_CPP_SRC_FILES}")
   add_library(AntTweakBar STATIC "${ANTTWEAKBAR_SRC_FILES}")
   target_include_directories(AntTweakBar PUBLIC "${ANTTWEAKBAR_INCLUDE_DIR}")
   if(APPLE)
     set_target_properties(
-      AntTweakBar 
+      AntTweakBar
       PROPERTIES
       COMPILE_FLAGS
       "-fPIC -fno-strict-aliasing -x objective-c++")
@@ -286,7 +286,7 @@ if(LIBIGL_WITH_LIM)
   add_subdirectory("${LIM_DIR}" "lim")
   list(APPEND LIBIGL_INCLUDE_DIRS ${LIM_DIR})
   ## it depends on ligigl, so placing it here solve linking problems
-  #list(APPEND LIBIGL_LIBRARIES "lim") 
+  #list(APPEND LIBIGL_LIBRARIES "lim")
   # ^--- Alec: I don't understand this comment. Does lim need to come before
   # libigl libraries? Why can't lim be placed where it belongs in
   # LIBIGL_EXTRA_LIBRARIES?
@@ -340,6 +340,9 @@ endif()
 
 ### Compile the png parts ###
 if(LIBIGL_WITH_PNG)
+  if(LIBIGL_WITH_NANOGUI)
+    set(STBI_LOAD OFF CACHE BOOL " " FORCE)
+  endif()
   set(STB_IMAGE_DIR "${LIBIGL_EXTERNAL}/stb_image")
   add_subdirectory("${STB_IMAGE_DIR}" "stb_image")
   list(APPEND LIBIGL_INCLUDE_DIRS ${STB_IMAGE_DIR})
@@ -389,7 +392,12 @@ if(LIBIGL_WITH_VIEWER)
   if(LIBIGL_WITH_NANOGUI)
     list(APPEND LIBIGL_DEFINITIONS "-DIGL_VIEWER_WITH_NANOGUI")
 
-    set(NANOGUI_BUILD_PYTHON OFF CACHE BOOL " " FORCE)
+    if (LIBIGL_WITH_PYTHON)
+      set(NANOGUI_BUILD_PYTHON ON CACHE BOOL " " FORCE)
+    else()
+      set(NANOGUI_BUILD_PYTHON OFF CACHE BOOL " " FORCE)
+    endif()
+
     set(NANOGUI_BUILD_EXAMPLE OFF CACHE BOOL " " FORCE)
     set(NANOGUI_BUILD_SHARED  OFF CACHE BOOL " " FORCE)
     add_subdirectory("${NANOGUI_DIR}" "nanogui")

style-guidelines.html

@@ -185,8 +185,8 @@ Eigen::SparseMatrix<Atype> adjacency_matrix(const ... & F);
 <h2 id="templatingwitheigen">Templating with Eigen</h2>
 
 <p>Functions taking Eigen dense matrices/arrays as inputs and outputs (but <strong>not</strong>
-return arguments), should template on top of <code>Eigen::PlainObjectBase</code>. **Each
-parameter** should be derived using its own template.</p>
+return arguments), should template on top of <code>Eigen::PlainObjectBase</code>. <strong>Each
+parameter</strong> should be derived using its own template.</p>
 
 <p>For example,</p>
 
@@ -264,8 +264,8 @@ unnecessary prefaces. For example, instead of <code>compute_adjacency_matrix</co
 
 <h2 id="variablenamingconventions">Variable naming conventions</h2>
 
-<p>Libigl prefers short (even single character) variable names _with heavy
-documentation_ in the comments in the header file or above the declaration of
+<p>Libigl prefers short (even single character) variable names <em>with heavy
+documentation</em> in the comments in the header file or above the declaration of
 the function. When possible use <code>V</code> to mean a list of vertex positions and <code>F</code>
 to mean a list of faces/triangles.</p>
 

tutorial/401_BiharmonicDeformation/main.cpp

@@ -34,7 +34,7 @@ bool pre_draw(igl::viewer::Viewer & viewer)
     U = V+D;
   }else
   {
-    igl::harmonic(V,F,b,U_bc_anim,2,U);
+    igl::harmonic(V,F,b,U_bc_anim,2.,U);
   }
   viewer.data.set_vertices(U);
   viewer.data.compute_normals();
@@ -66,7 +66,7 @@ int main(int argc, char *argv[])
   VectorXi S;
   igl::readDMAT(TUTORIAL_SHARED_PATH "/decimated-max-selection.dmat",S);
   igl::colon<int>(0,V.rows()-1,b);
-  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
+  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(),
    [&S](int i)->bool{return S(i)>=0;})-b.data());
 
   // Boundary conditions directly on deformed positions

tutorial/402_PolyharmonicDeformation/main.cpp

@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
   VectorXb is_inner = (V.rowwise().norm().array()-0.15)<1e-15;
   VectorXb in_b = is_outer.array() || is_inner.array();
   igl::colon<int>(0,V.rows()-1,b);
-  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(), 
+  b.conservativeResize(stable_partition( b.data(), b.data()+b.size(),
    [&in_b](int i)->bool{return in_b(i);})-b.data());
   bc.resize(b.size(),1);
   for(int bi = 0;bi<b.size();bi++)

tutorial/403_BoundedBiharmonicWeights/main.cpp

@@ -32,7 +32,7 @@
 
 #include "tutorial_shared_path.h"
 
-typedef 
+typedef
   std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
   RotationList;
 
@@ -78,7 +78,7 @@ bool pre_draw(igl::viewer::Viewer & viewer)
     MatrixXd CT;
     MatrixXi BET;
     igl::deform_skeleton(C,BE,T,CT,BET);
-    
+
     viewer.data.set_vertices(U);
     viewer.data.set_edges(CT,BET,sea_green);
     viewer.data.compute_normals();

tutorial/501_HarmonicParam/main.cpp

@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
 
   // Harmonic parametrization for the internal vertices
   igl::harmonic(V,F,bnd,bnd_uv,1,V_uv);
-  
+
   // Scale UV to make the texture more clear
   V_uv *= 5;
 

tutorial/709_VectorFieldVisualizer/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(709_VectorFieldVisualizer)
+
+add_executable(${PROJECT_NAME}_bin
+        main.cpp)
+target_include_directories(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_INCLUDE_DIRS})
+target_compile_definitions(${PROJECT_NAME}_bin PRIVATE ${LIBIGL_DEFINITIONS})
+target_link_libraries(${PROJECT_NAME}_bin ${LIBIGL_LIBRARIES} ${LIBIGL_EXTRA_LIBRARIES})

tutorial/709_VectorFieldVisualizer/main.cpp

@@ -0,0 +1,162 @@
+#include <igl/barycenter.h>
+#include <igl/edge_topology.h>
+#include <igl/local_basis.h>
+#include <igl/parula.h>
+#include <igl/per_face_normals.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/polyvector_field_matchings.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/readOFF.h>
+#include <igl/slice.h>
+#include <igl/sort_vectors_ccw.h>
+#include <igl/streamlines.h>
+#include <igl/triangle_triangle_adjacency.h>
+#include <igl/copyleft/comiso/nrosy.h>
+#include <igl/viewer/Viewer.h>
+
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+#include <fstream>
+
+
+// Mesh
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+igl::StreamlineData sl_data;
+igl::StreamlineState sl_state;
+
+int degree;         // degree of the vector field
+int half_degree;    // degree/2 if treat_as_symmetric
+bool treat_as_symmetric = true;
+
+int anim_t = 0;
+int anim_t_dir = 1;
+
+
+void representative_to_nrosy(
+        const Eigen::MatrixXd &V,
+        const Eigen::MatrixXi &F,
+        const Eigen::MatrixXd &R,
+        const int N,
+        Eigen::MatrixXd &Y)
+{
+    using namespace Eigen;
+    using namespace std;
+    MatrixXd B1, B2, B3;
+
+    igl::local_basis(V, F, B1, B2, B3);
+
+    Y.resize(F.rows(), 3 * N);
+    for (unsigned i = 0; i < F.rows(); ++i)
+    {
+        double x = R.row(i) * B1.row(i).transpose();
+        double y = R.row(i) * B2.row(i).transpose();
+        double angle = atan2(y, x);
+
+        for (unsigned j = 0; j < N; ++j)
+        {
+            double anglej = angle + M_PI * double(j) / double(N);
+            double xj = cos(anglej);
+            double yj = sin(anglej);
+            Y.block(i, j * 3, 1, 3) = xj * B1.row(i) + yj * B2.row(i);
+        }
+    }
+}
+
+bool pre_draw(igl::viewer::Viewer &viewer)
+{
+    using namespace Eigen;
+    using namespace std;
+
+    if (!viewer.core.is_animating)
+        return false;
+
+    igl::streamlines_next(V, F, sl_data, sl_state);
+    Eigen::RowVector3d color = Eigen::RowVector3d::Zero();
+    double value = ((anim_t) % 100) / 100.;
+
+    if (value > 0.5)
+        value = 1 - value;
+    value = value / 0.5;
+    igl::parula(value, color[0], color[1], color[2]);
+
+    viewer.data.add_edges(sl_state.start_point, sl_state.end_point, color);
+
+    anim_t += anim_t_dir;
+
+    return false;
+}
+
+bool key_down(igl::viewer::Viewer &viewer, unsigned char key, int modifier)
+{
+    if (key == ' ')
+    {
+        viewer.core.is_animating = !viewer.core.is_animating;
+        return true;
+    }
+    return false;
+}
+
+int main(int argc, char *argv[])
+{
+    using namespace Eigen;
+    using namespace std;
+
+
+    // Load a mesh in OFF format
+    igl::readOFF(TUTORIAL_SHARED_PATH "/bumpy.off", V, F);
+    // Create a Vector Field
+    Eigen::VectorXi b;
+    Eigen::MatrixXd bc;
+    Eigen::VectorXd S; // unused
+
+    b.resize(1);
+    b << 0;
+    bc.resize(1, 3);
+    bc << 1, 1, 1;
+
+    half_degree = 3;
+    treat_as_symmetric = true;
+
+    Eigen::MatrixXd temp_field, temp_field2;
+    igl::copyleft::comiso::nrosy(V, F, b, bc, VectorXi(), VectorXd(), MatrixXd(), 1, 0.5, temp_field, S);
+    representative_to_nrosy(V, F, temp_field, half_degree, temp_field2);
+
+
+    igl::streamlines_init(V, F, temp_field2, treat_as_symmetric, sl_data, sl_state);
+
+
+    // Viewer Settings
+    igl::viewer::Viewer viewer;
+    viewer.data.set_mesh(V, F);
+    viewer.callback_pre_draw = &pre_draw;
+    viewer.callback_key_down = &key_down;
+
+    viewer.core.show_lines = false;
+
+    viewer.core.is_animating = false;
+    viewer.core.animation_max_fps = 30.;
+
+    // Paint mesh grayish
+    Eigen::MatrixXd C;
+    C.setConstant(viewer.data.V.rows(), 3, .9);
+    viewer.data.set_colors(C);
+
+
+    // Draw vector field on sample points
+    igl::StreamlineState sl_state0;
+    sl_state0 = sl_state;
+    igl::streamlines_next(V, F, sl_data, sl_state0);
+    Eigen::MatrixXd v = sl_state0.end_point - sl_state0.start_point;
+    v.rowwise().normalize();
+
+    viewer.data.add_edges(sl_state0.start_point,
+                          sl_state0.start_point + 0.059 * v,
+                          Eigen::RowVector3d::Constant(1.0f));
+
+    cout <<
+    "Press [space] to toggle animation" << endl;
+    viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -180,4 +180,6 @@ if(TUTORIALS_CHAPTER7)
   endif()
   add_subdirectory("707_SweptVolume")
   add_subdirectory("708_Picking")
+  add_subdirectory("709_VectorFieldVisualizer")
+
 endif()

tutorial/images/streamlines.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+8a716f5201965b5b5b9b8ce8cfee1aa0ef8e970f

tutorial/shared/README.md

@@ -0,0 +1,83 @@
+This directory contains a slew of models and data needed to run the tutorial
+examples. This README.md contains an attempt to track down original sources of
+this data for purposes of attribution. In some cases, only the oldest known
+use is listed.
+
+| Filename                           | Source             |
+|------------------------------------|--------------------|
+| 2triangles.off                     | [#libigl][]        |
+| 3holes.off                         | ?                  |
+| arm-weights.dmat                   | [#libigl][]        |
+| arm.obj                            | ?                  |
+| arm.tgf                            | [#libigl][]        |
+| armadillo-weights.dmat             | [#libigl][]        |
+| armadillo.obj                      | [#stanford][]      |
+| beetle.off                         | [#ivan][]          |
+| big-sigcat.mesh                    | [#jacobson_2013][] |
+| bump-domain.obj                    | ?                  |
+| bumpy-cube.dmat                    | ?                  |
+| bumpy-cube.obj                     | ?                  |
+| bumpy.off                          | ?                  |
+| bunny.mesh                         | [#stanford][]      |
+| bunny.off                          | [#stanford][]      |
+| camelhead.off                      | [#sorkine_2004][]  |
+| cheburashka-scalar.dmat            | [#libigl][]        |
+| cheburashka.off                    | [#cosmic_blobs][]  |
+| cow.off                            | ?                  |
+| cube.obj                           | [#libigl][]        |
+| cube.off                           | [#libigl][]        |
+| decimated-knight-selection.dmat    | [#libigl][]        |
+| decimated-knight.off               | [#cosmic_blobs][]  |
+| decimated-max-selection.dmat       | [#libigl][]        |
+| decimated-max.obj                  | [#mpi][]           |
+| fandisk.off                        | [#aim_at_shape][]  |
+| fertility.off                      | [#aim_at_shape][]  |
+| grid.off                           | [#libigl][]        |
+| hand-pose.dmat                     | [#libigl][]        |
+| hand.mesh                          | ?                  |
+| hand.tgf                           | [#libigl][]        |
+| horse_quad.obj                     | ?                  |
+| inspired_mesh.dmat                 | [#libigl][]        |
+| inspired_mesh.obj                  | ?                  |
+| inspired_mesh_b.dmat               | [#libigl][]        |
+| inspired_mesh_bc.dmat              | [#libigl][]        |
+| inspired_mesh_quads_Conjugate.off  | ?                  |
+| inspired_mesh_quads_Smooth.off     | ?                  |
+| lilium.crossfield                  | [#libigl][]        |
+| lilium.obj                         | ?                  |
+| lilium.samples.0.2                 | [#libigl][]        |
+| lilium_b.dmat                      | [#libigl][]        |
+| lilium_bc.dmat                     | [#libigl][]        |
+| lion.off                           | ?                  |
+| octopus-high.mesh                  | [#pauly][]         |
+| octopus-low.mesh                   | [#pauly][]         |
+| planexy.off                        | ?                  |
+| screwdriver.off                    | ?                  |
+| snail.dmat                         | [#libigl][]        |
+| snail.obj                          | ?                  |
+| snail.samples.0.2                  | [#libigl][]        |
+| snail1.dmat                        | [#libigl][]        |
+| snail2.dmat                        | [#libigl][]        |
+| snail3.dmat                        | [#libigl][]        |
+| snail4.dmat                        | [#libigl][]        |
+| sphere.obj                         | ?                  |
+| truck.obj                          | [#shrec][]         |
+| xcylinder.obj                      | [#libigl][]        |
+| ycylinder.obj                      | [#libigl][]        |
+| zcylinder.obj                      | [#libigl][]        |
+
+[#aim_at_shape]: AIM@SHAPE,IAL_SHARED_PATH "/cow.off".cnr.it:8080/ontologies/shapes/viewgroup.jsp?id=225-Fandisk_MC
+[#cosmic_blobs]: Cosmic blobs, http://www.mit.edu/~ibaran/autorig/
+[#ivan]: Hand measured from physical object, http://www.cs.utah.edu/docs/misc/Uteapot03.pdf
+[#jacobson_2013]: Alec Jacobson, Ladislav Kavan, and Olga Sorkine.
+  [Robust Inside-Outside Segmentation using Generalized Winding
+  Numbers](https://www.google.com/search?q=Robust+Inside-Outside+Segmentation+using+Generalized+Winding+Numbers),
+  2013.
+[#libigl]: Original data produced by libigl authors, http://libigl.github.io/libigl/
+[#mpi]: Max Planck Institute at Saarbrucken.
+[#pauly]: Modeled by Mark Pauly, "Shape Modeling with Point-Sampled Geometry"
+[#shrec]: SHREC 2009 Dataset,
+  http://www.itl.nist.gov/iad/vug/sharp/benchmark/shrecPartial/data.html
+[#sorkine_2004]: Olga Sorkine, Daniel Cohen-Or, [Least-squares
+  Meshes](https://www.google.com/search?q=Least+squares+meshes), 2004.
+[#stanford]: The Stanford 3D Scanning Repository, http://graphics.stanford.edu/data/3Dscanrep/

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-47cfcb9ea6dd7363e69106910fe3eaaf0aacc40b
+b4d291ea3ce115a6b9af9ffe3559ffd5f157bea6

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-129cef4e745b13b5c567053630e25f92dfead7e1
+d0e101aea39668e374ff776b29cdd79667b7d8bd