Robust isolines (#2251)

* robust isolines, test

* documentation

include/igl/cat.cpp

@@ -236,6 +236,10 @@ IGL_INLINE void igl::cat(const int dim, const std::vector<T> & A, Eigen::PlainOb
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::cat<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int, std::vector<Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::allocator<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::cat<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(int, std::vector<Eigen::Matrix<double, -1, -1, 0, -1, -1>, std::allocator<Eigen::Matrix<double, -1, -1, 0, -1, -1> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::cat<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(int, std::vector<Eigen::Matrix<float, -1, -1, 0, -1, -1>, std::allocator<Eigen::Matrix<float, -1, -1, 0, -1, -1> > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::cat<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(int, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&);

include/igl/edge_crossings.cpp

@@ -0,0 +1,43 @@
+#include "edge_crossings.h"
+#include "list_to_matrix.h"
+
+template <
+  typename DeriveduE,
+  typename DerivedS,
+  typename DerivedT>
+void igl::edge_crossings(
+  const Eigen::MatrixBase<DeriveduE> & uE,
+  const Eigen::MatrixBase<DerivedS> & S,
+  const typename DerivedS::Scalar val,
+  std::unordered_map<int,int> & uE2I,
+  Eigen::PlainObjectBase<DerivedT> & T)
+{
+  using Scalar = typename DerivedS::Scalar;
+  {
+    std::vector<Scalar> vT;
+    for(int u = 0;u<uE.rows();u++)
+    {
+      // Does edge cross isovalue?
+      const int i = uE(u,0);
+      const int j = uE(u,1);
+      const Scalar Si = S(i);
+      const Scalar Sj = S(j);
+      // Specifically do not include endpoints
+      if( (Si-val)*(Sj-val) >= 0)
+      {
+        continue;
+      }
+      // find crossing point
+      const Scalar t = (val-Si)/(Sj-Si);
+      uE2I[u] = vT.size();
+      vT.emplace_back(t);
+    }
+    igl::list_to_matrix(vT,T);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::edge_crossings<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::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>::Scalar, std::unordered_map<int, int, std::hash<int>, std::equal_to<int>, std::allocator<std::pair<int const, int> > >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/edge_crossings.h

@@ -0,0 +1,36 @@
+#ifndef IGL_EDGE_CROSSINGS_H
+#define IGL_EDGE_CROSSINGS_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <unordered_map>
+
+namespace igl
+{
+  /// Compute the each point that a scalar field crosses a specified value along
+  /// an edge of a mesh.
+  ///
+  /// @param[in] uE  #E by 2 list of edge indices
+  /// @param[in] S  #V by 1 list of scalar field values
+  /// @param[in] val  value to check for crossings
+  /// @param[out] uE2I #T map from edge index to index in T
+  /// @param[out] T  #T by 1 list of parametric coordinates of crossings
+  ///
+  /// \see isolines, isolines_intrinsic
+  template <
+    typename DeriveduE,
+    typename DerivedS,
+    typename DerivedT>
+  void edge_crossings(
+    const Eigen::MatrixBase<DeriveduE> & uE,
+    const Eigen::MatrixBase<DerivedS> & S,
+    const typename DerivedS::Scalar val,
+    std::unordered_map<int,int> & uE2I,
+    Eigen::PlainObjectBase<DerivedT> & T);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "edge_crossings.cpp"
+#endif
+
+#endif

include/igl/isolines.cpp

@@ -1,116 +1,49 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
-// Copyright (C) 2017 Oded Stein <[email protected]>
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
 //
 // 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 "isolines.h"
-
-#include <vector>
-#include <array>
-#include <iostream>
-
-#include "remove_duplicate_vertices.h"
-
-
-template <typename DerivedV,
-typename DerivedF,
-typename DerivedZ,
-typename DerivedIsoV,
-typename DerivedIsoE>
-IGL_INLINE void igl::isolines(
-                              const Eigen::MatrixBase<DerivedV>& V,
-                              const Eigen::MatrixBase<DerivedF>& F,
-                              const Eigen::MatrixBase<DerivedZ>& z,
-                              const int n,
-                              Eigen::PlainObjectBase<DerivedIsoV>& isoV,
-                              Eigen::PlainObjectBase<DerivedIsoE>& isoE)
+#include "isolines_intrinsic.h"
+
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedS,
+  typename Derivedvals,
+  typename DerivediV,
+  typename DerivediE,
+  typename DerivedI>
+void igl::isolines(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedS> & S,
+  const Eigen::MatrixBase<Derivedvals> & vals,
+  Eigen::PlainObjectBase<DerivediV> & iV,
+  Eigen::PlainObjectBase<DerivediE> & iE,
+  Eigen::PlainObjectBase<DerivedI> & I)
 {
-    //Constants
-    const int dim = V.cols();
-    assert(dim==2 || dim==3);
-    const int nVerts = V.rows();
-    assert(z.rows() == nVerts &&
-           "There must be as many function entries as vertices");
-    const int nFaces = F.rows();
-    const int np1 = n+1;
-    const double min = z.minCoeff(), max = z.maxCoeff();
-
-
-    //Following http://www.alecjacobson.com/weblog/?p=2529
-    typedef typename DerivedZ::Scalar Scalar;
-    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vec;
-    Vec iso(np1);
-    for(int i=0; i<np1; ++i)
-        iso(i) = Scalar(i)/Scalar(n)*(max-min) + min;
-
-    typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix;
-    std::array<Matrix,3> t{{Matrix(nFaces, np1),
-        Matrix(nFaces, np1), Matrix(nFaces, np1)}};
-    for(int i=0; i<nFaces; ++i) {
-        for(int k=0; k<3; ++k) {
-            const Scalar z1=z(F(i,k), 0), z2=z(F(i,(k+1)%3), 0);
-            for(int j=0; j<np1; ++j) {
-                t[k](i,j) = (iso(j)-z1) / (z2-z1);
-                if(t[k](i,j)<0 || t[k](i,j)>1)
-                    t[k](i,j) = std::numeric_limits<Scalar>::quiet_NaN();
-            }
-        }
-    }
-
-    std::array<std::vector<int>,3> Fij, Iij;
-    for(int i=0; i<nFaces; ++i) {
-        for(int j=0; j<np1; ++j) {
-            for(int k=0; k<3; ++k) {
-                const int kp1=(k+1)%3, kp2=(k+2)%3;
-                if(std::isfinite(t[kp1](i,j)) && std::isfinite(t[kp2](i,j))) {
-                    Fij[k].push_back(i);
-                    Iij[k].push_back(j);
-                }
-            }
-        }
-    }
-
-    const int K = Fij[0].size()+Fij[1].size()+Fij[2].size();
-    isoV.resize(2*K, dim);
-    int b = 0;
-    for(int k=0; k<3; ++k) {
-        const int kp1=(k+1)%3, kp2=(k+2)%3;
-        for(int i=0; i<Fij[k].size(); ++i) {
-            isoV.row(b+i) = (1.-t[kp1](Fij[k][i],Iij[k][i]))*
-            V.row(F(Fij[k][i],kp1)) +
-            t[kp1](Fij[k][i],Iij[k][i])*V.row(F(Fij[k][i],kp2));
-            isoV.row(K+b+i) = (1.-t[kp2](Fij[k][i],Iij[k][i]))*
-            V.row(F(Fij[k][i],kp2)) +
-            t[kp2](Fij[k][i],Iij[k][i])*V.row(F(Fij[k][i],k));
-        }
-        b += Fij[k].size();
-    }
-
-    isoE.resize(K,2);
-    for(int i=0; i<K; ++i)
-        isoE.row(i) << i, K+i;
-
-
-    //Remove double entries
-    typedef typename DerivedIsoV::Scalar LScalar;
-    typedef typename DerivedIsoE::Scalar LInt;
-    typedef Eigen::Matrix<LInt, Eigen::Dynamic, 1> LIVec;
-    typedef Eigen::Matrix<LScalar, Eigen::Dynamic, Eigen::Dynamic> LMat;
-    typedef Eigen::Matrix<LInt, Eigen::Dynamic, Eigen::Dynamic> LIMat;
-    LIVec dummy1, dummy2;
-    igl::remove_duplicate_vertices(LMat(isoV), LIMat(isoE),
-                                   2.2204e-15, isoV, dummy1, dummy2, isoE);
-
+  using Scalar = typename DerivedS::Scalar;
+  using MatrixXS = Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>;
+
+  MatrixXS iB;
+  Eigen::VectorXi iFI;
+  isolines_intrinsic(F,S,vals,iB,iFI,iE,I);
+  iV.resize(iB.rows(),V.cols());
+  for(int i = 0;i<iB.rows();i++)
+  {
+    iV.row(i) = 
+      iB(i,0)*V.row(F(iFI(i),0))+
+      iB(i,1)*V.row(F(iFI(i),1))+
+      iB(i,2)*V.row(F(iFI(i),2));
+  }
 }
 
-
-
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::isolines<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::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int const, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > &);
+// generated by autoexplicit.sh
+template void igl::isolines<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::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<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<int, -1, 1, 0, -1, 1> >&);
 #endif
-

include/igl/isolines.h

@@ -1,6 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
-// Copyright (C) 2017 Oded Stein <[email protected]>
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
 //
 // 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
@@ -11,33 +11,47 @@
 #define IGL_ISOLINES_H
 #include "igl_inline.h"
 
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
+#include <Eigen/Core>
 
 
 namespace igl
 {
-  /// Constructs isolines for a function z given on a mesh (V,F)
+  /// Compute isolines of a scalar field on a triangle mesh.
+  ///
+  /// Isolines may cross perfectly at vertices. The output should not contain
+  /// degenerate segments (so long as the input does not contain degenerate
+  /// faces). The output segments are *oriented* so that isolines curl
+  /// counter-clockwise around local maxima (i.e., for 2D scalar fields). Unless
+  /// an isoline hits a boundary, it should be a closed loop. Isolines may run
+  /// perfectly along boundaries. Isolines should appear just "above" constants
+  /// regions.
   ///
   /// @param[in] V  #V by dim list of mesh vertex positions
-  /// @param[in] F  #F by 3 list of mesh faces (must be triangles)
-  /// @param[in] z  #V by 1 list of function values evaluated at vertices
-  /// @param[in] n  the number of desired isolines
-  /// @param[out] isoV  #isoV by dim list of isoline vertex positions
-  /// @param[out] isoE  #isoE by 2 list of isoline edge positions
+  /// @param[in] F  #F by 3 list of mesh triangle indices into V
+  /// @param[in] S  #S by 1 list of per-vertex scalar values
+  /// @param[in] vals  #vals by 1 list of values to compute isolines for
+  /// @param[out] iV  #iV by dim list of isoline vertex positions
+  /// @param[out] iE  #iE by 2 list of edge indices into iV
+  /// @param[out] I  #iE by 1 list of indices into vals indicating which value
+  ///   each segment belongs to
   ///
-  template <typename DerivedV,
-  typename DerivedF,
-  typename DerivedZ,
-  typename DerivedIsoV,
-  typename DerivedIsoE>
-  IGL_INLINE void isolines(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedF>& F,
-    const Eigen::MatrixBase<DerivedZ>& z,
-    const int n,
-    Eigen::PlainObjectBase<DerivedIsoV>& isoV,
-    Eigen::PlainObjectBase<DerivedIsoE>& isoE);
+  /// \see isolines_intrinsic, edge_crossings
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedS,
+    typename Derivedvals,
+    typename DerivediV,
+    typename DerivediE,
+    typename DerivedI>
+  void isolines(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedS> & S,
+    const Eigen::MatrixBase<Derivedvals> & vals,
+    Eigen::PlainObjectBase<DerivediV> & iV,
+    Eigen::PlainObjectBase<DerivediE> & iE,
+    Eigen::PlainObjectBase<DerivedI> & I);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/isolines_intrinsic.cpp

@@ -0,0 +1,300 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
+//
+// 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 "isolines_intrinsic.h"
+#include "edge_crossings.h"
+#include "cat.h"
+#include "unique_edge_map.h"
+#ifndef NDEBUG
+#  include "is_edge_manifold.h"
+#  include "is_vertex_manifold.h"
+#endif
+#include <unordered_map>
+#include <vector>
+
+template <
+  typename DerivedF,
+  typename DerivedS,
+  typename Derivedvals,
+  typename DerivediB,
+  typename DerivediFI,
+  typename DerivediE,
+  typename DerivedI>
+void igl::isolines_intrinsic(
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedS> & S,
+  const Eigen::MatrixBase<Derivedvals> & vals,
+  Eigen::PlainObjectBase<DerivediB> & iB,
+  Eigen::PlainObjectBase<DerivediFI> & iFI,
+  Eigen::PlainObjectBase<DerivediE> & iE,
+  Eigen::PlainObjectBase<DerivedI> & I)
+{
+  using Scalar = typename DerivedS::Scalar;
+
+  Eigen::MatrixXi uE;
+  Eigen::VectorXi EMAP,uEC,uEE;
+
+  {
+    Eigen::MatrixXi E;
+    igl::unique_edge_map(F,E,uE,EMAP,uEC,uEE);
+  }
+
+  {
+    std::vector<DerivediB> viB(vals.size());
+    std::vector<DerivediFI> viFI(vals.size());
+    std::vector<DerivediE> viE(vals.size());
+    std::vector<DerivedI> vI(vals.size());
+    int num_vertices = 0;
+    for(int j = 0;j<vals.size();j++)
+    {
+      isolines_intrinsic(F,S,uE,EMAP,uEC,uEE,vals[j],viB[j],viFI[j],viE[j]);
+      viE[j].array() += num_vertices;
+      num_vertices += viB[j].rows();
+      vI[j] = Eigen::VectorXi::Constant(viE[j].rows(),j);
+    }
+    igl::cat(1,viB,iB);
+    igl::cat(1,viFI,iFI);
+    igl::cat(1,viE,iE);
+    igl::cat(1,vI,I);
+  }
+}
+
+template <
+  typename DerivedF,
+  typename DerivedS,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename DeriveduEC,
+  typename DeriveduEE,
+  typename DerivediB,
+  typename DerivediFI,
+  typename DerivediE>
+void igl::isolines_intrinsic(
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedS> & S,
+  const Eigen::MatrixBase<DeriveduE> & uE,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  const Eigen::MatrixBase<DeriveduEC> & uEC,
+  const Eigen::MatrixBase<DeriveduEE> & uEE,
+  const typename DerivedS::Scalar val,
+  Eigen::PlainObjectBase<DerivediB> & iB,
+  Eigen::PlainObjectBase<DerivediFI> & iFI,
+  Eigen::PlainObjectBase<DerivediE> & iE)
+{
+  using Scalar = typename DerivedS::Scalar;
+  using VectorXS = Eigen::Matrix<Scalar,Eigen::Dynamic,1>;
+
+  std::unordered_map<int,int> uE2I;
+  Eigen::Matrix<Scalar,Eigen::Dynamic,1> T;
+  igl::edge_crossings(uE,S,val,uE2I,T);
+
+  iB.resize(uE2I.size(),F.cols());
+  iFI.resize(uE2I.size());
+  Eigen::VectorXi U(uE2I.size());
+  const auto flipped = 
+    [&uE,&F,&EMAP,&uEE,&uEC](const int u, const int f, const int k)->bool
+  {
+    return uE(u,0) != F(f,(k+1)%3);
+  };
+  // for each value in uE2I
+  {
+    for(auto & pair : uE2I)
+    {
+      const int u = pair.first;
+      const int w = pair.second;
+      // first face incident on uE(u,:)
+      const int e = uEE(uEC(u));
+      const int f = e % F.rows();
+      const int k = e / F.rows();
+      const bool flip = uE(u,0) != F(f,(k+1)%3);
+      const double t = T(w);
+      iB(w,k) = 0;
+      iB(w,(k+1)%3) = flip?  t:1-t;
+      iB(w,(k+2)%3) = flip?1-t:t;
+      iFI(w) = f;
+      U(w) = u;
+    }
+  }
+  
+  const int num_edge_crossings = iB.rows();
+  // Vertex crossings
+  std::unordered_map<int,int> V2I;
+  {
+    const auto add_vertex_crossing = [&iB,&iFI](const int k, const int i, const int j)
+    {
+      if(k >= iB.rows())
+      {
+        iB.conservativeResize(2*iB.rows()+1,Eigen::NoChange);
+        iFI.conservativeResize(2*iFI.rows()+1,Eigen::NoChange);
+      }
+      iFI(k) = i;
+      iB.row(k) << 0,0,0;
+      iB(k,j) = 1;
+    };
+    int k = iB.rows();
+    for(int i = 0;i<F.rows();i++)
+    {
+      for(int j = 0;j<3;j++)
+      {
+        const int v = F(i,j);
+        if(S(v) == val)
+        {
+          if(V2I.find(v) == V2I.end())
+          {
+            V2I[v] = k;
+            add_vertex_crossing(k++,i,j);
+          }
+        }
+      }
+    }
+    iB.conservativeResize(k,Eigen::NoChange);
+    iFI.conservativeResize(k,Eigen::NoChange);
+  }
+  const int num_vertex_crossings = iB.rows()-num_edge_crossings;
+
+  iE.resize(uE2I.size(),2);
+  const auto set_row = [&iE](const int k, const int i, const int j)
+  {
+    if(k >= iE.rows())
+    {
+      iE.conservativeResize(2*iE.rows()+1,Eigen::NoChange);
+    }
+    iE.row(k) << i,j;
+  };
+  {
+    int r = 0;
+    for(int f = 0;f < F.rows();f++)
+    {
+      // find first crossing edge
+      int i;
+      for(i = 0;i<3;i++)
+      {
+        if(uE2I.find(EMAP(f+F.rows()*i)) != uE2I.end())
+        {
+          break;
+        }
+      }
+      int j;
+      for(j = i+1;j<3;j++)
+      {
+        if(uE2I.find(EMAP(f+F.rows()*j)) != uE2I.end())
+        {
+          break;
+        }
+      }
+      if(j<3)
+      {
+        // Connect two edge crossings.
+
+        // other vertex
+        const int k = 3-i-j;
+        const int wi = uE2I[EMAP(f+F.rows()*i)];
+        const int wj = uE2I[EMAP(f+F.rows()*j)];
+        // flip orientation based on triangle gradient
+        bool flip = S(F(f,k)) < val;
+        flip = k%2? !flip:flip;
+        if(flip)
+        {
+          set_row(r++,wi,wj);
+        }else
+        {
+          set_row(r++,wj,wi);
+        }
+      }else if(i<3)
+      {
+        // The only valid vertex crossing is the opposite vertex
+        const int v = F(f,i);
+        // Is it a crossing?
+        assert(V2I.find(v) != V2I.end());
+        //if(V2I.find(v) != V2I.end())
+        {
+          const int wv = V2I[v];
+          const int wi = uE2I[EMAP(f+F.rows()*i)];
+          const bool flip = S(F(f,(i+1)%3)) > val;
+          if(flip)
+          {
+            set_row(r++,wi,wv);
+          }else
+          {
+            set_row(r++,wv,wi);
+          }
+        }
+      }else
+      {
+        // Could have 2 vertex crossings. We're only interested if there're exactly two and if the other vertex is "above".
+        int i = 0;
+        for(i = 0;i<3;i++)
+        {
+          if(S(F(f,i)) == val)
+          {
+            break;
+          }
+        }
+        int j;
+        for(j = i+1;j<3;j++)
+        {
+          if(S(F(f,j)) == val)
+          {
+            break;
+          }
+        }
+        if(j<3)
+        {
+          // check if the third is a crossing.
+          const int k = 3-i-j;
+          // Triangle is constant on the val. Skip.
+          if(S(F(f,k)) == val){ continue; }
+          // Is this a boundary edge?
+          const int u = EMAP(f+F.rows()*k);
+          const int count = uEC(u+1)-uEC(u);
+          if( count == 1 || S(F(f,k)) > val)
+          {
+            const int wi = V2I[F(f,i)];
+            const int wj = V2I[F(f,j)];
+            bool flip = S(F(f,k)) < val;
+            flip = k%2 ? !flip:flip;
+            if(flip)
+            {
+              set_row(r++,wj,wi);
+            }else
+            {
+              set_row(r++,wi,wj);
+            }
+          }
+        }
+      }
+    }
+    iE.conservativeResize(r,Eigen::NoChange);
+  }
+
+#ifndef NDEBUG
+  if(igl::is_vertex_manifold(F) && igl::is_edge_manifold(F))
+  {
+    // Check that every vertex has one in one out
+    Eigen::VectorXi in_count = Eigen::VectorXi::Zero(iB.rows());
+    Eigen::VectorXi out_count = Eigen::VectorXi::Zero(iB.rows());
+    for(int e = 0;e<iE.rows();e++)
+    {
+      const int i = iE(e,0);
+      out_count(i)++;
+      const int j = iE(e,1);
+      in_count(j)++;
+    }
+    for(int i = 0;i<iB.rows();i++)
+    {
+      assert(in_count(i) <= 1);
+      assert(out_count(i) <= 1);
+    }
+  }
+#endif
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::isolines_intrinsic<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::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::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<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<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/isolines_intrinsic.h

@@ -0,0 +1,97 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
+//
+// 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_ISOLINES_INTRINSIC_H
+#define IGL_ISOLINES_INTRINSIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Core>
+
+
+namespace igl
+{
+  /// Compute isolines of a scalar field on a triangle mesh intrinsically.
+  ///
+  /// See isolines.h for details.
+  ///
+  /// @param[in] F  #F by 3 list of mesh triangle indices into some V
+  /// @param[in] S  #S by 1 list of per-vertex scalar values
+  /// @param[in] vals  #vals by 1 list of values to compute isolines for
+  /// @param[out] iB  #iB by 3 list of barycentric coordinates so that 
+  ///   iV.row(i) = iB(i,0)*V.row(F(iFI(i,0)) +
+  ///               iB(i,1)*V.row(F(iFI(i,1)) +
+  ///               iB(i,2)*V.row(F(iFI(i,2))
+  /// @param[out] iF  #iB list of triangle indices for each row of iB (all
+  ///   points will either lie on an edge or vertex: an arbitrary incident face
+  ///   will be given).
+  /// @param[out] iE  #iE by 2 list of edge indices into iB
+  /// @param[out] I  #iE by 1 list of indices into vals indicating which value
+  ///   each segment belongs to
+  ///
+  /// \see isolines, edge_crossings
+  template <
+    typename DerivedF,
+    typename DerivedS,
+    typename Derivedvals,
+    typename DerivediB,
+    typename DerivediFI,
+    typename DerivediE,
+    typename DerivedI>
+  void isolines_intrinsic(
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedS> & S,
+    const Eigen::MatrixBase<Derivedvals> & vals,
+    Eigen::PlainObjectBase<DerivediB> & iB,
+    Eigen::PlainObjectBase<DerivediFI> & iFI,
+    Eigen::PlainObjectBase<DerivediE> & iE,
+    Eigen::PlainObjectBase<DerivedI> & I);
+  /// \overload
+  ///
+  /// @param[in] val  scalar value to compute isoline at
+  /// @param[in] uE  #uE by 2 list of unique undirected edges
+  /// @param[in] EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  ///     undirected edge so that uE(EMAP(f+#F*c)) is the unique edge
+  ///     corresponding to E.row(f+#F*c)
+  /// @param[in] uEC  #uE+1 list of cumulative counts of directed edges sharing each
+  ///     unique edge so the uEC(i+1)-uEC(i) is the number of directed edges
+  ///     sharing the ith unique edge.
+  /// @param[in] uEE  #E list of indices into E, so that the consecutive segment of
+  ///     indices uEE.segment(uEC(i),uEC(i+1)-uEC(i)) lists all directed edges
+  ///     sharing the ith unique edge.
+  ///
+  /// \see unique_edge_map
+  template <
+    typename DerivedF,
+    typename DerivedS,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename DeriveduEC,
+    typename DeriveduEE,
+    typename DerivediB,
+    typename DerivediFI,
+    typename DerivediE>
+  void isolines_intrinsic(
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedS> & S,
+    const Eigen::MatrixBase<DeriveduE> & uE,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+    const Eigen::MatrixBase<DeriveduEC> & uEC,
+    const Eigen::MatrixBase<DeriveduEE> & uEE,
+    const typename DerivedS::Scalar val,
+    Eigen::PlainObjectBase<DerivediB> & iB,
+    Eigen::PlainObjectBase<DerivediFI> & iFI,
+    Eigen::PlainObjectBase<DerivediE> & iE);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "isolines_intrinsic.cpp"
+#endif
+
+#endif
+

tests/include/igl/isolines.cpp

@@ -0,0 +1,185 @@
+#include <test_common.h>
+#include <igl/isolines.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+TEST_CASE("isolines: broken-sphere", "[igl]" )
+{
+  // Test a case with edge-edge, edge-vertex, vertex-vertex, boundary edge, and
+  // degenerate boundaries
+  Eigen::MatrixXd V(42,3);
+  V<<
+        0,0.52573111211913359,0.85065080835203999,
+    0,-0.52573111211913359,0.85065080835203999,
+    0,0.52573111211913359,-0.85065080835203999,
+    0,-0.52573111211913359,-0.85065080835203999,
+    0.52573111211913359,0.85065080835203999,0,
+    -0.52573111211913359,0.85065080835203999,0,
+    0.52573111211913359,-0.85065080835203999,0,
+    -0.52573111211913359,-0.85065080835203999,0,
+    0.85065080835203999,0,0.52573111211913359,
+    -0.85065080835203999,0,0.52573111211913359,
+    0.85065080835203999,0,-0.52573111211913359,
+    -0.85065080835203999,0,-0.52573111211913359,
+    0,0,1,
+    0.3090169943749474,0.80901699437494734,0.5,
+    -0.3090169943749474,0.80901699437494734,0.5,
+    0.5,0.3090169943749474,0.80901699437494734,
+    -0.5,0.3090169943749474,0.80901699437494734,
+    0.3090169943749474,-0.80901699437494734,0.5,
+    -0.3090169943749474,-0.80901699437494734,0.5,
+    0.5,-0.3090169943749474,0.80901699437494734,
+    -0.5,-0.3090169943749474,0.80901699437494734,
+    0,0,-1,
+    0.3090169943749474,0.80901699437494734,-0.5,
+    -0.3090169943749474,0.80901699437494734,-0.5,
+    0.5,0.3090169943749474,-0.80901699437494734,
+    -0.5,0.3090169943749474,-0.80901699437494734,
+    0.3090169943749474,-0.80901699437494734,-0.5,
+    -0.3090169943749474,-0.80901699437494734,-0.5,
+    0.5,-0.3090169943749474,-0.80901699437494734,
+    -0.5,-0.3090169943749474,-0.80901699437494734,
+    0,1,0,
+    0.80901699437494734,0.5,0.3090169943749474,
+    0.80901699437494734,0.5,-0.3090169943749474,
+    -0.80901699437494734,0.5,0.3090169943749474,
+    -0.80901699437494734,0.5,-0.3090169943749474,
+    0,-1,0,
+    0.80901699437494734,-0.5,0.3090169943749474,
+    0.80901699437494734,-0.5,-0.3090169943749474,
+    -0.80901699437494734,-0.5,0.3090169943749474,
+    -0.80901699437494734,-0.5,-0.3090169943749474,
+    1,0,0,
+    -1,0,0;
+  Eigen::MatrixXi F(58,3);
+  F <<
+    0,12,15,
+    0,15,13,
+    0,13,14,
+    0,14,16,
+    0,16,12,
+    8,40,31,
+    4,32,22,
+    4,22,30,
+    5,23,34,
+    5,34,33,
+    1,18,17,
+    3,26,27,
+    8,31,15,
+    4,30,13,
+    5,33,14,
+    9,20,16,
+    6,36,17,
+    6,37,36,
+    10,32,40,
+    10,24,32,
+    2,23,22,
+    2,25,23,
+    11,41,34,
+    7,18,38,
+    7,35,18,
+    6,35,26,
+    7,39,27,
+    11,25,29,
+    2,24,21,
+    8,15,19,
+    4,13,31,
+    5,14,30,
+    9,16,33,
+    4,31,32,
+    2,22,24,
+    5,30,23,
+    11,34,25,
+    9,33,41,
+    7,38,39,
+    6,17,35,
+    7,27,35,
+    2,21,25,
+    10,28,24,
+    6,26,37,
+    19,15,12,
+    31,13,15,
+    30,14,13,
+    33,16,14,
+    20,12,16,
+    32,31,40,
+    24,22,32,
+    23,30,22,
+    25,34,23,
+    41,33,34,
+    35,17,18,
+    35,27,26,
+    25,21,29,
+    24,28,21;
+  Eigen::MatrixXd iV;
+  Eigen::MatrixXi iE;
+  Eigen::VectorXi I;
+  Eigen::VectorXd vals(3);
+  vals<<-0.5,0,0.5;
+  igl::isolines(V,F,V.col(1).eval(),vals,iV,iE,I);
+  Eigen::MatrixXd iV_gt(30,3);
+  iV_gt<<
+      0.809016994374947,-0.5, 0.309016994374947,
+      0.809016994374947,-0.5,-0.309016994374947,
+     -0.809016994374947,-0.5, 0.309016994374947,
+     -0.809016994374947,-0.5,-0.309016994374947,
+                    0.5,   0, 0.809016994374947,
+                   -0.5,   0, 0.809016994374947,
+                    0.5,   0,-0.809016994374947,
+                   -0.5,   0,-0.809016994374947,
+                      0,   0,                 1,
+       0.85065080835204,   0, 0.525731112119134,
+                      1,   0,                 0,
+      -0.85065080835204,   0, 0.525731112119134,
+       0.85065080835204,   0,-0.525731112119134,
+      -0.85065080835204,   0,-0.525731112119134,
+                     -1,   0,                 0,
+                      0,   0,                -1,
+                      0, 0.5, 0.857960478079794,
+     0.0593664879495925, 0.5, 0.845707501720507,
+    -0.0593664879495925, 0.5, 0.845707501720507,
+                      0, 0.5,-0.857960478079794,
+     0.0593664879495925, 0.5,-0.845707501720507,
+    -0.0593664879495925, 0.5,-0.845707501720507,
+      0.427050983124842, 0.5, 0.690983005625053,
+     -0.427050983124842, 0.5, 0.690983005625052,
+      0.427050983124842, 0.5,-0.690983005625053,
+     -0.427050983124842, 0.5,-0.690983005625053,
+      0.809016994374947, 0.5, 0.309016994374947,
+      0.809016994374947, 0.5,-0.309016994374947,
+     -0.809016994374947, 0.5,-0.309016994374947,
+     -0.809016994374947, 0.5, 0.309016994374947;
+  Eigen::MatrixXi iE_gt(28,2);
+  iE_gt<<
+     0, 1,
+     3, 2,
+     9,10,
+    11, 5,
+    10,12,
+    13,14,
+     7,13,
+     4, 9,
+    14,11,
+    12, 6,
+     8, 4,
+     5, 8,
+    15, 7,
+     6,15,
+    16,17,
+    17,22,
+    23,18,
+    18,16,
+    28,29,
+    21,25,
+    20,19,
+    26,27,
+    24,20,
+    19,21,
+    22,26,
+    29,23,
+    27,24,
+    25,28;
+
+  test_common::assert_near(iV,iV_gt,1e-15);
+  test_common::assert_eq(iE,iE_gt);
+}

tutorial/905_Isolines/main.cpp

@@ -0,0 +1,51 @@
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/material_colors.h>
+#include <igl/isolines.h>
+
+int main(int argc, char *argv[])
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(
+      argc>1?argv[1]:TUTORIAL_SHARED_PATH "/fertility.off",V,F);
+  // Use y-coordinate as scalar field
+  Eigen::VectorXd f = V.col(1);
+
+  Eigen::MatrixXd iV;
+  Eigen::MatrixXi iE;
+  {
+    // How many isolines in the range (min,max)?
+    const int n = argc>2?atoi(argv[2]):128;
+    // This is actually unnecessary since isolines will not output degenerate
+    // segments.
+    //Eigen::VectorXd vals = Eigen::VectorXd::LinSpaced(n+2,f.minCoeff(),f.maxCoeff()).segment(1,n);
+    // Instead just use all n+2 and if the min-,max-value isolines are
+    // non-degneerate then we'll compute them, too.
+    Eigen::VectorXd vals = Eigen::VectorXd::LinSpaced(n+2,f.minCoeff(),f.maxCoeff());
+    {
+      Eigen::VectorXi I;
+      igl::isolines(V,F,f,vals,iV,iE,I);
+    }
+  }
+
+  // Plot the mesh
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V, F);
+  viewer.data().label_size = 10;
+  viewer.data().set_face_based(true);
+  viewer.data().show_faces = true;
+  viewer.data().show_lines = false;
+  viewer.data().uniform_colors(
+    Eigen::Vector3d(0.94*viewer.core().background_color.head<3>().cast<double>()),
+    Eigen::Vector3d(0.05*viewer.core().background_color.head<3>().cast<double>()),
+    Eigen::Vector3d(0.01*viewer.core().background_color.head<3>().cast<double>()));
+
+  viewer.core().lighting_factor = 0.5;
+  viewer.data().set_edges(iV, iE, 
+      Eigen::RowVector3d(igl::GOLD_DIFFUSE[0], igl::GOLD_DIFFUSE[1], igl::GOLD_DIFFUSE[2]));
+  viewer.data().line_width = 1;
+  
+  viewer.launch();
+}
+

tutorial/CMakeLists.txt

@@ -126,4 +126,5 @@ if(LIBIGL_TUTORIALS_CHAPTER9)
     igl_add_tutorial(902_VectorParallelTransport igl::glfw)
     igl_add_tutorial(903_FastFindSelfIntersections igl::imgui igl::glfw)
     igl_add_tutorial(904_FastFindIntersections igl::imgui igl::glfw)
+    igl_add_tutorial(905_Isolines igl::imgui igl::glfw)
 endif()