Merge branch 'master' of github.com:libigl/libigl into alecjacobson

Former-commit-id: 859ad326e6bf75d98edf3cdb271330774e200845

.gitignore

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

include/igl/AABB.cpp

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

include/igl/copyleft/cgal/assign_scalar.cpp

@@ -38,10 +38,10 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
-      const float next = nextafter(d, interval.second);
+      const float next = nextafter(d, float(interval.second));
       if (CGAL::abs(cgal-d) < CGAL::abs(cgal-next)) break;
       d = next;
-  } while (d < interval.second);
+  } while (d < float(interval.second));
 }
 
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(

include/igl/copyleft/tetgen/tetrahedralize.cpp

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

include/igl/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_ids.cpp

@@ -0,0 +1,45 @@
+// 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_ids.h"
+
+#include <igl/edge_topology.h>
+
+template <typename Scalar, typename Index>
+IGL_INLINE Eigen::VectorXi flipped_triangles_ids(
+  const Eigen::PlainObjectBase<Scalar> & V,
+  const Eigen::PlainObjectBase<Index> & F
+)
+{
+  assert(V.cols() == 2);
+  std::vector<int> flip_idx;
+  for (int i = 0; i < F.rows(); i++) {
+
+    Eigen::Vector2d v1_n = V.row(F(i,0)); Eigen::Vector2d v2_n = V.row(F(i,1)); Eigen::Vector2d v3_n = V.row(F(i,2));
+
+    Eigen::MatrixXd T2_Homo(3,3);
+    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);
+    if (det < 0) {
+      flip_idx.push_back(i);
+    }
+  }
+
+  Eigen::VectorXi ret(flip_idx.size());
+  for (unsigned i=0; i<flip_idx.size();++i)
+    ret[i] = flip_idx[i];
+
+  return ret;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+#endif

include/igl/flipped_triangles_ids.h

@@ -0,0 +1,38 @@
+// 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_IDS_H
+#define IGL_FLIPPED_TRIANGLES_IDS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <vector>
+namespace igl
+{
+  // Finds the ids of the flipped triangles of the mesh V,F in the UV mapping uv
+  // Templates:
+  //   Scalar  should be a floating point number type
+  //   Index   should be an integer type
+  // Inputs:
+  //   V       #V by 2 list of mesh vertex positions
+  //   F       #F by dim list of mesh faces (must be triangles)
+  // Outputs:
+  //   A vector containing the ids of the flipped triangles
+  template <typename Scalar, typename Index>
+  IGL_INLINE Eigen::VectorXi flipped_triangles_ids(
+    const Eigen::PlainObjectBase<Scalar> & V,
+    const Eigen::PlainObjectBase<Index> & F
+  );
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flipped_triangles_ids.cpp"
+#endif
+
+#endif

include/igl/loop.cpp

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

include/igl/loop.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Oded Stein <oded.stein@columbia.edu>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_LOOP_H
+#define IGL_LOOP_H
+
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+    // LOOP Given the triangle mesh [V, F], where n_verts = V.rows(), computes newV and a sparse matrix S s.t. [newV, newF] is the subdivided mesh where newV = S*V.
+    //
+    // Inputs:
+    //  n_verts an integer (number of mesh vertices)
+    //  F an m by 3 matrix of integers of triangle faces
+    // Outputs:
+    //  S a sparse matrix (will become the subdivision matrix)
+    //  newF a matrix containing the new faces
+    IGL_INLINE void loop(const int n_verts,
+                         const Eigen::MatrixXi& F,
+                         Eigen::SparseMatrix<double>& S,
+                         Eigen::MatrixXi& newF);
+    
+    // LOOP Given the triangle mesh [V, F], computes number_of_subdivs steps of loop subdivision and outputs the new mesh [newV, newF]
+    //
+    // Inputs:
+    //  V an n by 3 matrix of vertices
+    //  F an m by 3 matrix of integers of triangle faces
+    //  number_of_subdivs an integer that specifies how many subdivision steps to do
+    // Outputs:
+    //  newV a matrix containing the new vertices
+    //  newF a matrix containing the new faces
+    IGL_INLINE void loop(const Eigen::MatrixXd& V,
+                         const Eigen::MatrixXi& F,
+                         Eigen::MatrixXd& newV,
+                         Eigen::MatrixXi& newF,
+                         const int number_of_subdivs = 1);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "loop.cpp"
+#endif
+
+#endif

include/igl/lscm.cpp

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

include/igl/min_quad_with_fixed.cpp

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

include/igl/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/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/viewer/Viewer.cpp

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

include/igl/viewer/Viewer.h

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

index.html

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

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

shared/cmake/CMakeLists.txt

@@ -340,6 +340,9 @@ endif()
 
 ### Compile the png parts ###
 if(LIBIGL_WITH_PNG)
+  if(LIBIGL_WITH_NANOGUI)
+    set(STBI_LOAD OFF CACHE BOOL " " FORCE)
+  endif()
   set(STB_IMAGE_DIR "${LIBIGL_EXTERNAL}/stb_image")
   add_subdirectory("${STB_IMAGE_DIR}" "stb_image")
   list(APPEND LIBIGL_INCLUDE_DIRS ${STB_IMAGE_DIR})

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

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

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

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

tutorial/tutorial.html.REMOVED.git-id

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

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-81c8a0d29caab90b96cb6372b610a0ea4c48cb11
+d0e101aea39668e374ff776b29cdd79667b7d8bd