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

Former-commit-id: 77ecfa6a7b0d7174d7b5f56811588a7ef85eab61

.gitignore

@@ -96,3 +96,5 @@ python/build4
 python/scripts/generated
 python/builddebug
 tutorial/.idea
+python/buildstatic
+tutorial/cmake-build-debug

.travis.yml

@@ -14,7 +14,7 @@ matrix:
         - cmake -DCMAKE_CXX_COMPILER=g++-4.8 -DCMAKE_C_COMPILER=gcc-4.8 -DLIBIGL_WITH_EMBREE=OFF -DLIBIGL_USE_STATIC_LIBRARY=ON ../
         - make -j 2
         - cd ../tutorial
-        - python 101_FileIO.py
+        - python 101_FileIO.py || { cd ../; mkdir build2; cd build2; cmake -DCMAKE_CXX_COMPILER=g++-4.8 -DCMAKE_C_COMPILER=gcc-4.8 -DLIBIGL_WITH_EMBREE=OFF -DLIBIGL_USE_STATIC_LIBRARY=ON - DCHECK_UNDEFINED=ON; make -j 2; }
         - cd ../../
         - cd tutorial
         - mkdir build

README.md

@@ -107,13 +107,13 @@ For more information see our [tutorial](tutorial/tutorial.html).
 Libigl compartmentalizes its **optional** dependences via its directory
 organization in the `include/` folder. All header files located _directly_ in
 the `include/igl/` folder have only stl and Eigen as dependencies. For example,
-all of the headers that depend on CGAL are located in `include/igl/cgal`. For a
-full list of _optional_ dependencies check `optional/CMakeLists.txt`.
+all of the headers that depend on CGAL are located in `include/igl/copyleft/cgal`.
+For a full list of _optional_ dependencies check `optional/CMakeLists.txt`.
 
 ### GCC and the Optional CGAL Dependency
-The `include/igl/cgal/*.h` headers depend on CGAL. It has come to our attention
-that CGAL does not work properly with GCC 4.8. To the best of our knowledge,
-GCC 4.7 and clang will work correctly.
+The `include/igl/copyleft/cgal/*.h` headers depend on CGAL. It has come to
+our attention that CGAL does not work properly with GCC 4.8. To the best of
+our knowledge, GCC 4.7 and clang will work correctly.
 
 ### OpenMP and Windows
 Some of our functions will take advantage of OpenMP if available. However, it

include/igl/HalfEdgeIterator.h

@@ -100,9 +100,9 @@ namespace igl
     bool reverse;
 
     // All the same type? This is likely to break.
-    const DerivedF & F;
-    const DerivedFF & FF;
-    const DerivedFFi & FFi;
+    const Eigen::PlainObjectBase<DerivedF> & F;
+    const Eigen::PlainObjectBase<DerivedFF> & FF;
+    const Eigen::PlainObjectBase<DerivedFFi> & FFi;
   };
 
 }

include/igl/copyleft/marching_cubes.cpp

@@ -25,39 +25,42 @@
 
 #include "marching_cubes.h"
 #include "marching_cubes_tables.h"
-#include <map>
+
+#include <unordered_map>
+
 
 extern const int edgeTable[256];
 extern const int triTable[256][2][17];
 extern const int polyTable[8][16];
 
-class EdgeKey
+struct EdgeKey
 {
-public:
-
-  EdgeKey(unsigned i0, unsigned i1) {
-    if (i0 < i1)  { i0_ = i0;  i1_ = i1; }
-    else            { i0_ = i1;  i1_ = i0; }
-  }
+  EdgeKey(unsigned i0, unsigned i1) : i0_(i0), i1_(i1) {}
 
-  bool operator<(const EdgeKey& _rhs) const
+  bool operator==(const EdgeKey& _rhs) const
   {
-    if (i0_ != _rhs.i0_)
-      return (i0_ < _rhs.i0_);
-    else
-      return (i1_ < _rhs.i1_);
+    return i0_ == _rhs.i0_ && i1_ == _rhs.i1_;
   }
 
-private:
   unsigned i0_, i1_;
 };
 
+struct EdgeHash
+{
+    std::size_t operator()(const EdgeKey& key) const {
+        std::size_t seed = 0;
+        seed ^= key.i0_ + 0x9e3779b9 + (seed<<6) + (seed>>2); // Copied from boost::hash_combine
+        seed ^= key.i1_ + 0x9e3779b9 + (seed<<6) + (seed>>2);
+        return std::hash<std::size_t>()(seed);
+    }
+};
+
 
 template <typename Derivedvalues, typename Derivedpoints,typename Derivedvertices, typename DerivedF>
 class MarchingCubes
 {
-  typedef std::map<EdgeKey, unsigned>  MyMap;
-  typedef typename MyMap::const_iterator   MyMapIterator;
+  typedef std::unordered_map<EdgeKey, unsigned, EdgeHash> MyMap;
+  typedef typename MyMap::const_iterator                  MyMapIterator;
 
 public:
   MarchingCubes(

include/igl/cut_mesh.cpp

@@ -23,16 +23,16 @@ namespace igl {
   public:
     // Input
     //mesh
-    const DerivedV &V;
-    const DerivedF &F;
+    const Eigen::PlainObjectBase<DerivedV> &V;
+    const Eigen::PlainObjectBase<DerivedF> &F;
     // TT is the same type as TTi? This is likely to break at some point
-    const DerivedTT &TT;
-    const DerivedTT &TTi;
+    const Eigen::PlainObjectBase<DerivedTT> &TT;
+    const Eigen::PlainObjectBase<DerivedTT> &TTi;
     const std::vector<std::vector<VFType> >& VF;
     const std::vector<std::vector<VFType> >& VFi;
     const std::vector<bool> &V_border; // bool
     //edges to cut
-    const DerivedC &Handle_Seams; // 3 bool
+    const Eigen::PlainObjectBase<DerivedC> &Handle_Seams; // 3 bool
 
     // total number of scalar variables
     int num_scalar_variables;

include/igl/cut_mesh_from_singularities.cpp

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

include/igl/find_cross_field_singularities.cpp

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

include/igl/hessian.cpp

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//  and 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 "hessian.h"
+#include <vector>
+
+#include "grad.h"
+#include "igl/doublearea.h"
+#include "igl/repdiag.h"
+
+
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::hessian(
+                             const Eigen::MatrixBase<DerivedV> & V,
+                             const Eigen::MatrixBase<DerivedF> & F,
+                             Eigen::SparseMatrix<Scalar>& H)
+{
+    typedef typename DerivedV::Scalar denseScalar;
+    typedef typename Eigen::Matrix<denseScalar, Eigen::Dynamic, 1> VecXd;
+    typedef typename Eigen::SparseMatrix<Scalar> SparseMat;
+    typedef typename Eigen::DiagonalMatrix
+                       <Scalar, Eigen::Dynamic, Eigen::Dynamic> DiagMat;
+    
+    int dim = V.cols();
+    assert((dim==2 || dim==3) &&
+           "The dimension of the vertices should be 2 or 3");
+    
+    //Construct the combined gradient matric
+    SparseMat G;
+    igl::grad(Eigen::PlainObjectBase<DerivedV>(V),
+              Eigen::PlainObjectBase<DerivedF>(F),
+              G, false);
+    SparseMat GG(F.rows(), dim*V.rows());
+    GG.reserve(G.nonZeros());
+    for(int i=0; i<dim; ++i)
+        GG.middleCols(i*G.cols(),G.cols()) = G.middleRows(i*F.rows(),F.rows());
+    SparseMat D;
+    igl::repdiag(GG,dim,D);
+    
+    //Compute area matrix
+    VecXd areas;
+    igl::doublearea(V, F, areas);
+    DiagMat A = (0.5*areas).replicate(dim,1).asDiagonal();
+    
+    //Compute FEM Hessian
+    H = D.transpose()*A*G;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::hessian<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/hessian.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//  and 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_FEM_HESSIAN_H
+#define IGL_FEM_HESSIAN_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+    // Constructs the finite element Hessian matrix
+    // as described in https://arxiv.org/abs/1707.04348,
+    // Natural Boundary Conditions for Smoothing in Geometry Processing
+    // (Oded Stein, Eitan Grinspun, Max Wardetzky, Alec Jacobson)
+    // The interior vertices are NOT set to zero yet.
+    //
+    // Inputs:
+    //   V  #V by dim list of mesh vertex positions
+    //   F  #F by 3 list of mesh faces (must be triangles)
+    // Outputs:
+    //   H  #V by #V Hessian energy matrix, each column i
+    //     corresponding to V(i,:)
+    //
+    //
+    //
+    template <typename DerivedV, typename DerivedF, typename Scalar>
+    IGL_INLINE void hessian(
+                            const Eigen::MatrixBase<DerivedV> & V,
+                            const Eigen::MatrixBase<DerivedF> & F,
+                            Eigen::SparseMatrix<Scalar>& H);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "hessian.cpp"
+#endif
+
+#endif

include/igl/hessian_energy.cpp

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//   and 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 "hessian_energy.h"
+#include <vector>
+
+#include "hessian.h"
+#include "massmatrix.h"
+#include "boundary_loop.h"
+
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::hessian_energy(
+                                    const Eigen::MatrixBase<DerivedV> & V,
+                                    const Eigen::MatrixBase<DerivedF> & F,
+                                    Eigen::SparseMatrix<Scalar>& Q)
+{
+    typedef typename DerivedV::Scalar denseScalar;
+    typedef typename Eigen::Matrix<denseScalar, Eigen::Dynamic, 1> VecXd;
+    typedef typename Eigen::SparseMatrix<Scalar> SparseMat;
+    typedef typename Eigen::DiagonalMatrix
+                       <Scalar, Eigen::Dynamic, Eigen::Dynamic> DiagMat;
+    
+    int dim = V.cols();
+    assert((dim==2 || dim==3) &&
+           "The dimension of the vertices should be 2 or 3");
+    
+    SparseMat M;
+    igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
+    
+    //Kill non-interior DOFs
+    VecXd Mint = M.diagonal();
+    std::vector<std::vector<int> > bdryLoop;
+    igl::boundary_loop(Eigen::PlainObjectBase<DerivedF>(F),bdryLoop);
+    for(const std::vector<int>& loop : bdryLoop)
+        for(const int& bdryVert : loop)
+            Mint(bdryVert) = 0.;
+    
+    //Invert Mint
+    for(int i=0; i<Mint.rows(); ++i)
+        if(Mint(i) > 0)
+            Mint(i) = 1./Mint(i);
+    
+    //Repeat Mint to form diaginal matrix
+    DiagMat stackedMinv = Mint.replicate(dim*dim,1).asDiagonal();
+    
+    //Compute squared Hessian
+    SparseMat H;
+    igl::hessian(V,F,H);
+    Q = H.transpose()*stackedMinv*H;
+    
+}
+
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::hessian_energy<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+#endif

include/igl/hessian_energy.h

@@ -0,0 +1,45 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Alec Jacobson <alecjacobson@gmail.com>
+//  and 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_HESSIAN_ENERGY_H
+#define IGL_HESSIAN_ENERGY_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+    // Constructs the Hessian energy matrix using mixed FEM
+    // as described in https://arxiv.org/abs/1707.04348
+    // Natural Boundary Conditions for Smoothing in Geometry Processing
+    // (Oded Stein, Eitan Grinspun, Max Wardetzky, Alec Jacobson)
+    //
+    // Inputs:
+    //   V  #V by dim list of mesh vertex positions
+    //   F  #F by 3 list of mesh faces (must be triangles)
+    // Outputs:
+    //   Q  #V by #V Hessian energy matrix, each row/column i
+    //     corresponding to V(i,:)
+    //
+    //
+    //
+    template <typename DerivedV, typename DerivedF, typename Scalar>
+    IGL_INLINE void hessian_energy(
+                                   const Eigen::MatrixBase<DerivedV> & V,
+                                   const Eigen::MatrixBase<DerivedF> & F,
+                                   Eigen::SparseMatrix<Scalar>& Q);
+    
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "hessian_energy.cpp"
+#endif
+
+#endif

include/igl/isolines.cpp

@@ -0,0 +1,116 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 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 "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)
+{
+    //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)), z2=z(F(i,(k+1)%3));
+            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);
+    
+}
+
+
+
+#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> > &);
+#endif
+

include/igl/isolines.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 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_ISOLINES_H
+#define IGL_ISOLINES_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+
+namespace igl
+{
+    // Constructs isolines for a function z given on a mesh (V,F)
+    //
+    //
+    // Inputs:
+    //   V  #V by dim list of mesh vertex positions
+    //   F  #F by 3 list of mesh faces (must be triangles)
+    //   z  #V by 1 list of function values evaluated at vertices
+    //   n  the number of desired isolines
+    // Outputs:
+    //   isoV  #isoV by dim list of isoline vertex positions
+    //   isoE  #isoE by 2 list of isoline edge positions
+    //
+    //
+    
+    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);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "isolines.cpp"
+#endif
+
+#endif

include/igl/signed_angle.cpp

@@ -61,6 +61,8 @@ template Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, fal
 template Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>::Scalar igl::signed_angle<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&);
 // generated by autoexplicit.sh
 template Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>::Scalar igl::signed_angle<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&);
+template Eigen::Block<Eigen::Matrix<float, -1, 3, 0, -1, 3> const, 1, 3, false>::Scalar igl::signed_angle<Eigen::Block<Eigen::Matrix<float, -1, 3, 0, -1, 3> const, 1, 3, false>, Eigen::Block<Eigen::Matrix<float, -1, 3, 0, -1, 3> const, 1, 3, false>, Eigen::Matrix<float, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<float, -1, 3, 0, -1, 3> const, 1, 3, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<float, -1, 3, 0, -1, 3> const, 1, 3, false> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 2, 1, 1, 2> > const&);
+template Eigen::Block<Eigen::Matrix<float, -1, 3, 1, -1, 3> const, 1, 3, true>::Scalar igl::signed_angle<Eigen::Block<Eigen::Matrix<float, -1, 3, 1, -1, 3> const, 1, 3, true>, Eigen::Block<Eigen::Matrix<float, -1, 3, 1, -1, 3> const, 1, 3, true>, Eigen::Matrix<float, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<float, -1, 3, 1, -1, 3> const, 1, 3, true> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<float, -1, 3, 1, -1, 3> const, 1, 3, true> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 2, 1, 1, 2> > const&);
 #ifdef WIN32
 template float igl::signed_angle<class Eigen::Block<class Eigen::Matrix<float,-1,3,1,-1,3> const ,1,3,1>,class Eigen::Block<class Eigen::Matrix<float,-1,3,1,-1,3> const ,1,3,1>,class Eigen::Matrix<float,1,2,1,1,2> >(class Eigen::MatrixBase<class Eigen::Block<class Eigen::Matrix<float,-1,3,1,-1,3> const ,1,3,1> > const &,class Eigen::MatrixBase<class Eigen::Block<class Eigen::Matrix<float,-1,3,1,-1,3> const ,1,3,1> > const &,class Eigen::MatrixBase<class Eigen::Matrix<float,1,2,1,1,2> > const &);
 template float igl::signed_angle<class Eigen::Block<class Eigen::Matrix<float,-1,3,0,-1,3> const ,1,3,0>,class Eigen::Block<class Eigen::Matrix<float,-1,3,0,-1,3> const ,1,3,0>,class Eigen::Matrix<float,1,2,1,1,2> >(class Eigen::MatrixBase<class Eigen::Block<class Eigen::Matrix<float,-1,3,0,-1,3> const ,1,3,0> > const &,class Eigen::MatrixBase<class Eigen::Block<class Eigen::Matrix<float,-1,3,0,-1,3> const ,1,3,0> > const &,class Eigen::MatrixBase<class Eigen::Matrix<float,1,2,1,1,2> > const &);

include/igl/winding_number.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "winding_number.h"
 #include "WindingNumberAABB.h"
@@ -40,10 +40,10 @@ IGL_INLINE void igl::winding_number(
     }
     case 3:
     {
-      WindingNumberAABB< 
-        Eigen::Matrix<typename DerivedV::Scalar,1,3>, 
+      WindingNumberAABB<
+        Eigen::Matrix<typename DerivedV::Scalar,1,3>,
         DerivedV,
-        DerivedF> 
+        DerivedF>
         hier(V,F);
       hier.grow();
       // loop over origins

tutorial/505_MIQ/main.cpp

@@ -15,6 +15,7 @@
 #include <sstream>
 
 #include "tutorial_shared_path.h"
+#include <igl/serialize.h>
 
 // Input mesh
 Eigen::MatrixXd V;
@@ -38,13 +39,13 @@ Eigen::MatrixXd BIS1, BIS2;
 Eigen::MatrixXd BIS1_combed, BIS2_combed;
 
 // Per-corner, integer mismatches
-Eigen::MatrixXi MMatch;
+Eigen::Matrix<int, Eigen::Dynamic, 3> MMatch;
 
 // Field singularities
-Eigen::VectorXi isSingularity, singularityIndex;
+Eigen::Matrix<int, Eigen::Dynamic, 1> isSingularity, singularityIndex;
 
 // Per corner seams
-Eigen::MatrixXi Seams;
+Eigen::Matrix<int, Eigen::Dynamic, 3> Seams;
 
 // Combed field
 Eigen::MatrixXd X1_combed, X2_combed;
@@ -234,50 +235,50 @@ int main(int argc, char *argv[])
   // Load a mesh in OFF format
   igl::readOFF(TUTORIAL_SHARED_PATH "/3holes.off", V, F);
 
+  double gradient_size = 50;
+  double iter = 0;
+  double stiffness = 5.0;
+  bool direct_round = 0;
+  std::string filename = "/Users/daniele/x.dat";
+
   // Compute face barycenters
   igl::barycenter(V, F, B);
 
   // Compute scale for visualizing fields
   global_scale =  .5*igl::avg_edge_length(V, F);
 
+    // Contrain one face
+    VectorXi b(1);
+    b << 0;
+    MatrixXd bc(1, 3);
+    bc << 1, 0, 0;
 
-  // Contrain one face
-  VectorXi b(1);
-  b << 0;
-  MatrixXd bc(1,3);
-  bc << 1, 0, 0;
-
-  // Create a smooth 4-RoSy field
-  VectorXd S;
-  igl::copyleft::comiso::nrosy(V,F,b,bc,VectorXi(),VectorXd(),MatrixXd(),4,0.5,X1,S);
+    // Create a smooth 4-RoSy field
+    VectorXd S;
+    igl::copyleft::comiso::nrosy(V, F, b, bc, VectorXi(), VectorXd(), MatrixXd(), 4, 0.5, X1, S);
 
-  // Find the the orthogonal vector
-  MatrixXd B1,B2,B3;
-  igl::local_basis(V,F,B1,B2,B3);
-  X2 = igl::rotate_vectors(X1, VectorXd::Constant(1,M_PI/2), B1, B2);
-
-  double gradient_size = 50;
-  double iter = 0;
-  double stiffness = 5.0;
-  bool direct_round = 0;
+    // Find the the orthogonal vector
+    MatrixXd B1, B2, B3;
+    igl::local_basis(V, F, B1, B2, B3);
+    X2 = igl::rotate_vectors(X1, VectorXd::Constant(1, M_PI / 2), B1, B2);
 
-  // Always work on the bisectors, it is more general
-  igl::compute_frame_field_bisectors(V, F, X1, X2, BIS1, BIS2);
+    // Always work on the bisectors, it is more general
+    igl::compute_frame_field_bisectors(V, F, X1, X2, BIS1, BIS2);
 
-  // Comb the field, implicitly defining the seams
-  igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
+    // Comb the field, implicitly defining the seams
+    igl::comb_cross_field(V, F, BIS1, BIS2, BIS1_combed, BIS2_combed);
 
-  // Find the integer mismatches
-  igl::cross_field_missmatch(V, F, BIS1_combed, BIS2_combed, true, MMatch);
+    // Find the integer mismatches
+    igl::cross_field_missmatch(V, F, BIS1_combed, BIS2_combed, true, MMatch);
 
-  // Find the singularities
-  igl::find_cross_field_singularities(V, F, MMatch, isSingularity, singularityIndex);
+    // Find the singularities
+    igl::find_cross_field_singularities(V, F, MMatch, isSingularity, singularityIndex);
 
-  // Cut the mesh, duplicating all vertices on the seams
-  igl::cut_mesh_from_singularities(V, F, MMatch, Seams);
+    // Cut the mesh, duplicating all vertices on the seams
+    igl::cut_mesh_from_singularities(V, F, MMatch, Seams);
 
-  // Comb the frame-field accordingly
-  igl::comb_frame_field(V, F, X1, X2, BIS1_combed, BIS2_combed, X1_combed, X2_combed);
+    // Comb the frame-field accordingly
+    igl::comb_frame_field(V, F, X1, X2, BIS1_combed, BIS2_combed, X1_combed, X2_combed);
 
   // Global parametrization
   igl::copyleft::comiso::miq(V,

tutorial/712_DataSmoothing/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(712_DataSmoothing)
+
+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_VIEWER_EXTRA_LIBRARIES} ${LIBIGL_OPENGL_EXTRA_LIBRARIES} ${LIBIGL_OPENGL_GLFW_EXTRA_LIBRARIES})

tutorial/712_DataSmoothing/main.cpp

@@ -0,0 +1,109 @@
+#include <igl/read_triangle_mesh.h>
+#include <igl/hessian_energy.h>
+#include <igl/massmatrix.h>
+#include <igl/cotmatrix.h>
+#include <igl/jet.h>
+#include <igl/edges.h>
+#include <igl/components.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/viewer/Viewer.h>
+
+#include <Eigen/Core>
+#include <Eigen/SparseCholesky>
+
+#include <iostream>
+#include <set>
+#include <limits>
+#include <stdlib.h>
+
+#include "tutorial_shared_path.h"
+
+#include <igl/isolines.h>
+
+
+int main(int argc, char * argv[])
+{
+    typedef Eigen::SparseMatrix<double> SparseMat;
+
+    //Read our mesh
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F, E;
+    if(!igl::read_triangle_mesh(
+        argc>1?argv[1]: TUTORIAL_SHARED_PATH "/beetle.off",V,F)) {
+        std::cout << "Failed to load mesh." << std::endl;
+    }
+    igl::edges(F,E);
+
+    //Constructing an exact function to smooth
+    Eigen::VectorXd zexact = V.block(0,2,V.rows(),1).array()
+        + 0.5*V.block(0,1,V.rows(),1).array()
+        + V.block(0,1,V.rows(),1).array().pow(2)
+        + V.block(0,2,V.rows(),1).array().pow(3);
+
+    //Make the exact function noisy
+    srand(5);
+    const double s = 0.2*(zexact.maxCoeff() - zexact.minCoeff());
+    Eigen::VectorXd znoisy = zexact + s*Eigen::VectorXd::Random(zexact.size());
+
+    //Constructing the squared Laplacian and squared Hessian energy
+    SparseMat L, M;
+    igl::cotmatrix(V, F, L);
+    igl::massmatrix(V, F, igl::MASSMATRIX_TYPE_BARYCENTRIC, M);
+    Eigen::SimplicialLDLT<SparseMat> solver(M);
+    SparseMat MinvL = solver.solve(L);
+    SparseMat QL = L.transpose()*MinvL;
+    SparseMat QH;
+    igl::hessian_energy(V, F, QH);
+
+    //Solve to find Laplacian-smoothed and Hessian-smoothed solutions
+    const double al = 8e-4;
+    Eigen::SimplicialLDLT<SparseMat> lapSolver(al*QL + (1.-al)*M);
+    Eigen::VectorXd zl = lapSolver.solve(al*M*znoisy);
+    const double ah = 5e-6;
+    Eigen::SimplicialLDLT<SparseMat> hessSolver(ah*QH + (1.-ah)*M);
+    Eigen::VectorXd zh = hessSolver.solve(ah*M*znoisy);
+
+    //Viewer that shows all functions: zexact, znoisy, zl, zh
+    igl::viewer::Viewer viewer;
+    viewer.data.set_mesh(V,F);
+    viewer.core.show_lines = false;
+    viewer.callback_key_down =
+      [&](igl::viewer::Viewer & viewer, unsigned char key, int mod)->bool {
+        //Graduate result to show isolines, then compute color matrix
+        const Eigen::VectorXd* z;
+        switch(key) {
+            case '1':
+                z = &zexact;
+                break;
+            case '2':
+                z = &znoisy;
+                break;
+            case '3':
+                z = &zl;
+                break;
+            case '4':
+                z = &zh;
+                break;
+            default:
+                return false;
+        }
+        Eigen::MatrixXd isoV;
+        Eigen::MatrixXi isoE;
+        if(key!='2')
+            igl::isolines(V, F, *z, 30, isoV, isoE);
+        viewer.data.set_edges(isoV,isoE,Eigen::RowVector3d(0,0,0));
+        Eigen::MatrixXd colors;
+        igl::jet(*z, true, colors);
+        viewer.data.set_colors(colors);
+    };
+    std::cout << R"(Usage:
+1  Show original function
+2  Show noisy function
+3  Biharmonic smoothing (zero Neumann boundary)
+4  Biharmonic smoothing (natural Hessian boundary)
+
+)";
+    viewer.launch();
+
+    return 0;
+}

tutorial/CMakeLists.txt

@@ -186,4 +186,5 @@ if(TUTORIALS_CHAPTER7)
   endif()
   add_subdirectory("710_SLIM")
   add_subdirectory("711_Subdivision")
+  add_subdirectory("712_DataSmoothing")
 endif()

tutorial/images/712_beetles.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+5513240871aa3e7a212ac8c2ae8cb466ffac3793

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-5ae3f8bf55895f273b2e1417314d6ed544096ea6
+e0514ec101376a87462a079591da84d17750f2fa

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-4868a892c78947f78ff6985ff7ea6cc069dd1b4d
+58d89e2b385d7c8849f9e3d7d1ec66e970cad92d