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

Former-commit-id: f89f0e0a69dc781ab6cfa0bf5ae45d6092c0aabc

.gitignore

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

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

@@ -946,4 +946,10 @@ template double igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_
 template double igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, 1, 2, 1, 1, 2> const&, int&, Eigen::Matrix<double, 1, 2, 1, 1, 2>&) const;
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&) const;
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&) const;
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::init<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::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::init<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::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template std::vector<int, std::allocator<int> > igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::find<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, bool) const;
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::serialize<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, int) const;
+template std::vector<int, std::allocator<int> > igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::find<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, bool) const;
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::serialize<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, int) const;
 #endif

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

@@ -0,0 +1,55 @@
+// 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 "connect_boundary_to_infinity.h"
+#include "boundary_facets.h"
+
+template <typename DerivedF, typename DerivedFO>
+IGL_INLINE void igl::connect_boundary_to_infinity(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedFO> & FO)
+{
+  return connect_boundary_to_infinity(F,F.maxCoeff(),FO);
+}
+template <typename DerivedF, typename DerivedFO>
+IGL_INLINE void igl::connect_boundary_to_infinity(
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const typename DerivedF::Scalar inf_index,
+  Eigen::PlainObjectBase<DerivedFO> & FO)
+{
+  // Determine boundary edges
+  Eigen::Matrix<typename DerivedFO::Scalar,Eigen::Dynamic,Eigen::Dynamic> O;
+  boundary_facets(F,O);
+  FO.resize(F.rows()+O.rows(),F.cols());
+  typedef Eigen::Matrix<typename DerivedFO::Scalar,Eigen::Dynamic,1> VectorXI;
+  FO.topLeftCorner(F.rows(),F.cols()) = F;
+  FO.bottomLeftCorner(O.rows(),O.cols()) = O.rowwise().reverse();
+  FO.bottomRightCorner(O.rows(),1).setConstant(inf_index);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedVO, 
+  typename DerivedFO>
+IGL_INLINE void igl::connect_boundary_to_infinity(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedVO> & VO,
+  Eigen::PlainObjectBase<DerivedFO> & FO)
+{
+  typename DerivedV::Index inf_index = V.rows();
+  connect_boundary_to_infinity(F,inf_index,FO);
+  VO.resize(V.rows()+1,V.cols());
+  VO.topLeftCorner(V.rows(),V.cols()) = V;
+  auto inf = std::numeric_limits<typename DerivedVO::Scalar>::infinity();
+  VO.row(V.rows()).setConstant(inf);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::connect_boundary_to_infinity<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/connect_boundary_to_infinity.h

@@ -0,0 +1,56 @@
+// 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_CONNECT_BOUNDARY_TO_INFINITY_H
+#define IGL_CONNECT_BOUNDARY_TO_INFINITY_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Connect all boundary edges to a ficticious point at infinity.
+  //
+  // Inputs:
+  //   F  #F by 3 list of face indices into some V
+  // Outputs:
+  //   FO  #F+#O by 3 list of face indices into [V;inf inf inf], original F are
+  //     guaranteed to come first. If (V,F) was a manifold mesh, now it is
+  //     closed with a possibly non-manifold vertex at infinity (but it will be
+  //     edge-manifold).
+  template <typename DerivedF, typename DerivedFO>
+  IGL_INLINE void connect_boundary_to_infinity(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedFO> & FO);
+  // Inputs:
+  //   inf_index  index of point at infinity (usually V.rows() or F.maxCoeff())
+  template <typename DerivedF, typename DerivedFO>
+  IGL_INLINE void connect_boundary_to_infinity(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const typename DerivedF::Scalar inf_index,
+    Eigen::PlainObjectBase<DerivedFO> & FO);
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of face indices into some V
+  // Outputs:
+  //   VO  #V+1 by 3 list of vertex positions, original V are guaranteed to
+  //     come first. Last point is inf, inf, inf
+  //   FO  #F+#O by 3 list of face indices into VO
+  // 
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedVO, 
+    typename DerivedFO>
+  IGL_INLINE void connect_boundary_to_infinity(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedVO> & VO,
+    Eigen::PlainObjectBase<DerivedFO> & FO);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "connect_boundary_to_infinity.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -381,8 +381,12 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   log_time("box_and_bind");
 #endif
   // Run the self intersection algorithm with all defaults
+  CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  log_time("box_intersection_d");
+#endif
   try{
-    CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
+    process_intersecting_boxes();
   }catch(int e)
   {
     // Rethrow if not IGL_FIRST_HIT_EXCEPTION
@@ -392,10 +396,6 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
-#ifdef IGL_SELFINTERSECTMESH_DEBUG
-  log_time("box_intersection_d");
-#endif
-  process_intersecting_boxes();
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("resolve_intersection");
 #endif
@@ -432,18 +432,6 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("remesh_intersection");
 #endif
-
-  // Q: Does this give the same result as TETGEN?
-  // A: For the cow and beast, yes.
-
-  // Q: Is tetgen faster than this CGAL implementation?
-  // A: Well, yes. But Tetgen is only solving the detection (predicates)
-  // problem. This is also appending the intersection objects (construction).
-  // But maybe tetgen is still faster for the detection part. For example, this
-  // CGAL implementation on the beast takes 98 seconds but tetgen detection
-  // takes 14 seconds
-
-
 }
 
 
@@ -512,6 +500,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   {
     throw IGL_FIRST_HIT_EXCEPTION;
   }
+
 }
 
 template <
@@ -869,105 +858,124 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   std::vector<std::mutex> triangle_locks(T.size());
   std::vector<std::mutex> vertex_locks(V.rows());
   std::mutex index_lock;
-  auto process_chunk = [&](const size_t first, const size_t last) -> void{
-    assert(last >= first);
+  std::mutex exception_mutex;
+  bool exception_fired = false;
+  int exception = -1;
+  auto process_chunk = 
+    [&](
+      const size_t first, 
+      const size_t last) -> void
+  {
+    try
+    {
+      assert(last >= first);
 
-    for (size_t i=first; i<last; i++) {
-      Index fa=T.size(), fb=T.size();
+      for (size_t i=first; i<last; i++) 
       {
-        // Before knowing which triangles are involved, we need to lock
-        // everything to prevent race condition in updating reference counters.
-        std::lock_guard<std::mutex> guard(index_lock);
-        const auto& box_pair = candidate_box_pairs[i];
-        const auto& a = box_pair.first;
-        const auto& b = box_pair.second;
-        fa = a.handle()-T.begin();
-        fb = b.handle()-T.begin();
-      }
-      assert(fa < T.size());
-      assert(fb < T.size());
+        if(exception_fired) return;
+        Index fa=T.size(), fb=T.size();
+        {
+          // Before knowing which triangles are involved, we need to lock
+          // everything to prevent race condition in updating reference counters.
+          std::lock_guard<std::mutex> guard(index_lock);
+          const auto& box_pair = candidate_box_pairs[i];
+          const auto& a = box_pair.first;
+          const auto& b = box_pair.second;
+          fa = a.handle()-T.begin();
+          fb = b.handle()-T.begin();
+        }
+        assert(fa < T.size());
+        assert(fb < T.size());
 
-      // Lock triangles
-      std::lock_guard<std::mutex> guard_A(triangle_locks[fa]);
-      std::lock_guard<std::mutex> guard_B(triangle_locks[fb]);
+        // Lock triangles
+        std::lock_guard<std::mutex> guard_A(triangle_locks[fa]);
+        std::lock_guard<std::mutex> guard_B(triangle_locks[fb]);
 
-      // Lock vertices
-      std::list<std::lock_guard<std::mutex> > guard_vertices;
-      {
-        std::vector<typename DerivedF::Scalar> unique_vertices;
-        std::vector<size_t> tmp1, tmp2;
-        igl::unique({F(fa,0), F(fa,1), F(fa,2), F(fb,0), F(fb,1), F(fb,2)},
-            unique_vertices, tmp1, tmp2);
-        std::for_each(unique_vertices.begin(), unique_vertices.end(),
-            [&](const typename DerivedF::Scalar& vi) {
-            guard_vertices.emplace_back(vertex_locks[vi]);
-            });
-      }
-
-      const Triangle_3& A = T[fa];
-      const Triangle_3& B = T[fb];
-
-      // Number of combinatorially shared vertices
-      Index comb_shared_vertices = 0;
-      // Number of geometrically shared vertices (*not* including combinatorially
-      // shared)
-      Index geo_shared_vertices = 0;
-      // Keep track of shared vertex indices
-      std::vector<std::pair<Index,Index> > shared;
-      Index ea,eb;
-      for(ea=0;ea<3;ea++)
-      {
-        for(eb=0;eb<3;eb++)
+        // Lock vertices
+        std::list<std::lock_guard<std::mutex> > guard_vertices;
         {
-          if(F(fa,ea) == F(fb,eb))
-          {
-            comb_shared_vertices++;
-            shared.emplace_back(ea,eb);
-          }else if(A.vertex(ea) == B.vertex(eb))
+          std::vector<typename DerivedF::Scalar> unique_vertices;
+          std::vector<size_t> tmp1, tmp2;
+          igl::unique({F(fa,0), F(fa,1), F(fa,2), F(fb,0), F(fb,1), F(fb,2)},
+              unique_vertices, tmp1, tmp2);
+          std::for_each(unique_vertices.begin(), unique_vertices.end(),
+              [&](const typename DerivedF::Scalar& vi) {
+              guard_vertices.emplace_back(vertex_locks[vi]);
+              });
+        }
+        if(exception_fired) return;
+
+        const Triangle_3& A = T[fa];
+        const Triangle_3& B = T[fb];
+
+        // Number of combinatorially shared vertices
+        Index comb_shared_vertices = 0;
+        // Number of geometrically shared vertices (*not* including combinatorially
+        // shared)
+        Index geo_shared_vertices = 0;
+        // Keep track of shared vertex indices
+        std::vector<std::pair<Index,Index> > shared;
+        Index ea,eb;
+        for(ea=0;ea<3;ea++)
+        {
+          for(eb=0;eb<3;eb++)
           {
-            geo_shared_vertices++;
-            shared.emplace_back(ea,eb);
+            if(F(fa,ea) == F(fb,eb))
+            {
+              comb_shared_vertices++;
+              shared.emplace_back(ea,eb);
+            }else if(A.vertex(ea) == B.vertex(eb))
+            {
+              geo_shared_vertices++;
+              shared.emplace_back(ea,eb);
+            }
           }
         }
-      }
-      const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
+        const Index total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
+        if(exception_fired) return;
 
-      if(comb_shared_vertices== 3)
-      {
-        assert(shared.size() == 3);
-        //// Combinatorially duplicate face, these should be removed by preprocessing
-        //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
-        continue;
-      }
-      if(total_shared_vertices== 3)
-      {
-        assert(shared.size() == 3);
-        //// Geometrically duplicate face, these should be removed by preprocessing
-        //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
-        continue;
-      }
-      if(total_shared_vertices == 2)
-      {
-        assert(shared.size() == 2);
-        // Q: What about coplanar?
-        //
-        // o    o
-        // |\  /|
-        // | \/ |
-        // | /\ |
-        // |/  \|
-        // o----o
-        double_shared_vertex(A,B,fa,fb,shared);
-        continue;
-      }
-      assert(total_shared_vertices<=1);
-      if(total_shared_vertices==1)
-      {
-        single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
-      }else
-      {
-        intersect(A,B,fa,fb);
+        if(comb_shared_vertices== 3)
+        {
+          assert(shared.size() == 3);
+          //// Combinatorially duplicate face, these should be removed by preprocessing
+          //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
+          continue;
+        }
+        if(total_shared_vertices== 3)
+        {
+          assert(shared.size() == 3);
+          //// Geometrically duplicate face, these should be removed by preprocessing
+          //cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
+          continue;
+        }
+        if(total_shared_vertices == 2)
+        {
+          assert(shared.size() == 2);
+          // Q: What about coplanar?
+          //
+          // o    o
+          // |\  /|
+          // | \/ |
+          // | /\ |
+          // |/  \|
+          // o----o
+          double_shared_vertex(A,B,fa,fb,shared);
+          continue;
+        }
+        assert(total_shared_vertices<=1);
+        if(total_shared_vertices==1)
+        {
+          single_shared_vertex(A,B,fa,fb,shared[0].first,shared[0].second);
+        }else
+        {
+          intersect(A,B,fa,fb);
+        }
       }
+    }catch(int e)
+    {
+      std::lock_guard<std::mutex> exception_lock(exception_mutex);
+      exception_fired = true;
+      exception = e;
     }
   };
   size_t num_threads=0;
@@ -982,14 +990,17 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   const size_t num_pairs = candidate_box_pairs.size();
   const size_t chunk_size = num_pairs / num_threads;
   std::vector<std::thread> threads;
-  for (size_t i=0; i<num_threads-1; i++) {
+  for (size_t i=0; i<num_threads-1; i++) 
+  {
     threads.emplace_back(process_chunk, i*chunk_size, (i+1)*chunk_size);
   }
   // Do some work in the master thread.
   process_chunk((num_threads-1)*chunk_size, num_pairs);
-  for (auto& t : threads) {
+  for (auto& t : threads) 
+  {
     if (t.joinable()) t.join();
   }
+  if(exception_fired) throw exception;
   //process_chunk(0, candidate_box_pairs.size());
 }
 

include/igl/copyleft/tetgen/tetrahedralize.cpp

@@ -100,7 +100,120 @@ IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
   return e;
 }
 
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedVM, 
+  typename DerivedFM, 
+  typename DerivedTV, 
+  typename DerivedTT, 
+  typename DerivedTF,
+  typename DerivedTM>
+IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedVM>& VM,
+  const Eigen::PlainObjectBase<DerivedFM>& FM,
+  const std::string switches,
+  Eigen::PlainObjectBase<DerivedTV>& TV,
+  Eigen::PlainObjectBase<DerivedTT>& TT,
+  Eigen::PlainObjectBase<DerivedTF>& TF,
+  Eigen::PlainObjectBase<DerivedTM>& TM)
+{
+  using namespace std;
+  vector<vector<REAL> > vV,vTV;
+  vector<vector<int> > vF,vTT,vTF;
+  vector<int> vTM;
+	
+  matrix_to_list(V,vV);
+  matrix_to_list(F,vF);
+	vector<int> vVM = matrix_to_list(VM);
+	vector<int> vFM = matrix_to_list(FM);
+  int e = tetrahedralize(vV,vF,vVM,vFM,switches,vTV,vTT,vTF,vTM);
+  if(e == 0)
+  {
+    bool TV_rect = list_to_matrix(vTV,TV);
+    if(!TV_rect)
+    {
+      return false;
+    }
+    bool TT_rect = list_to_matrix(vTT,TT);
+    if(!TT_rect)
+    {
+      return false;
+    }
+    bool TF_rect = list_to_matrix(vTF,TF);
+    if(!TF_rect)
+    {
+      return false;
+    }
+    bool TM_rect = list_to_matrix(vTM,TM);
+    if(!TM_rect)
+    {
+      return false;
+    }
+  }
+  return e;
+}
+IGL_INLINE int igl::copyleft::tetgen::tetrahedralize(
+  const std::vector<std::vector<REAL > > & V, 
+  const std::vector<std::vector<int> > & F, 
+  const std::vector<int> & VM, 
+	const std::vector<int> & FM,
+  const std::string switches,
+  std::vector<std::vector<REAL > > & TV, 
+  std::vector<std::vector<int > > & TT, 
+  std::vector<std::vector<int> > & TF,
+  std::vector<int> & TM)
+{
+  using namespace std;
+  tetgenio in,out;
+  bool success;
+  success = mesh_to_tetgenio(V,F,in);
+  if(!success)
+  {
+    return -1;
+  }
+	in.pointmarkerlist = new int[VM.size()];
+	for (int i = 0; i < VM.size(); ++i) {
+		in.pointmarkerlist[i] = VM[i];
+	}
+  // These have already been created in mesh_to_tetgenio.
+  // Reset them here.
+	for (int i = 0; i < FM.size(); ++i) {
+		in.facetmarkerlist[i] = FM[i];
+	}
+  try
+  {
+    char * cswitches = new char[switches.size() + 1];
+    std::strcpy(cswitches,switches.c_str());
+    ::tetrahedralize(cswitches,&in, &out);
+    delete[] cswitches;
+  }catch(int e)
+  {
+    cerr<<"^"<<__FUNCTION__<<": TETGEN CRASHED... KABOOOM!!!"<<endl;
+    return 1;
+  }
+  if(out.numberoftetrahedra == 0)
+  {
+    cerr<<"^"<<__FUNCTION__<<": Tetgen failed to create tets"<<endl;
+    return 2;
+  }
+  success = tetgenio_to_tetmesh(out,TV,TT,TF);
+  if(!success)
+  {
+    return -1;
+  }
+	TM.resize(out.numberofpoints);
+	for (int i = 0; i < out.numberofpoints; ++i) {
+		TM[i] = out.pointmarkerlist[i];
+	}
+  //boundary_facets(TT,TF);
+  return 0;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template int igl::copyleft::tetgen::tetrahedralize<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::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&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template int igl::copyleft::tetgen::tetrahedralize<Eigen::Matrix<double, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, 1, 0, -1, 1>,Eigen::Matrix<int, -1, 1, 0, -1, 1>,Eigen::Matrix<double, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, -1, 0, -1, -1>,Eigen::Matrix<int, -1, 1, 0, -1, 1> >(const Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &,const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &,const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > &,const Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > &,const std::basic_string<char, std::char_traits<char>, std::allocator<char> >,Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &,Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &,Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > &);
 #endif

include/igl/copyleft/tetgen/tetrahedralize.h

@@ -67,6 +67,62 @@ namespace igl
         Eigen::PlainObjectBase<DerivedTV>& TV,
         Eigen::PlainObjectBase<DerivedTT>& TT,
         Eigen::PlainObjectBase<DerivedTF>& TF);
+      
+			// Mesh the interior of a surface mesh (V,F) using tetgen
+      //
+      // Inputs:
+      //   V  #V by 3 vertex position list
+      //   F  #F list of polygon face indices into V (0-indexed)
+			//   M  #V list of markers for vertices
+      //   switches  string of tetgen options (See tetgen documentation) e.g.
+      //     "pq1.414a0.01" tries to mesh the interior of a given surface with
+      //       quality and area constraints
+      //     "" will mesh the convex hull constrained to pass through V (ignores F)
+      // Outputs:
+      //   TV  #V by 3 vertex position list
+      //   TT  #T by 4 list of tet face indices
+      //   TF  #F by 3 list of triangle face indices
+			//   TM  #V list of markers for vertices
+      // Returns status:
+      //   0 success
+      //   1 tetgen threw exception
+      //   2 tetgen did not crash but could not create any tets (probably there are
+      //     holes, duplicate faces etc.)
+      //   -1 other error
+      IGL_INLINE int tetrahedralize(
+        const std::vector<std::vector<REAL > > & V, 
+        const std::vector<std::vector<int> > & F, 
+				const std::vector<int> & VM,
+				const std::vector<int> & FM,
+        const std::string switches,
+        std::vector<std::vector<REAL > > & TV, 
+        std::vector<std::vector<int > > & TT, 
+        std::vector<std::vector<int> > & TF,
+				std::vector<int> & TM);
+      
+      // Wrapper with Eigen types
+      // Templates:
+      //   DerivedV  real-value: i.e. from MatrixXd
+      //   DerivedF  integer-value: i.e. from MatrixXi
+      template <
+        typename DerivedV, 
+        typename DerivedF, 
+				typename DerivedVM,
+				typename DerivedFM,
+        typename DerivedTV, 
+        typename DerivedTT, 
+        typename DerivedTF, 
+        typename DerivedTM>
+      IGL_INLINE int tetrahedralize(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedVM>& VM,
+        const Eigen::PlainObjectBase<DerivedFM>& FM,
+        const std::string switches,
+        Eigen::PlainObjectBase<DerivedTV>& TV,
+        Eigen::PlainObjectBase<DerivedTT>& TT,
+        Eigen::PlainObjectBase<DerivedTF>& TF,
+        Eigen::PlainObjectBase<DerivedTM>& TM);
     }
   }
 }

include/igl/decimate.cpp

@@ -1,6 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 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 
@@ -9,7 +9,8 @@
 #include "collapse_edge.h"
 #include "edge_flaps.h"
 #include "remove_unreferenced.h"
-#include "max_faces_stopping_condition.h"
+#include "slice_mask.h"
+#include "connect_boundary_to_infinity.h"
 #include <iostream>
 
 IGL_INLINE bool igl::decimate(
@@ -20,6 +21,9 @@ IGL_INLINE bool igl::decimate(
   Eigen::MatrixXi & G,
   Eigen::VectorXi & J)
 {
+  // Original number of faces
+  const int orig_m = F.rows();
+  // Tracking number of faces
   int m = F.rows();
   const auto & shortest_edge_and_midpoint = [](
     const int e,
@@ -35,14 +39,47 @@ IGL_INLINE bool igl::decimate(
     cost = (V.row(E(e,0))-V.row(E(e,1))).norm();
     p = 0.5*(V.row(E(e,0))+V.row(E(e,1)));
   };
-  return decimate(
-    V,
-    F,
+  typedef Eigen::MatrixXd DerivedV;
+  typedef Eigen::MatrixXi DerivedF;
+  DerivedV VO;
+  DerivedF FO;
+  igl::connect_boundary_to_infinity(V,F,VO,FO);
+  const auto & max_non_infinite_faces_stopping_condition = 
+    [max_m,orig_m,&m](
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const std::set<std::pair<double,int> > &,
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,
+      const Eigen::MatrixXd &,
+      const int,
+      const int,
+      const int,
+      const int f1,
+      const int f2) -> bool
+    {
+      // Only subtract if we're collapsing a real face
+      if(f1 < orig_m) m-=1;
+      if(f2 < orig_m) m-=1;
+      return m<=(int)max_m;
+    };
+  bool ret = decimate(
+    VO,
+    FO,
     shortest_edge_and_midpoint,
-    max_faces_stopping_condition(m,max_m),
+    max_non_infinite_faces_stopping_condition,
     U,
     G,
     J);
+  const Eigen::Array<bool,Eigen::Dynamic,1> keep = (J.array()<orig_m);
+  igl::slice_mask(Eigen::MatrixXi(G),keep,1,G);
+  igl::slice_mask(Eigen::VectorXi(J),keep,1,J);
+  Eigen::VectorXi _1;
+  igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1);
+  return ret;
 }
 
 IGL_INLINE bool igl::decimate(

include/igl/decimate.h

@@ -1,6 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+// 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 
@@ -13,10 +13,9 @@
 #include <set>
 namespace igl
 {
-  // Assumes (V,F) is a closed manifold mesh 
-  // Collapses edges until desired number of faces is achieved. This uses
-  // default edge cost and merged vertex placement functions {edge length, edge
-  // midpoint}.
+  // Assumes (V,F) is a manifold mesh (possibly with boundary) Collapses edges
+  // until desired number of faces is achieved. This uses default edge cost and
+  // merged vertex placement functions {edge length, edge midpoint}.
   //
   // Inputs:
   //   V  #V by dim list of vertex positions
@@ -34,6 +33,11 @@ namespace igl
     Eigen::MatrixXd & U,
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J);
+  // Assumes a **closed** manifold mesh. See igl::connect_boundary_to_infinity
+  // and igl::decimate in decimate.cpp
+  // is handling meshes with boundary by connecting all boundary edges with
+  // dummy facets to infinity **and** modifying the stopping criteria.
+  //
   // Inputs:
   //   cost_and_placement  function computing cost of collapsing an edge and 3d
   //     position where it should be placed:
@@ -47,30 +51,32 @@ namespace igl
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
     const std::function<void(
-      const int,
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      double &,
-      Eigen::RowVectorXd &)> & cost_and_placement,
+      const int              /*e*/,
+      const Eigen::MatrixXd &/*V*/,
+      const Eigen::MatrixXi &/*F*/,
+      const Eigen::MatrixXi &/*E*/,
+      const Eigen::VectorXi &/*EMAP*/,
+      const Eigen::MatrixXi &/*EF*/,
+      const Eigen::MatrixXi &/*EI*/,
+      double &               /*cost*/,
+      Eigen::RowVectorXd &   /*p*/
+      )> & cost_and_placement,
     const std::function<bool(
-      const Eigen::MatrixXd &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::VectorXi &,
-      const Eigen::MatrixXi &,
-      const Eigen::MatrixXi &,
-      const std::set<std::pair<double,int> > &,
-      const std::vector<std::set<std::pair<double,int> >::iterator > &,
-      const Eigen::MatrixXd &,
-      const int,
-      const int,
-      const int,
-      const int,
-      const int)> & stopping_condition,
+      const Eigen::MatrixXd &                                         ,/*V*/
+      const Eigen::MatrixXi &                                         ,/*F*/
+      const Eigen::MatrixXi &                                         ,/*E*/
+      const Eigen::VectorXi &                                         ,/*EMAP*/
+      const Eigen::MatrixXi &                                         ,/*EF*/
+      const Eigen::MatrixXi &                                         ,/*EI*/
+      const std::set<std::pair<double,int> > &                        ,/*Q*/
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
+      const Eigen::MatrixXd &                                         ,/*C*/
+      const int                                                       ,/*e*/
+      const int                                                       ,/*e1*/
+      const int                                                       ,/*e2*/
+      const int                                                       ,/*f1*/
+      const int                                                        /*f2*/
+      )> & stopping_condition,
     Eigen::MatrixXd & U,
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J);

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

@@ -118,4 +118,5 @@ template void igl::per_face_normals_stable<Eigen::Matrix<double, -1, 3, 0, -1, 3
 template void igl::per_face_normals_stable<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_face_normals_stable<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> >&);
 #endif

include/igl/png/render_to_png.cpp

@@ -17,7 +17,7 @@ IGL_INLINE bool igl::png::render_to_png(
   const bool alpha,
   const bool fast)
 {
-  unsigned char * data = new unsigned char[width*height];
+  unsigned char * data = new unsigned char[4*width*height];
   glReadPixels(
     0,
     0,
@@ -35,7 +35,7 @@ IGL_INLINE bool igl::png::render_to_png(
       data[4*(i+j*width)+3] = 255;
     }
   }
-  bool ret = stbi_write_png(png_file.c_str(), width, height, 4, data, width*sizeof(unsigned char));
+  bool ret = stbi_write_png(png_file.c_str(), width, height, 4, data, 4*width*sizeof(unsigned char));
   delete [] data;
   return ret;
 }

include/igl/seam_edges.cpp

@@ -0,0 +1,180 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.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 "seam_edges.h"
+#include <unordered_map>
+#include <unordered_set>
+#include <cassert>
+
+namespace {
+    // Computes the orientation of `c` relative to the line between `a` and `b`.
+    // Assumes 2D vector input.
+    // Based on: https://www.cs.cmu.edu/~quake/robust.html
+    template< typename scalar_t >
+    inline scalar_t orientation(
+        const Eigen::Matrix< scalar_t, 2, 1 >& a,
+        const Eigen::Matrix< scalar_t, 2, 1 >& b,
+        const Eigen::Matrix< scalar_t, 2, 1 >& c
+        ) {
+        typedef Eigen::Matrix< scalar_t, 2, 1 > Vector2S;
+        const Vector2S row0 = a - c;
+        const Vector2S row1 = b - c;
+        return row0(0)*row1(1) - row1(0)*row0(1);
+    }
+}
+
+// I have verified that this function produces the exact same output as
+// `find_seam_fast.py` for `cow_triangled.obj`.
+template <typename DerivedV, typename DerivedF, typename DerivedT>
+IGL_INLINE void seam_edges(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& TC,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedF>& FTC,
+    Eigen::PlainObjectBase<DerivedF>& seams,
+    Eigen::PlainObjectBase<DerivedF>& boundaries,
+    Eigen::PlainObjectBase<DerivedF>& foldovers
+    )
+{
+    // Assume triangles.
+    assert( F.cols() == 3 );
+    assert( F.cols() == FTC.cols() );
+    assert( F.rows() == FTC.rows() );
+    
+    // Assume 2D texture coordinates (foldovers tests).
+    assert( TC.cols() == 2 );
+    typedef Eigen::Matrix< typename DerivedT::Scalar, 2, 1 > Vector2S;
+    
+    seams     .setZero( 3*F.rows(), 4 );
+    boundaries.setZero( 3*F.rows(), 2 );
+    foldovers .setZero( 3*F.rows(), 4 );
+    
+    int num_seams = 0;
+    int num_boundaries = 0;
+    int num_foldovers = 0;
+    
+    // A map from a pair of vertex indices to the index (face and endpoints)
+    // into face_position_indices.
+    // The following should be true for every key, value pair:
+    //    key == face_position_indices[ value ]
+    // This gives us a "reverse map" so that we can look up other face
+    // attributes based on position edges.
+    // The value are written in the format returned by numpy.where(),
+    // which stores multi-dimensional indices such as array[a0,b0], array[a1,b1]
+    // as ( (a0,a1), (b0,b1) ).
+    
+    // We need to make a hash function for our directed edges.
+    // We'll use i*V.rows() + j.
+    typedef std::pair< typename DerivedF::Scalar, typename DerivedF::Scalar > directed_edge;
+	const int numV = V.rows();
+	const int numF = F.rows();
+	auto edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
+	// When we pass a hash function object, we also need to specify the number of buckets.
+	// The Euler characteristic says that the number of undirected edges is numV + numF -2*genus.
+	std::unordered_map< directed_edge, std::pair< int, int >, decltype( edge_hasher ) > directed_position_edge2face_position_index( 2*( numV + numF ), edge_hasher );
+    for( int fi = 0; fi < F.rows(); ++fi ) {
+        for( int i = 0; i < 3; ++i ) {
+            const int j = ( i+1 ) % 3;
+            directed_position_edge2face_position_index[ std::make_pair( F(fi,i), F(fi,j) ) ] = std::make_pair( fi, i );
+        }
+    }
+    
+    // First find all undirected position edges (collect both orientations of
+    // the directed edges).
+    std::unordered_set< directed_edge, decltype( edge_hasher ) > undirected_position_edges( numV + numF, edge_hasher );
+    for( const auto& el : directed_position_edge2face_position_index ) {
+        undirected_position_edges.insert( el.first );
+        undirected_position_edges.insert( std::make_pair( el.first.second, el.first.first ) );
+    }
+    
+    // Now we will iterate over all position edges.
+    // Seam edges are the edges whose two opposite directed edges have different
+    // texcoord indices (or one doesn't exist at all in the case of a mesh
+    // boundary).
+    for( const auto& vp_edge : undirected_position_edges ) {
+        const auto vp_edge_reverse = std::make_pair( vp_edge.second, vp_edge.first );
+        // If it and its opposite exist as directed edges, check if their
+        // texture coordinate indices match.
+        if( directed_position_edge2face_position_index.count( vp_edge ) &&
+            directed_position_edge2face_position_index.count( vp_edge_reverse ) ) {
+            const auto forwards = directed_position_edge2face_position_index[ vp_edge ];
+            const auto backwards = directed_position_edge2face_position_index[ vp_edge_reverse ];
+            
+            // NOTE: They should never be equal.
+            assert( forwards != backwards );
+            
+            // NOTE: Non-matching seam edges will show up twice, once as
+            //       ( forwards, backwards ) and once as
+            //       ( backwards, forwards ). We don't need to examine both,
+            //       so continue only if forwards < backwards.
+            if( forwards < backwards ) continue;
+
+            // If the texcoord indices match (are similarly flipped),
+            // this edge is not a seam. It could be a foldover.
+            if( std::make_pair( FTC( forwards.first, forwards.second ), FTC( forwards.first, ( forwards.second+1 ) % 3 ) )
+                ==
+                std::make_pair( FTC( backwards.first, ( backwards.second+1 ) % 3 ), FTC( backwards.first, backwards.second ) )
+                ) {
+                
+                // Check for foldovers in UV space.
+                // Get the edge (a,b) and the two opposite vertices's texture coordinates.
+                const Vector2S a = TC.row( FTC( forwards.first,  forwards.second ) );
+                const Vector2S b = TC.row( FTC( forwards.first, (forwards.second+1) % 3 ) );
+                const Vector2S c_forwards  = TC.row( FTC( forwards .first, (forwards .second+2) % 3 ) );
+                const Vector2S c_backwards = TC.row( FTC( backwards.first, (backwards.second+2) % 3 ) );
+                // If the opposite vertices' texture coordinates fall on the same side
+                // of the edge, we have a UV-space foldover.
+                const auto orientation_forwards = orientation( a, b, c_forwards );
+                const auto orientation_backwards = orientation( a, b, c_backwards );
+                if( ( orientation_forwards > 0 && orientation_backwards > 0 ) ||
+                    ( orientation_forwards < 0 && orientation_backwards < 0 )
+                    ) {
+                    foldovers( num_foldovers, 0 ) = forwards.first;
+                    foldovers( num_foldovers, 1 ) = forwards.second;
+                    foldovers( num_foldovers, 2 ) = backwards.first;
+                    foldovers( num_foldovers, 3 ) = backwards.second;
+                    num_foldovers += 1;
+                }
+            }
+            
+            // Otherwise, we have a non-matching seam edge.
+            else {
+                seams( num_seams, 0 ) = forwards.first;
+                seams( num_seams, 1 ) = forwards.second;
+                seams( num_seams, 2 ) = backwards.first;
+                seams( num_seams, 3 ) = backwards.second;
+                num_seams += 1;
+            }
+        }
+        // Otherwise, the edge and its opposite aren't both in the directed
+        // edges. One of them should be.
+        else if( directed_position_edge2face_position_index.count( vp_edge ) ) {
+            const auto forwards = directed_position_edge2face_position_index[ vp_edge ];
+            boundaries( num_boundaries, 0 ) = forwards.first;
+            boundaries( num_boundaries, 1 ) = forwards.second;
+            num_boundaries += 1;
+        }
+        else if( directed_position_edge2face_position_index.count( vp_edge_reverse ) ) {
+            const auto backwards = directed_position_edge2face_position_index[ vp_edge_reverse ];
+            boundaries( num_boundaries, 0 ) = backwards.first;
+            boundaries( num_boundaries, 1 ) = backwards.second;
+            num_boundaries += 1;
+        }
+        else {
+            // This should never happen! One of these two must have been seen.
+            assert(
+                directed_position_edge2face_position_index.count( vp_edge ) ||
+                directed_position_edge2face_position_index.count( vp_edge_reverse )
+                );
+        }
+    }
+    
+    seams     .conservativeResize( num_seams,      Eigen::NoChange_t() );
+    boundaries.conservativeResize( num_boundaries, Eigen::NoChange_t() );
+    foldovers .conservativeResize( num_foldovers,  Eigen::NoChange_t() );
+}

include/igl/seam_edges.h

@@ -0,0 +1,57 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Yotam Gingold <yotam@yotamgingold.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_SEAM_EDGES_H
+#define IGL_SEAM_EDGES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+    /*
+    Returns all UV-space boundaries of a mesh.
+    Inputs:
+        V: The positions of the input mesh.
+        TC: The 2D texture coordinates of the input mesh (2 columns).
+        F: A manifold-with-boundary triangle mesh, as vertex indices into `V` for the three vertices of each triangle.
+        FTC: Indices into `TC` for the three vertices of each triangle.
+    Outputs:
+        seams: Edges where the forwards and backwards directions have different texture
+               coordinates, as a #seams-by-4 matrix of indices.
+               Each row is organized as [ forward_face_index, forward_face_vertex_index,
+               backwards_face_index, backwards_face_vertex_index ]
+               such that one side of the seam is the edge:
+                    F[ seams( i, 0 ), seams( i, 1 ) ], F[ seams( i, 0 ), (seams( i, 1 ) + 1) % 3 ]
+               and the other side is the edge:
+                    F[ seams( i, 2 ), seams( i, 3 ) ], F[ seams( i, 2 ), (seams( i, 3 ) + 1) % 3 ]
+        boundaries: Edges with only one incident triangle, as a #boundaries-by-2 matrix of
+                    indices. Each row is organized as [ face_index, face_vertex_index ]
+                    such that the edge is:
+                        F[ boundaries( i, 0 ), boundaries( i, 1 ) ], F[ boundaries( i, 0 ), (boundaries( i, 1 ) + 1) % 3 ]
+        foldovers: Edges where the two incident triangles fold over each other in UV-space,
+                   as a #foldovers-by-4 matrix of indices.
+                   Each row is organized as [ forward_face_index, forward_face_vertex_index,
+                   backwards_face_index, backwards_face_vertex_index ]
+                   such that one side of the foldover is the edge:
+                        F[ foldovers( i, 0 ), foldovers( i, 1 ) ], F[ foldovers( i, 0 ), (foldovers( i, 1 ) + 1) % 3 ]
+                   and the other side is the edge:
+                        F[ foldovers( i, 2 ), foldovers( i, 3 ) ], F[ foldovers( i, 2 ), (foldovers( i, 3 ) + 1) % 3 ]
+    */
+    template <typename DerivedV, typename DerivedF, typename DerivedT>
+    IGL_INLINE void seam_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedT>& TC,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedF>& FTC,
+        Eigen::PlainObjectBase<DerivedF>& seams,
+        Eigen::PlainObjectBase<DerivedF>& boundaries,
+        Eigen::PlainObjectBase<DerivedF>& foldovers
+        );
+}
+#ifndef IGL_STATIC_LIBRARY
+#   include "seam_edges.cpp"
+#endif
+#endif // IGL_SEAM_EDGES_H

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

@@ -120,5 +120,5 @@ template void igl::slice_mask<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::Pla
 template void igl::slice_mask<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice_mask<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice_mask<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-
+template void igl::slice_mask<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::Array<bool, -1, 1, 0, -1, 1> const&, int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

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/triangle/triangulate.cpp

@@ -129,6 +129,107 @@ IGL_INLINE void igl::triangle::triangulate(
 
 }
 
+// Allows use of markers
+IGL_INLINE void igl::triangle::triangulate(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& E,
+  const Eigen::MatrixXd& H,
+  const Eigen::VectorXi& VM,
+  const Eigen::VectorXi& EM,
+  const std::string flags,
+  Eigen::MatrixXd& V2,
+  Eigen::MatrixXi& F2,
+  Eigen::VectorXi& M2)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // Prepare the flags
+  string full_flags = flags + "pz";
+
+  // Prepare the input struct
+  triangulateio in;
+
+  assert(V.cols() == 2);
+
+  in.numberofpoints = V.rows();
+  in.pointlist = (double*)calloc(V.rows()*2,sizeof(double));
+  for (unsigned i=0;i<V.rows();++i)
+    for (unsigned j=0;j<2;++j)
+      in.pointlist[i*2+j] = V(i,j);
+
+  in.numberofpointattributes = 0;
+  in.pointmarkerlist = (int*)calloc(VM.rows(),sizeof(int));
+
+   for (unsigned i=0;i<VM.rows();++i) {
+     in.pointmarkerlist[i] = VM(i);
+   }
+
+  in.trianglelist = NULL;
+  in.numberoftriangles = 0;
+  in.numberofcorners = 0;
+  in.numberoftriangleattributes = 0;
+  in.triangleattributelist = NULL;
+
+  in.numberofsegments = E.rows();
+  in.segmentlist = (int*)calloc(E.rows()*2,sizeof(int));
+  for (unsigned i=0;i<E.rows();++i)
+    for (unsigned j=0;j<2;++j)
+      in.segmentlist[i*2+j] = E(i,j);
+  in.segmentmarkerlist = (int*)calloc(E.rows(),sizeof(int));
+  for (unsigned i=0;i<E.rows();++i)
+    in.segmentmarkerlist[i] = EM(i);
+
+  in.numberofholes = H.rows();
+  in.holelist = (double*)calloc(H.rows()*2,sizeof(double));
+  for (unsigned i=0;i<H.rows();++i)
+    for (unsigned j=0;j<2;++j)
+      in.holelist[i*2+j] = H(i,j);
+  in.numberofregions = 0;
+
+  // Prepare the output struct
+  triangulateio out;
+
+  out.pointlist = NULL;
+  out.trianglelist = NULL;
+  out.segmentlist = NULL;
+  out.segmentmarkerlist = NULL;
+  out.pointmarkerlist = NULL;
+
+  // Call triangle
+  ::triangulate(const_cast<char*>(full_flags.c_str()), &in, &out, 0);
+
+  // Return the mesh
+  V2.resize(out.numberofpoints,2);
+  for (unsigned i=0;i<V2.rows();++i)
+    for (unsigned j=0;j<2;++j)
+      V2(i,j) = out.pointlist[i*2+j];
+
+  F2.resize(out.numberoftriangles,3);
+  for (unsigned i=0;i<F2.rows();++i)
+    for (unsigned j=0;j<3;++j)
+      F2(i,j) = out.trianglelist[i*3+j];
+
+  M2.resize(out.numberofpoints);
+  for (unsigned int i = 0; i < out.numberofpoints; ++i) {
+    M2(i) = out.pointmarkerlist[i];
+  }
+
+  // Cleanup in
+  free(in.pointlist);
+  free(in.pointmarkerlist);
+  free(in.segmentlist);
+  free(in.segmentmarkerlist);
+  free(in.holelist);
+
+  // Cleanup out
+  free(out.pointlist);
+  free(out.trianglelist);
+  free(out.segmentlist);
+  free(out.segmentmarkerlist);
+  free(out.pointmarkerlist);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 #endif

include/igl/triangle/triangulate.h

@@ -35,6 +35,32 @@ namespace igl
       const std::string flags,
       Eigen::MatrixXd& V2,
       Eigen::MatrixXi& F2);
+		
+		// Triangulate the interior of a polygon using the triangle library.
+    //
+    // Inputs:
+    //   V #V by 2 list of 2D vertex positions
+    //   E #E by 2 list of vertex ids forming unoriented edges of the boundary of the polygon
+    //   H #H by 2 coordinates of points contained inside holes of the polygon
+		//   M #V list of markers for input vertices
+    //   flags  string of options pass to triangle (see triangle documentation)
+    // Outputs:
+    //   V2  #V2 by 2  coordinates of the vertives of the generated triangulation
+    //   F2  #F2 by 3  list of indices forming the faces of the generated triangulation
+		//   M2  #V2 list of markers for output vertices
+    //
+    // TODO: expose the option to prevent Steiner points on the boundary
+    //
+		IGL_INLINE void triangulate(
+			const Eigen::MatrixXd& V,
+			const Eigen::MatrixXi& E,
+			const Eigen::MatrixXd& H,
+			const Eigen::VectorXi& VM,
+			const Eigen::VectorXi& EM,
+			const std::string flags,
+			Eigen::MatrixXd& V2,
+			Eigen::MatrixXi& F2,
+			Eigen::VectorXi& M2);
   }
 }
 

include/igl/viewer/ViewerData.cpp

@@ -12,6 +12,7 @@
 
 #include <igl/per_face_normals.h>
 #include <igl/material_colors.h>
+#include <igl/parula.h>
 #include <igl/per_vertex_normals.h>
 
 #ifdef ENABLE_SERIALIZATION
@@ -154,6 +155,12 @@ IGL_INLINE void igl::viewer::ViewerData::set_colors(const Eigen::MatrixXd &C)
 {
   using namespace std;
   using namespace Eigen;
+  if(C.rows()>0 && C.cols() == 1)
+  {
+    Eigen::MatrixXd C3;
+    igl::parula(C,true,C3);
+    return set_colors(C3);
+  }
   // Ambient color should be darker color
   const auto ambient = [](const MatrixXd & C)->MatrixXd
   {

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

@@ -378,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
     //
@@ -475,12 +482,12 @@ const char *__doc_igl_get_seconds = R"igl_Qu8mg5v7(// Return the current time in
   //    ... // 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
+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".
   //
@@ -981,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:

python/py_doc.h

@@ -26,6 +26,7 @@ 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;
@@ -80,6 +81,9 @@ 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;

python/py_igl.cpp

@@ -25,6 +25,7 @@
 #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>
@@ -68,6 +69,8 @@
 #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>
@@ -101,6 +104,7 @@ 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"
@@ -144,6 +148,8 @@ 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"

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_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_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/python_shared.cpp

@@ -77,6 +77,7 @@ PYBIND11_PLUGIN(pyigl) {
            cut_mesh_from_singularities
            doublearea
            edge_lengths
+           edge_topology
            eigs
            embree_ambient_occlusion
            embree_reorient_facets_raycast
@@ -124,6 +125,9 @@ PYBIND11_PLUGIN(pyigl) {
            slice_mask
            slice_tets
            sortrows
+           streamlines_init
+           streamlines_next
+           triangle_triangle_adjacency
            triangle_triangulate
            unique
            unproject_onto_mesh

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

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

@@ -181,4 +181,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/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-85350ebcb5f876161287c77568f3b6900d02b707
+b4d291ea3ce115a6b9af9ffe3559ffd5f157bea6

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-24049c9ad2cc3eff90b021f5f1d017c8c042e6da
+d0e101aea39668e374ff776b29cdd79667b7d8bd