* adding streamline code and tutorial

* adding python bindings for streamline and python tutorial


Former-commit-id: 9abb24e0bc9c9f548270846f874e738d80aa11cc

include/igl/segment_segment_intersect.cpp

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// 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
+// TODO: Explicit template specialization
+#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) 2013 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// 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,156 @@
+#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 <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,84 @@
+//
+// Created by Francisca Gil Ureta on 7/1/16.
+//
+
+#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

python/py_doc.cpp

@@ -988,6 +988,24 @@ 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).
   //

python/py_doc.h

@@ -81,6 +81,8 @@ 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;

python/py_igl.cpp

@@ -69,6 +69,7 @@
 #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>
@@ -147,6 +148,7 @@ 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"

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

@@ -125,6 +125,8 @@ PYBIND11_PLUGIN(pyigl) {
            slice_mask
            slice_tets
            sortrows
+           streamlines_init
+           streamlines_next
            triangle_triangle_adjacency
            triangle_triangulate
            unique

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

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,165 @@
+#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 data;
+igl::StreamlineState 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, data, 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(state.start_point, 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, data, 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 state0;
+    state0 = state;
+    igl::streamlines_next(V, F, data, state0);
+    Eigen::MatrixXd v = state0.end_point - state0.start_point;
+    v.rowwise().normalize();
+
+    viewer.data.add_edges(state0.start_point,
+                          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()