Merge remote-tracking branch 'upstream/master'

Former-commit-id: 231b46547274ce78b36c52ecc76401725d3a110e

.travis.yml

@@ -45,9 +45,9 @@ matrix:
     - os: osx
       compiler: clang
       script:
-        - brew update
-        - brew upgrade cmake
-        - brew upgrade cgal
+        # - brew update
+        # - brew upgrade cmake
+        # - brew upgrade cgal
         - git submodule update --init --recursive
         - cd python
         - mkdir build

README.md

@@ -43,6 +43,10 @@ and Windows with Visual Studio 2015 Community Edition.
 As of version 1.0, libigl includes an introductory
 [tutorial](http://libigl.github.io/libigl/tutorial/tutorial.html) that covers many functionalities.
 
+## libigl example project
+
+We provide a [blank project example](https://github.com/libigl/libigl-example-project) showing how to use libigl and cmake. Feel free and encouraged to copy or fork this project as a way of starting a new personal project using libigl.
+
 ## Installation
 
 Libigl is a **header-only** library. You do **not** need to build anything to
@@ -209,6 +213,7 @@ few labs/companies/institutions using libigl:
  - ETH Zurich, [Interactive Geometry Lab](http://igl.ethz.ch/) and [Advanced Technologies Lab](http://ait.inf.ethz.ch/), Swizterland
  - George Mason University, [CraGL](http://cs.gmu.edu/~ygingold/), USA
  - [Hong Kong University of Science and Technology](http://www.ust.hk/), Hong Kong
+ - [Inria](Université Grenoble Alpes), France
  - [National Institute of Informatics](http://www.nii.ac.jp/en/), Japan
  - New York University, [Media Research Lab](http://mrl.nyu.edu/), USA
  - NYUPoly, [Game Innovation Lab](http://game.engineering.nyu.edu/), USA
@@ -216,6 +221,7 @@ few labs/companies/institutions using libigl:
  - [TU Delft](http://www.tudelft.nl/en/), Netherlands
  - [TU Wien](https://www.tuwien.ac.at/en/tuwien_home/), Austria
  - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
+ - [Utrecht University](http://www.staff.science.uu.nl/~vaxma001/), The Netherlands
  - [Universidade Federal de Santa Catarina](http://mtm.ufsc.br/~leo/), Brazil
  - [University College London](http://vecg.cs.ucl.ac.uk/), England
  - [University of California Berkeley](http://vis.berkeley.edu/), USA
@@ -224,6 +230,7 @@ few labs/companies/institutions using libigl:
  - [University of Texas at Austin](http://www.cs.utexas.edu/users/evouga/), USA
  - [University of Toronto](http://dgp.toronto.edu), Canada
  - [University of Victoria](https://www.csc.uvic.ca/Research/graphics/), Canada
+ - [University of Wisconsin-Eau Claire](http://www.uwec.edu/computer-science/), USA
  - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
  - [Université Toulouse III Paul Sabatier](http://www.univ-tlse3.fr/), France
  - [Zhejiang University](http://www.math.zju.edu.cn/cagd/), China
@@ -248,7 +255,7 @@ page](https://github.com/libigl/libigl/issues).
 
 ## Copyright
 2016 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
-Zhou, Sebastian Koch, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
+Zhou, Sebastian Koch, Amir Vaxman, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.
 
 Please see individual files for appropriate copyright notices.

include/igl/AABB.cpp

@@ -803,7 +803,14 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // Actually process elements
     assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
     // Cheesecake way of hitting element
-    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hits);
+    bool ret = ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hits);
+    // Since we only gave ray_mesh_intersect a single face, it will have set
+    // any hits to id=0. Set these to this primitive's id
+    for(auto & hit : hits)
+    {
+      hit.id = m_primitive;
+    }
+    return ret;
   }
   std::vector<igl::Hit> left_hits;
   std::vector<igl::Hit> right_hits;
@@ -853,6 +860,8 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
         ray_mesh_intersect(origin,dir,V,Ele.row(tree->m_primitive),leaf_hit)&&
         leaf_hit.t < hit.t)
       {
+        // correct the id
+        leaf_hit.id = tree->m_primitive;
         hit = leaf_hit;
       }
       continue;
@@ -904,7 +913,9 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
     // Actually process elements
     assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
     // Cheesecake way of hitting element
-    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hit);
+    bool ret = ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive),hit);
+    hit.id = m_primitive;
+    return ret;
   }
 
   // Doesn't seem like smartly choosing left before/after right makes a

include/igl/active_set.cpp

@@ -77,7 +77,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   {
     lx = p_lx;
   }
-  if(ux.size() == 0)
+  if(p_ux.size() == 0)
   {
     ux = Eigen::PlainObjectBase<Derivedux>::Constant(
       n,1,numeric_limits<typename Derivedux::Scalar>::max());

include/igl/components.cpp

@@ -1,9 +1,9 @@
 // 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 
+//
+// 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 "components.h"
 #include "adjacency_matrix.h"
@@ -90,4 +90,5 @@ IGL_INLINE void igl::components(
 template void igl::components<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::components<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::components<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::components<double, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/copyleft/cgal/delaunay_triangulation.cpp

@@ -0,0 +1,24 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@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 "delaunay_triangulation.h"
+#include "../../delaunay_triangulation.h"
+#include "orient2D.h"
+#include "incircle.h"
+
+template<
+  typename DerivedV,
+  typename DerivedF>
+IGL_INLINE void igl::copyleft::cgal::delaunay_triangulation(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  igl::delaunay_triangulation(V, orient2D<Scalar>, incircle<Scalar>, F);
+}
+

include/igl/copyleft/cgal/delaunay_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_DELAUNAY_TRIANGULATION_H
+#define IGL_COPYLEFT_CGAL_DELAUNAY_TRIANGULATION_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+
+      // Given a set of points in 2D, return a Delaunay triangulation of these
+      // points.
+      //
+      // Inputs:
+      //   V  #V by 2 list of vertex positions
+      //
+      // Outputs:
+      //   F  #F by 3 of faces in Delaunay triangulation.
+      template<
+        typename DerivedV,
+        typename DerivedF
+        >
+      IGL_INLINE void delaunay_triangulation(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "delaunay_triangulation.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/incircle.cpp

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "incircle.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::incircle(
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2],
+    const Scalar pd[2])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::side_of_oriented_circle(
+        typename Kernel::Point_2(pa[0], pa[1]),
+        typename Kernel::Point_2(pb[0], pb[1]),
+        typename Kernel::Point_2(pc[0], pc[1]),
+        typename Kernel::Point_2(pd[0], pd[1]))) {
+    case CGAL::ON_POSITIVE_SIDE:
+      return 1;
+    case CGAL::ON_NEGATIVE_SIDE:
+      return -1;
+    case CGAL::ON_ORIENTED_BOUNDARY:
+      return 0;
+    default:
+      throw "Invalid incircle result";
+  }
+}

include/igl/copyleft/cgal/incircle.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_INCIRCLE_H
+#define IGL_COPYLEFT_CGAL_INCIRCLE_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd  2D points.
+      // Output:
+      //   1 if pd is inside of the oriented circle formed by pa,pb,pc.
+      //   0 if pd is co-circular with pa,pb,pc.
+      //  -1 if pd is outside of the oriented circle formed by pa,pb,pc.
+      template <typename Scalar>
+      IGL_INLINE short incircle(
+          const Scalar pa[2],
+          const Scalar pb[2],
+          const Scalar pc[2],
+          const Scalar pd[2]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "incircle.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/insphere.cpp

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "insphere.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::insphere(
+    const Scalar pa[3],
+    const Scalar pb[3],
+    const Scalar pc[3],
+    const Scalar pd[3],
+    const Scalar pe[3])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::side_of_oriented_sphere(
+        typename Kernel::Point_3(pa[0], pa[1], pa[2]),
+        typename Kernel::Point_3(pb[0], pb[1], pb[2]),
+        typename Kernel::Point_3(pc[0], pc[1], pc[2]),
+        typename Kernel::Point_3(pd[0], pd[1], pd[2]),
+        typename Kernel::Point_3(pe[0], pe[1], pe[2]))) {
+    case CGAL::ON_POSITIVE_SIDE:
+      return 1;
+    case CGAL::ON_NEGATIVE_SIDE:
+      return -1;
+    case CGAL::ON_ORIENTED_BOUNDARY:
+      return 0;
+    default:
+      throw "Invalid incircle result";
+  }
+}

include/igl/copyleft/cgal/insphere.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_INSPHERE_H
+#define IGL_COPYLEFT_CGAL_INSPHERE_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd,pe  3D points.
+      // Output:
+      //   1 if pe is inside of the oriented sphere formed by pa,pb,pc,pd.
+      //   0 if pe is co-spherical with pa,pb,pc,pd.
+      //  -1 if pe is outside of the oriented sphere formed by pa,pb,pc,pd.
+      template <typename Scalar>
+      IGL_INLINE short insphere(
+          const Scalar pa[3],
+          const Scalar pb[3],
+          const Scalar pc[3],
+          const Scalar pd[3],
+          const Scalar pe[3]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "insphere.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/lexicographic_triangulation.cpp

@@ -0,0 +1,24 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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 "lexicographic_triangulation.h"
+#include "../../lexicographic_triangulation.h"
+#include "orient2D.h"
+
+template<
+  typename DerivedP,
+  typename DerivedF
+  >
+IGL_INLINE void igl::copyleft::cgal::lexicographic_triangulation(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  igl::lexicographic_triangulation(P, orient2D<Scalar>, F);
+}

include/igl/copyleft/cgal/lexicographic_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_LEXICOGRAPHIC_TRIANGULATION_H
+#define IGL_COPYLEFT_CGAL_LEXICOGRAPHIC_TRIANGULATION_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+
+      // Given a set of points in 2D, return a lexicographic triangulation of these
+      // points.
+      //
+      // Inputs:
+      //   P  #P by 2 list of vertex positions
+      //
+      // Outputs:
+      //   F  #F by 3 of faces in lexicographic triangulation.
+      template<
+        typename DerivedP,
+        typename DerivedF
+        >
+      IGL_INLINE void lexicographic_triangulation(
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lexicographic_triangulation.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/orient2D.cpp

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "orient2D.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::orient2D(
+    const Scalar pa[2],
+    const Scalar pb[2],
+    const Scalar pc[2])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::orientation(
+        typename Kernel::Point_2(pa[0], pa[1]),
+        typename Kernel::Point_2(pb[0], pb[1]),
+        typename Kernel::Point_2(pc[0], pc[1]))) {
+    case CGAL::LEFT_TURN:
+      return 1;
+    case CGAL::RIGHT_TURN:
+      return -1;
+    case CGAL::COLLINEAR:
+      return 0;
+    default:
+      throw "Invalid orientation";
+  }
+}

include/igl/copyleft/cgal/orient2D.h

@@ -0,0 +1,38 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_ORIENT_2D_H
+#define IGL_COPYLEFT_CGAL_ORIENT_2D_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc   2D points.
+      // Output:
+      //   1 if pa,pb,pc are counterclockwise oriented.
+      //   0 if pa,pb,pc are counterclockwise oriented.
+      //  -1 if pa,pb,pc are clockwise oriented.
+      template <typename Scalar>
+      IGL_INLINE short orient2D(
+          const Scalar pa[2],
+          const Scalar pb[2],
+          const Scalar pc[2]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient2D.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/orient3D.cpp

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "orient3D.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+template<typename Scalar>
+IGL_INLINE short igl::copyleft::cgal::orient3D(
+    const Scalar pa[3],
+    const Scalar pb[3],
+    const Scalar pc[3],
+    const Scalar pd[3])
+{
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Epeck;
+  typedef CGAL::Exact_predicates_inexact_constructions_kernel Epick;
+  typedef typename std::conditional<std::is_same<Scalar, Epeck::FT>::value,
+          Epeck, Epick>::type Kernel;
+
+  switch(CGAL::orientation(
+        typename Kernel::Point_3(pa[0], pa[1], pa[2]),
+        typename Kernel::Point_3(pb[0], pb[1], pb[2]),
+        typename Kernel::Point_3(pc[0], pc[1], pc[2]),
+        typename Kernel::Point_3(pd[0], pd[1], pd[2]))) {
+    case CGAL::POSITIVE:
+      return 1;
+    case CGAL::NEGATIVE:
+      return -1;
+    case CGAL::COPLANAR:
+      return 0;
+    default:
+      throw "Invalid orientation";
+  }
+}

include/igl/copyleft/cgal/orient3D.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_COPYLEFT_CGAL_ORIENT_3D_H
+#define IGL_COPYLEFT_CGAL_ORIENT_3D_H
+
+#include "../../igl_inline.h"
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Inputs:
+      //   pa,pb,pc,pd  3D points.
+      // Output:
+      //   1 if pa,pb,pc,pd forms a tet of positive volume.
+      //   0 if pa,pb,pc,pd are coplanar.
+      //  -1 if pa,pb,pc,pd forms a tet of negative volume.
+      template <typename Scalar>
+      IGL_INLINE short orient3D(
+          const Scalar pa[3],
+          const Scalar pb[3],
+          const Scalar pc[3],
+          const Scalar pd[3]);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "orient3D.cpp"
+#endif
+#endif

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -43,7 +43,7 @@ void igl::crouzeix_raviart_massmatrix(
   using namespace std;
   assert(F.cols() == 3);
   // Mesh should be edge-manifold
-  assert(is_edge_manifold(V,F));
+  assert(is_edge_manifold(F));
   // number of elements (triangles)
   int m = F.rows();
   // Get triangle areas

include/igl/delaunay_triangulation.cpp

@@ -0,0 +1,83 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@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 "delaunay_triangulation.h"
+#include "flip_edge.h"
+#include "lexicographic_triangulation.h"
+#include "unique_edge_map.h"
+
+#include <vector>
+#include <sstream>
+
+template<
+  typename DerivedV,
+  typename Orient2D,
+  typename InCircle,
+  typename DerivedF>
+IGL_INLINE void igl::delaunay_triangulation(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Orient2D orient2D,
+    InCircle incircle,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  assert(V.cols() == 2);
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedV::Scalar Scalar;
+  igl::lexicographic_triangulation(V, orient2D, F);
+  const size_t num_faces = F.rows();
+  if (num_faces == 0) {
+    // Input points are degenerate.  No faces will be generated.
+    return;
+  }
+  assert(F.cols() == 3);
+
+  Eigen::MatrixXi E;
+  Eigen::MatrixXi uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<Index> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+  auto is_delaunay = [&V,&F,&uE2E,num_faces,&incircle](size_t uei) {
+    auto& half_edges = uE2E[uei];
+    if (half_edges.size() != 2) {
+      throw "Cannot flip non-manifold or boundary edge";
+    }
+
+    const size_t f1 = half_edges[0] % num_faces;
+    const size_t f2 = half_edges[1] % num_faces;
+    const size_t c1 = half_edges[0] / num_faces;
+    const size_t c2 = half_edges[1] / num_faces;
+    assert(c1 < 3);
+    assert(c2 < 3);
+    assert(f1 != f2);
+    const size_t v1 = F(f1, (c1+1)%3);
+    const size_t v2 = F(f1, (c1+2)%3);
+    const size_t v4 = F(f1, c1);
+    const size_t v3 = F(f2, c2);
+    const Scalar p1[] = {V(v1, 0), V(v1, 1)};
+    const Scalar p2[] = {V(v2, 0), V(v2, 1)};
+    const Scalar p3[] = {V(v3, 0), V(v3, 1)};
+    const Scalar p4[] = {V(v4, 0), V(v4, 1)};
+    auto orientation = incircle(p1, p2, p4, p3);
+    return orientation <= 0;
+  };
+
+  bool all_delaunay = false;
+  while(!all_delaunay) {
+    all_delaunay = true;
+    for (size_t i=0; i<uE2E.size(); i++) {
+      if (uE2E[i].size() == 2) {
+        if (!is_delaunay(i)) {
+          all_delaunay = false;
+          flip_edge(F, E, uE, EMAP, uE2E, i);
+        }
+      }
+    }
+  }
+}
+

include/igl/delaunay_triangulation.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_DELAUNAY_TRIANGULATION_H
+#define IGL_DELAUNAY_TRIANGULATION_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Given a set of points in 2D, return a Delaunay triangulation of these
+  // points.
+  //
+  // Inputs:
+  //   V  #V by 2 list of vertex positions
+  //   orient2D  A functor such that orient2D(pa, pb, pc) returns
+  //               1 if pa,pb,pc forms a conterclockwise triangle.
+  //              -1 if pa,pb,pc forms a clockwise triangle.
+  //               0 if pa,pb,pc are collinear.
+  //              where the argument pa,pb,pc are of type Scalar[2].
+  //   incircle  A functor such that incircle(pa, pb, pc, pd) returns
+  //               1 if pd is on the positive size of circumcirle of (pa,pb,pc)
+  //              -1 if pd is on the positive size of circumcirle of (pa,pb,pc)
+  //               0 if pd is cocircular with pa, pb, pc.
+  // Outputs:
+  //   F  #F by 3 of faces in Delaunay triangulation.
+  template<
+    typename DerivedV,
+    typename Orient2D,
+    typename InCircle,
+    typename DerivedF
+    >
+  IGL_INLINE void delaunay_triangulation(
+      const Eigen::PlainObjectBase<DerivedV>& V,
+      Orient2D orient2D,
+      InCircle incircle,
+      Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "delaunay_triangulation.cpp"
+#endif
+#endif

include/igl/doublearea.cpp

@@ -27,8 +27,6 @@ IGL_INLINE void igl::doublearea(
   const size_t m = F.rows();
   // Compute edge lengths
   Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> l;
-  // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
-  // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
 
   // Projected area helper
   const auto & proj_doublearea =
@@ -143,6 +141,11 @@ IGL_INLINE void igl::doublearea(
   const Index m = ul.rows();
   Eigen::Matrix<typename Derivedl::Scalar, Eigen::Dynamic, 3> l;
   MatrixXi _;
+  // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
+  // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
+  //
+  // "Miscalculating Area and Angles of a Needle-like Triangle"
+  // https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf
   sort(ul,2,false,l,_);
   // semiperimeters
   Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
@@ -198,6 +201,8 @@ Eigen::PlainObjectBase<DeriveddblA> & dblA)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::doublearea<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<double, -1, 1, 0, -1, 1> >(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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

include/igl/edge_flaps.cpp

@@ -17,8 +17,9 @@ IGL_INLINE void igl::edge_flaps(
   Eigen::MatrixXi & EF,
   Eigen::MatrixXi & EI)
 {
-  EF.resize(E.rows(),2);
-  EI.resize(E.rows(),2);
+  // Initialize to boundary value
+  EF.setConstant(E.rows(),2,-1);
+  EI.setConstant(E.rows(),2,-1);
   // loop over all faces
   for(int f = 0;f<F.rows();f++)
   {

include/igl/edge_topology.cpp

@@ -6,16 +6,16 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "edge_topology.h"
-#include <algorithm>
 #include "is_edge_manifold.h"
+#include <algorithm>
 
 template<typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::edge_topology(
-                                   const Eigen::PlainObjectBase<DerivedV>& V,
-                                   const Eigen::PlainObjectBase<DerivedF>& F,
-                                   Eigen::MatrixXi& EV,
-                                   Eigen::MatrixXi& FE,
-                                   Eigen::MatrixXi& EF)
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::MatrixXi& EV,
+  Eigen::MatrixXi& FE,
+  Eigen::MatrixXi& EF)
 {
   // Only needs to be edge-manifold
   if (V.rows() ==0 || F.rows()==0)
@@ -25,7 +25,7 @@ IGL_INLINE void igl::edge_topology(
     EF = Eigen::MatrixXi::Constant(0,2,-1);
     return;
   }
-  assert(igl::is_edge_manifold(V,F));
+  assert(igl::is_edge_manifold(F));
   std::vector<std::vector<int> > ETT;
   for(int f=0;f<F.rows();++f)
     for (int i=0;i<3;++i)

include/igl/edge_topology.h

@@ -15,14 +15,21 @@
 
 namespace igl 
 {
-  // Initialize Edges and their topological relations
+  // Initialize Edges and their topological relations (assumes an edge-manifold
+  // mesh)
   //
   // Output:
   // EV  : #Ex2, Stores the edge description as pair of indices to vertices
   // FE : #Fx3, Stores the Triangle-Edge relation
   // EF : #Ex2: Stores the Edge-Triangle relation
   //
-  // TODO: This seems to be a duplicate of edge_flaps.h
+  // TODO: This seems to be a inferior duplicate of edge_flaps.h:
+  //   - unused input parameter V
+  //   - roughly 2x slower than edge_flaps
+  //   - outputs less information: edge_flaps reveals corner opposite edge
+  //   - FE uses non-standard and ambiguous order: FE(f,c) is merely an edge
+  //     incident on corner c of face f. In contrast, edge_flaps's EMAP(f,c) reveals
+  //     the edge _opposite_ corner c of face f
 template <typename DerivedV, typename DerivedF>
   IGL_INLINE void edge_topology(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/euler_characteristic.cpp

@@ -29,4 +29,5 @@ IGL_INLINE int igl::euler_characteristic(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template int igl::euler_characteristic<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&);
 #endif

include/igl/euler_characteristic.h

@@ -1,6 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
-// Copyright (C) 2016 Michael Rabinovich <michaelrabinovich27@gmail.com@gmail.com>
+// Copyright (C) 2016 Michael Rabinovich <michaelrabinovich27@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

include/igl/extract_manifold_patches.cpp

@@ -13,7 +13,7 @@ IGL_INLINE size_t igl::extract_manifold_patches(
   const Eigen::PlainObjectBase<DerivedF>& F,
   const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
   const std::vector<std::vector<uE2EType> >& uE2E,
-  Eigen::PlainObjectBase<DerivedP>& P) 
+  Eigen::PlainObjectBase<DerivedP>& P)
 {
     assert(F.cols() == 3);
     const size_t num_faces = F.rows();
@@ -74,7 +74,7 @@ template<
   typename DerivedP>
 IGL_INLINE size_t igl::extract_manifold_patches(
   const Eigen::PlainObjectBase<DerivedF>& F,
-  Eigen::PlainObjectBase<DerivedP>& P) 
+  Eigen::PlainObjectBase<DerivedP>& P)
 {
   Eigen::MatrixXi E, uE;
   Eigen::VectorXi EMAP;
@@ -84,5 +84,7 @@ IGL_INLINE size_t igl::extract_manifold_patches(
 }
 
 #ifdef IGL_STATIC_LIBRARY
+#ifndef WIN32
 template unsigned long igl::extract_manifold_patches<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif
+#endif

include/igl/flip_avoiding_line_search.cpp

@@ -0,0 +1,272 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 "flip_avoiding_line_search.h"
+
+#include <Eigen/Dense>
+#include <vector>
+#include "line_search.h"
+#define TwoPi  2*M_PI//6.28318530717958648
+
+using namespace std;
+
+//---------------------------------------------------------------------------
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+//         2 real roots: x[0], x[1],          return 2
+//         1 real root : x[0], x[1] ± i*x[2], return 1
+// http://math.ivanovo.ac.ru/dalgebra/Khashin/poly/index.html
+int SolveP3(std::vector<double>& x,double a,double b,double c) { // solve cubic equation x^3 + a*x^2 + b*x + c
+  double a2 = a*a;
+    double q  = (a2 - 3*b)/9;
+  double r  = (a*(2*a2-9*b) + 27*c)/54;
+    double r2 = r*r;
+  double q3 = q*q*q;
+  double A,B;
+    if(r2<q3) {
+        double t=r/sqrt(q3);
+    if( t<-1) t=-1;
+    if( t> 1) t= 1;
+        t=acos(t);
+        a/=3; q=-2*sqrt(q);
+        x[0]=q*cos(t/3)-a;
+        x[1]=q*cos((t+TwoPi)/3)-a;
+        x[2]=q*cos((t-TwoPi)/3)-a;
+        return(3);
+    } else {
+        A =-pow(fabs(r)+sqrt(r2-q3),1./3);
+    if( r<0 ) A=-A;
+    B = A==0? 0 : B=q/A;
+
+    a/=3;
+    x[0] =(A+B)-a;
+        x[1] =-0.5*(A+B)-a;
+        x[2] = 0.5*sqrt(3.)*(A-B);
+    if(fabs(x[2])<1e-14) { x[2]=x[1]; return(2); }
+        return(1);
+    }
+}
+
+double get_smallest_pos_quad_zero(double a,double b, double c) {
+  double t1,t2;
+  if (a != 0) {
+    double delta_in = pow(b,2) - 4*a*c;
+    if (delta_in < 0) {
+      return INFINITY;
+    }
+    double delta = sqrt(delta_in);
+    t1 = (-b + delta)/ (2*a);
+    t2 = (-b - delta)/ (2*a);
+  } else {
+    t1 = t2 = -b/c;
+  }
+  assert (std::isfinite(t1));
+  assert (std::isfinite(t2));
+
+  double tmp_n = min(t1,t2);
+  t1 = max(t1,t2); t2 = tmp_n;
+  if (t1 == t2) {
+    return INFINITY; // means the orientation flips twice = doesn't flip
+  }
+  // return the smallest negative root if it exists, otherwise return infinity
+  if (t1 > 0) {
+    if (t2 > 0) {
+      return t2;
+    } else {
+      return t1;
+    }
+  } else {
+    return INFINITY;
+  }
+}
+
+ double get_min_pos_root_2D(const Eigen::MatrixXd& uv,const Eigen::MatrixXi& F,
+            Eigen::MatrixXd& d, int f) {
+/*
+      Finding the smallest timestep t s.t a triangle get degenerated (<=> det = 0)
+      The following code can be derived by a symbolic expression in matlab:
+
+      Symbolic matlab:
+      U11 = sym('U11');
+      U12 = sym('U12');
+      U21 = sym('U21');
+      U22 = sym('U22');
+      U31 = sym('U31');
+      U32 = sym('U32');
+
+      V11 = sym('V11');
+      V12 = sym('V12');
+      V21 = sym('V21');
+      V22 = sym('V22');
+      V31 = sym('V31');
+      V32 = sym('V32');
+
+      t = sym('t');
+
+      U1 = [U11,U12];
+      U2 = [U21,U22];
+      U3 = [U31,U32];
+
+      V1 = [V11,V12];
+      V2 = [V21,V22];
+      V3 = [V31,V32];
+
+      A = [(U2+V2*t) - (U1+ V1*t)];
+      B = [(U3+V3*t) - (U1+ V1*t)];
+      C = [A;B];
+
+      solve(det(C), t);
+      cf = coeffs(det(C),t); % Now cf(1),cf(2),cf(3) holds the coefficients for the polynom. at order c,b,a
+    */
+
+  int v1 = F(f,0); int v2 = F(f,1); int v3 = F(f,2);
+  // get quadratic coefficients (ax^2 + b^x + c)
+  #define U11 uv(v1,0)
+  #define U12 uv(v1,1)
+  #define U21 uv(v2,0)
+  #define U22 uv(v2,1)
+  #define U31 uv(v3,0)
+  #define U32 uv(v3,1)
+
+  #define V11 d(v1,0)
+  #define V12 d(v1,1)
+  #define V21 d(v2,0)
+  #define V22 d(v2,1)
+  #define V31 d(v3,0)
+  #define V32 d(v3,1)
+
+
+  double a = V11*V22 - V12*V21 - V11*V32 + V12*V31 + V21*V32 - V22*V31;
+  double b = U11*V22 - U12*V21 - U21*V12 + U22*V11 - U11*V32 + U12*V31 + U31*V12 - U32*V11 + U21*V32 - U22*V31 - U31*V22 + U32*V21;
+  double c = U11*U22 - U12*U21 - U11*U32 + U12*U31 + U21*U32 - U22*U31;
+
+  return get_smallest_pos_quad_zero(a,b,c);
+}
+
+double get_min_pos_root_3D(const Eigen::MatrixXd& uv,const Eigen::MatrixXi& F,
+            Eigen::MatrixXd& direc, int f) {
+  /*
+      Searching for the roots of:
+        +-1/6 * |ax ay az 1|
+                |bx by bz 1|
+                |cx cy cz 1|
+                |dx dy dz 1|
+      Every point ax,ay,az has a search direction a_dx,a_dy,a_dz, and so we add those to the matrix, and solve the cubic to find the step size t for a 0 volume
+      Symbolic matlab:
+        syms a_x a_y a_z a_dx a_dy a_dz % tetrahedera point and search direction
+        syms b_x b_y b_z b_dx b_dy b_dz
+        syms c_x c_y c_z c_dx c_dy c_dz
+        syms d_x d_y d_z d_dx d_dy d_dz
+        syms t % Timestep var, this is what we're looking for
+
+
+        a_plus_t = [a_x,a_y,a_z] + t*[a_dx,a_dy,a_dz];
+        b_plus_t = [b_x,b_y,b_z] + t*[b_dx,b_dy,b_dz];
+        c_plus_t = [c_x,c_y,c_z] + t*[c_dx,c_dy,c_dz];
+        d_plus_t = [d_x,d_y,d_z] + t*[d_dx,d_dy,d_dz];
+
+        vol_mat = [a_plus_t,1;b_plus_t,1;c_plus_t,1;d_plus_t,1]
+        //cf = coeffs(det(vol_det),t); % Now cf(1),cf(2),cf(3),cf(4) holds the coefficients for the polynom
+        [coefficients,terms] = coeffs(det(vol_det),t); % terms = [ t^3, t^2, t, 1], Coefficients hold the coeff we seek
+  */
+  int v1 = F(f,0); int v2 = F(f,1); int v3 = F(f,2); int v4 = F(f,3);
+  #define a_x uv(v1,0)
+  #define a_y uv(v1,1)
+  #define a_z uv(v1,2)
+  #define b_x uv(v2,0)
+  #define b_y uv(v2,1)
+  #define b_z uv(v2,2)
+  #define c_x uv(v3,0)
+  #define c_y uv(v3,1)
+  #define c_z uv(v3,2)
+  #define d_x uv(v4,0)
+  #define d_y uv(v4,1)
+  #define d_z uv(v4,2)
+
+  #define a_dx direc(v1,0)
+  #define a_dy direc(v1,1)
+  #define a_dz direc(v1,2)
+  #define b_dx direc(v2,0)
+  #define b_dy direc(v2,1)
+  #define b_dz direc(v2,2)
+  #define c_dx direc(v3,0)
+  #define c_dy direc(v3,1)
+  #define c_dz direc(v3,2)
+  #define d_dx direc(v4,0)
+  #define d_dy direc(v4,1)
+  #define d_dz direc(v4,2)
+
+  // Find solution for: a*t^3 + b*t^2 + c*d +d = 0
+  double a = a_dx*b_dy*c_dz - a_dx*b_dz*c_dy - a_dy*b_dx*c_dz + a_dy*b_dz*c_dx + a_dz*b_dx*c_dy - a_dz*b_dy*c_dx - a_dx*b_dy*d_dz + a_dx*b_dz*d_dy + a_dy*b_dx*d_dz - a_dy*b_dz*d_dx - a_dz*b_dx*d_dy + a_dz*b_dy*d_dx + a_dx*c_dy*d_dz - a_dx*c_dz*d_dy - a_dy*c_dx*d_dz + a_dy*c_dz*d_dx + a_dz*c_dx*d_dy - a_dz*c_dy*d_dx - b_dx*c_dy*d_dz + b_dx*c_dz*d_dy + b_dy*c_dx*d_dz - b_dy*c_dz*d_dx - b_dz*c_dx*d_dy + b_dz*c_dy*d_dx;
+  double b = a_dy*b_dz*c_x - a_dy*b_x*c_dz - a_dz*b_dy*c_x + a_dz*b_x*c_dy + a_x*b_dy*c_dz - a_x*b_dz*c_dy - a_dx*b_dz*c_y + a_dx*b_y*c_dz + a_dz*b_dx*c_y - a_dz*b_y*c_dx - a_y*b_dx*c_dz + a_y*b_dz*c_dx + a_dx*b_dy*c_z - a_dx*b_z*c_dy - a_dy*b_dx*c_z + a_dy*b_z*c_dx + a_z*b_dx*c_dy - a_z*b_dy*c_dx - a_dy*b_dz*d_x + a_dy*b_x*d_dz + a_dz*b_dy*d_x - a_dz*b_x*d_dy - a_x*b_dy*d_dz + a_x*b_dz*d_dy + a_dx*b_dz*d_y - a_dx*b_y*d_dz - a_dz*b_dx*d_y + a_dz*b_y*d_dx + a_y*b_dx*d_dz - a_y*b_dz*d_dx - a_dx*b_dy*d_z + a_dx*b_z*d_dy + a_dy*b_dx*d_z - a_dy*b_z*d_dx - a_z*b_dx*d_dy + a_z*b_dy*d_dx + a_dy*c_dz*d_x - a_dy*c_x*d_dz - a_dz*c_dy*d_x + a_dz*c_x*d_dy + a_x*c_dy*d_dz - a_x*c_dz*d_dy - a_dx*c_dz*d_y + a_dx*c_y*d_dz + a_dz*c_dx*d_y - a_dz*c_y*d_dx - a_y*c_dx*d_dz + a_y*c_dz*d_dx + a_dx*c_dy*d_z - a_dx*c_z*d_dy - a_dy*c_dx*d_z + a_dy*c_z*d_dx + a_z*c_dx*d_dy - a_z*c_dy*d_dx - b_dy*c_dz*d_x + b_dy*c_x*d_dz + b_dz*c_dy*d_x - b_dz*c_x*d_dy - b_x*c_dy*d_dz + b_x*c_dz*d_dy + b_dx*c_dz*d_y - b_dx*c_y*d_dz - b_dz*c_dx*d_y + b_dz*c_y*d_dx + b_y*c_dx*d_dz - b_y*c_dz*d_dx - b_dx*c_dy*d_z + b_dx*c_z*d_dy + b_dy*c_dx*d_z - b_dy*c_z*d_dx - b_z*c_dx*d_dy + b_z*c_dy*d_dx;
+  double c = a_dz*b_x*c_y - a_dz*b_y*c_x - a_x*b_dz*c_y + a_x*b_y*c_dz + a_y*b_dz*c_x - a_y*b_x*c_dz - a_dy*b_x*c_z + a_dy*b_z*c_x + a_x*b_dy*c_z - a_x*b_z*c_dy - a_z*b_dy*c_x + a_z*b_x*c_dy + a_dx*b_y*c_z - a_dx*b_z*c_y - a_y*b_dx*c_z + a_y*b_z*c_dx + a_z*b_dx*c_y - a_z*b_y*c_dx - a_dz*b_x*d_y + a_dz*b_y*d_x + a_x*b_dz*d_y - a_x*b_y*d_dz - a_y*b_dz*d_x + a_y*b_x*d_dz + a_dy*b_x*d_z - a_dy*b_z*d_x - a_x*b_dy*d_z + a_x*b_z*d_dy + a_z*b_dy*d_x - a_z*b_x*d_dy - a_dx*b_y*d_z + a_dx*b_z*d_y + a_y*b_dx*d_z - a_y*b_z*d_dx - a_z*b_dx*d_y + a_z*b_y*d_dx + a_dz*c_x*d_y - a_dz*c_y*d_x - a_x*c_dz*d_y + a_x*c_y*d_dz + a_y*c_dz*d_x - a_y*c_x*d_dz - a_dy*c_x*d_z + a_dy*c_z*d_x + a_x*c_dy*d_z - a_x*c_z*d_dy - a_z*c_dy*d_x + a_z*c_x*d_dy + a_dx*c_y*d_z - a_dx*c_z*d_y - a_y*c_dx*d_z + a_y*c_z*d_dx + a_z*c_dx*d_y - a_z*c_y*d_dx - b_dz*c_x*d_y + b_dz*c_y*d_x + b_x*c_dz*d_y - b_x*c_y*d_dz - b_y*c_dz*d_x + b_y*c_x*d_dz + b_dy*c_x*d_z - b_dy*c_z*d_x - b_x*c_dy*d_z + b_x*c_z*d_dy + b_z*c_dy*d_x - b_z*c_x*d_dy - b_dx*c_y*d_z + b_dx*c_z*d_y + b_y*c_dx*d_z - b_y*c_z*d_dx - b_z*c_dx*d_y + b_z*c_y*d_dx;
+  double d = a_x*b_y*c_z - a_x*b_z*c_y - a_y*b_x*c_z + a_y*b_z*c_x + a_z*b_x*c_y - a_z*b_y*c_x - a_x*b_y*d_z + a_x*b_z*d_y + a_y*b_x*d_z - a_y*b_z*d_x - a_z*b_x*d_y + a_z*b_y*d_x + a_x*c_y*d_z - a_x*c_z*d_y - a_y*c_x*d_z + a_y*c_z*d_x + a_z*c_x*d_y - a_z*c_y*d_x - b_x*c_y*d_z + b_x*c_z*d_y + b_y*c_x*d_z - b_y*c_z*d_x - b_z*c_x*d_y + b_z*c_y*d_x;
+
+  if (a==0) {
+    return get_smallest_pos_quad_zero(b,c,d);
+  }
+  b/=a; c/=a; d/=a; // normalize it all
+  std::vector<double> res(3);
+  int real_roots_num = SolveP3(res,b,c,d);
+  switch (real_roots_num) {
+    case 1:
+      return (res[0] >= 0) ? res[0]:INFINITY;
+    case 2: {
+      double max_root = max(res[0],res[1]); double min_root = min(res[0],res[1]);
+      if (min_root > 0) return min_root;
+      if (max_root > 0) return max_root;
+      return INFINITY;
+    }
+    case 3:
+    default: {
+      std::sort(res.begin(),res.end());
+      if (res[0] > 0) return res[0];
+      if (res[1] > 0) return res[1];
+      if (res[2] > 0) return res[2];
+      return INFINITY;
+    }
+  }
+
+}
+
+double compute_max_step_from_singularities(const Eigen::MatrixXd& uv,
+                                            const Eigen::MatrixXi& F,
+                                            Eigen::MatrixXd& d) {
+    double max_step = INFINITY;
+
+    // The if statement is outside the for loops to avoid branching/ease parallelizing
+    if (uv.cols() == 2) {
+      for (int f = 0; f < F.rows(); f++) {
+        double min_positive_root = get_min_pos_root_2D(uv,F,d,f);
+        max_step = min(max_step, min_positive_root);
+      }
+    } else { // volumetric deformation
+      for (int f = 0; f < F.rows(); f++) {
+        double min_positive_root = get_min_pos_root_3D(uv,F,d,f);
+        max_step = min(max_step, min_positive_root);
+      }
+    }
+    return max_step;
+ }
+
+IGL_INLINE double igl::flip_avoiding_line_search(
+  const Eigen::MatrixXi F,
+  Eigen::MatrixXd& cur_v,
+  Eigen::MatrixXd& dst_v,
+  std::function<double(Eigen::MatrixXd&)> energy,
+  double cur_energy)
+{
+    Eigen::MatrixXd d = dst_v - cur_v;
+
+    double min_step_to_singularity = compute_max_step_from_singularities(cur_v,F,d);
+    double max_step_size = min(1., min_step_to_singularity*0.8);
+
+    return igl::line_search(cur_v,d,max_step_size, energy, cur_energy);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/flip_avoiding_line_search.h

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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_FLIP_AVOIDING_LINE_SEARCH_H
+#define IGL_FLIP_AVOIDING_LINE_SEARCH_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // A bisection line search for a mesh based energy that avoids triangle flips as suggested in 
+  // 		"Bijective Parameterization with Free Boundaries" (Smith J. and Schaefer S., 2015).
+  //
+  // The user specifies an initial vertices position (that has no flips) and target one (that my have flipped triangles).
+  // This method first computes the largest step in direction of the destination vertices that does not incur flips,
+  // and then minimizes a given energy using this maximal step and a bisection linesearch (see igl::line_search).
+  //
+  // Supports both triangle and tet meshes.
+  //
+  // Inputs:
+  //   F  #F by 3/4 				list of mesh faces or tets
+  //   cur_v  						#V by dim list of variables
+  //   dst_v  						#V by dim list of target vertices. This mesh may have flipped triangles
+  //   energy       			    A function to compute the mesh-based energy (return an energy that is bigger than 0)
+  //   cur_energy(OPTIONAL)         The energy at the given point. Helps save redundant computations. 
+  //							    This is optional. If not specified, the function will compute it.
+  // Outputs:
+  //		cur_v  						#V by dim list of variables at the new location
+  // Returns the energy at the new point
+  IGL_INLINE double flip_avoiding_line_search(
+    const Eigen::MatrixXi F,
+    Eigen::MatrixXd& cur_v,
+    Eigen::MatrixXd& dst_v,
+    std::function<double(Eigen::MatrixXd&)> energy,
+    double cur_energy = -1);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flip_avoiding_line_search.cpp"
+#endif
+
+#endif

include/igl/flip_edge.cpp

@@ -0,0 +1,148 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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 "flip_edge.h"
+
+template <
+  typename DerivedF,
+  typename DerivedE,
+  typename DeriveduE,
+  typename DerivedEMAP,
+  typename uE2EType>
+IGL_INLINE void igl::flip_edge(
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DeriveduE> & uE,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+  std::vector<std::vector<uE2EType> > & uE2E,
+  const size_t uei)
+{
+  typedef typename DerivedF::Scalar Index;
+  const size_t num_faces = F.rows();
+  assert(F.cols() == 3);
+  //          v1                 v1
+  //          /|\                / \
+  //         / | \              /f1 \
+  //     v3 /f2|f1\ v4  =>  v3 /_____\ v4
+  //        \  |  /            \ f2  /
+  //         \ | /              \   /
+  //          \|/                \ /
+  //          v2                 v2
+  auto& half_edges = uE2E[uei];
+  if (half_edges.size() != 2) {
+    throw "Cannot flip non-manifold or boundary edge";
+  }
+
+  const size_t f1 = half_edges[0] % num_faces;
+  const size_t f2 = half_edges[1] % num_faces;
+  const size_t c1 = half_edges[0] / num_faces;
+  const size_t c2 = half_edges[1] / num_faces;
+  assert(c1 < 3);
+  assert(c2 < 3);
+
+  assert(f1 != f2);
+  const size_t v1 = F(f1, (c1+1)%3);
+  const size_t v2 = F(f1, (c1+2)%3);
+  const size_t v4 = F(f1, c1);
+  const size_t v3 = F(f2, c2);
+  assert(F(f2, (c2+2)%3) == v1);
+  assert(F(f2, (c2+1)%3) == v2);
+
+  const size_t e_12 = half_edges[0];
+  const size_t e_24 = f1 + ((c1 + 1) % 3) * num_faces;
+  const size_t e_41 = f1 + ((c1 + 2) % 3) * num_faces;
+  const size_t e_21 = half_edges[1];
+  const size_t e_13 = f2 + ((c2 + 1) % 3) * num_faces;
+  const size_t e_32 = f2 + ((c2 + 2) % 3) * num_faces;
+  assert(E(e_12, 0) == v1);
+  assert(E(e_12, 1) == v2);
+  assert(E(e_24, 0) == v2);
+  assert(E(e_24, 1) == v4);
+  assert(E(e_41, 0) == v4);
+  assert(E(e_41, 1) == v1);
+  assert(E(e_21, 0) == v2);
+  assert(E(e_21, 1) == v1);
+  assert(E(e_13, 0) == v1);
+  assert(E(e_13, 1) == v3);
+  assert(E(e_32, 0) == v3);
+  assert(E(e_32, 1) == v2);
+
+  const size_t ue_24 = EMAP[e_24];
+  const size_t ue_41 = EMAP[e_41];
+  const size_t ue_13 = EMAP[e_13];
+  const size_t ue_32 = EMAP[e_32];
+
+  F(f1, 0) = v1;
+  F(f1, 1) = v3;
+  F(f1, 2) = v4;
+  F(f2, 0) = v2;
+  F(f2, 1) = v4;
+  F(f2, 2) = v3;
+
+  uE(uei, 0) = v3;
+  uE(uei, 1) = v4;
+
+  const size_t new_e_34 = f1;
+  const size_t new_e_41 = f1 + num_faces;
+  const size_t new_e_13 = f1 + num_faces*2;
+  const size_t new_e_43 = f2;
+  const size_t new_e_32 = f2 + num_faces;
+  const size_t new_e_24 = f2 + num_faces*2;
+
+  E(new_e_34, 0) = v3;
+  E(new_e_34, 1) = v4;
+  E(new_e_41, 0) = v4;
+  E(new_e_41, 1) = v1;
+  E(new_e_13, 0) = v1;
+  E(new_e_13, 1) = v3;
+  E(new_e_43, 0) = v4;
+  E(new_e_43, 1) = v3;
+  E(new_e_32, 0) = v3;
+  E(new_e_32, 1) = v2;
+  E(new_e_24, 0) = v2;
+  E(new_e_24, 1) = v4;
+
+  EMAP[new_e_34] = uei;
+  EMAP[new_e_43] = uei;
+  EMAP[new_e_41] = ue_41;
+  EMAP[new_e_13] = ue_13;
+  EMAP[new_e_32] = ue_32;
+  EMAP[new_e_24] = ue_24;
+
+  auto replace = [](std::vector<Index>& array, Index old_v, Index new_v) {
+    std::replace(array.begin(), array.end(), old_v, new_v);
+  };
+  replace(uE2E[uei], e_12, new_e_34);
+  replace(uE2E[uei], e_21, new_e_43);
+  replace(uE2E[ue_13], e_13, new_e_13);
+  replace(uE2E[ue_32], e_32, new_e_32);
+  replace(uE2E[ue_24], e_24, new_e_24);
+  replace(uE2E[ue_41], e_41, new_e_41);
+
+#ifndef NDEBUG
+  auto sanity_check = [&](size_t ue) {
+    const auto& adj_faces = uE2E[ue];
+    if (adj_faces.size() == 2) {
+      const size_t first_f  = adj_faces[0] % num_faces;
+      const size_t first_c  = adj_faces[0] / num_faces;
+      const size_t second_f = adj_faces[1] % num_faces;
+      const size_t second_c = adj_faces[1] / num_faces;
+      const size_t vertex_0 = F(first_f, (first_c+1) % 3);
+      const size_t vertex_1 = F(first_f, (first_c+2) % 3);
+      assert(vertex_0 == F(second_f, (second_c+2) % 3));
+      assert(vertex_1 == F(second_f, (second_c+1) % 3));
+    }
+  };
+
+  sanity_check(uei);
+  sanity_check(ue_13);
+  sanity_check(ue_32);
+  sanity_check(ue_24);
+  sanity_check(ue_41);
+#endif
+}

include/igl/flip_edge.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingan Zhou <qnzhou@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_FLIP_EDGE_H
+#define IGL_FLIP_EDGE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Flip an edge in a triangle mesh.  The edge specified by uei must have
+  // exactly **two** adjacent faces.  Violation will result in exception.
+  // Another warning: edge flipping could convert manifold mesh into
+  // non-manifold.
+  //
+  // Inputs:
+  //   F    #F by 3 list of triangles.
+  //   E    #F*3 by 2 list of all of directed edges
+  //   uE   #uE by 2 list of unique undirected edges
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //        undirected edge
+  //   uE2E #uE list of lists of indices into E of coexisting edges
+  //   ue   index into uE the edge to be flipped.
+  //
+  // Output:
+  //   Updated F, E, uE, EMAP and uE2E.
+  template <
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType>
+  IGL_INLINE void flip_edge(
+    Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DeriveduE> & uE,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    std::vector<std::vector<uE2EType> > & uE2E,
+    const size_t uei);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "flip_edge.cpp"
+#endif
+#endif

include/igl/flipped_triangles.cpp

@@ -17,12 +17,12 @@ IGL_INLINE void igl::flipped_triangles(
 {
   assert(V.cols() == 2 && "V should contain 2D positions");
   std::vector<typename DerivedX::Scalar> flip_idx;
-  for (int i = 0; i < F.rows(); i++) 
+  for (int i = 0; i < F.rows(); i++)
   {
     // https://www.cs.cmu.edu/~quake/robust.html
     typedef Eigen::Matrix<typename DerivedV::Scalar,1,2> RowVector2S;
-    RowVector2S v1_n = V.row(F(i,0)); 
-    RowVector2S v2_n = V.row(F(i,1)); 
+    RowVector2S v1_n = V.row(F(i,0));
+    RowVector2S v2_n = V.row(F(i,1));
     RowVector2S v3_n = V.row(F(i,2));
     Eigen::Matrix<typename DerivedV::Scalar,3,3> T2_Homo;
     T2_Homo.col(0) << v1_n(0),v1_n(1),1.;
@@ -30,7 +30,7 @@ IGL_INLINE void igl::flipped_triangles(
     T2_Homo.col(2) << v3_n(0),v3_n(1),1.;
     double det = T2_Homo.determinant();
     assert(det == det && "det should not be NaN");
-    if (det < 0) 
+    if (det < 0)
     {
       flip_idx.push_back(i);
     }
@@ -51,4 +51,5 @@ IGL_INLINE Eigen::VectorXi igl::flipped_triangles(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::flipped_triangles<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template Eigen::Matrix<int, -1, 1, 0, -1, 1> igl::flipped_triangles<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&);
 #endif

include/igl/grad.cpp

@@ -9,12 +9,117 @@
 #include <Eigen/Geometry>
 #include <vector>
 
+#include "per_face_normals.h"
+#include "volume.h"
+#include "doublearea.h"
+
 template <typename DerivedV, typename DerivedF>
-IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
+IGL_INLINE void grad_tet(const Eigen::PlainObjectBase<DerivedV>&V,
+                     const Eigen::PlainObjectBase<DerivedF>&T,
+                            Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                            bool uniform) {
+  using namespace Eigen;
+  assert(T.cols() == 4);
+  const int n = V.rows(); int m = T.rows();
+
+  /*
+      F = [ ...
+      T(:,1) T(:,2) T(:,3); ...
+      T(:,1) T(:,3) T(:,4); ...
+      T(:,1) T(:,4) T(:,2); ...
+      T(:,2) T(:,4) T(:,3)]; */
+  MatrixXi F(4*m,3);
+  for (int i = 0; i < m; i++) {
+    F.row(0*m + i) << T(i,0), T(i,1), T(i,2);
+    F.row(1*m + i) << T(i,0), T(i,2), T(i,3);
+    F.row(2*m + i) << T(i,0), T(i,3), T(i,1);
+    F.row(3*m + i) << T(i,1), T(i,3), T(i,2);
+  }
+  // compute volume of each tet
+  VectorXd vol; igl::volume(V,T,vol);
+
+  VectorXd A(F.rows());
+  MatrixXd N(F.rows(),3);
+  if (!uniform) {
+    // compute tetrahedron face normals
+    igl::per_face_normals(V,F,N); int norm_rows = N.rows();
+    for (int i = 0; i < norm_rows; i++)
+      N.row(i) /= N.row(i).norm();
+    igl::doublearea(V,F,A); A/=2.;
+  } else {
+    // Use a uniform tetrahedra as a reference, with the same volume as the original one:
+    //
+    // Use normals of the uniform tet (V = h*[0,0,0;1,0,0;0.5,sqrt(3)/2.,0;0.5,sqrt(3)/6.,sqrt(2)/sqrt(3)])
+    //         0         0    1.0000
+    //         0.8165   -0.4714   -0.3333
+    //         0          0.9428   -0.3333
+    //         -0.8165   -0.4714   -0.3333
+    for (int i = 0; i < m; i++) {
+      N.row(0*m+i) << 0,0,1;
+      double a = sqrt(2)*std::cbrt(3*vol(i)); // area of a face in a uniform tet with volume = vol(i)
+      A(0*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(1*m+i) << 0.8165,-0.4714,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(1*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(2*m+i) << 0,0.9428,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(2*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+    for (int i = 0; i < m; i++) {
+      N.row(3*m+i) << -0.8165,-0.4714,-0.3333;
+      double a = sqrt(2)*std::cbrt(3*vol(i));
+      A(3*m+i) = (pow(a,2)*sqrt(3))/4.;
+    }
+
+  }
+
+  /*  G = sparse( ...
+      [0*m + repmat(1:m,1,4) ...
+       1*m + repmat(1:m,1,4) ...
+       2*m + repmat(1:m,1,4)], ...
+      repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1), ...
+      repmat(A./(3*repmat(vol,4,1)),3,1).*N(:), ...
+      3*m,n);*/
+  std::vector<Triplet<double> > G_t;
+  for (int i = 0; i < 4*m; i++) {
+    int T_j; // j indexes : repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1)
+    switch (i/m) {
+      case 0:
+        T_j = 3;
+        break;
+      case 1:
+        T_j = 1;
+        break;
+      case 2:
+        T_j = 2;
+        break;
+      case 3:
+        T_j = 0;
+        break;
+    }
+    int i_idx = i%m;
+    int j_idx = T(i_idx,T_j);
+
+    double val_before_n = A(i)/(3*vol(i_idx));
+    G_t.push_back(Triplet<double>(0*m+i_idx, j_idx, val_before_n * N(i,0)));
+    G_t.push_back(Triplet<double>(1*m+i_idx, j_idx, val_before_n * N(i,1)));
+    G_t.push_back(Triplet<double>(2*m+i_idx, j_idx, val_before_n * N(i,2)));
+  }
+  G.resize(3*m,n);
+  G.setFromTriplets(G_t.begin(), G_t.end());
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void grad_tri(const Eigen::PlainObjectBase<DerivedV>&V,
                      const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G)
+                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                    bool uniform)
 {
-  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> 
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3>
     eperp21(F.rows(),3), eperp13(F.rows(),3);
 
   for (int i=0;i<F.rows();++i)
@@ -28,17 +133,33 @@ IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
     Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v32 = V.row(i3) - V.row(i2);
     Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v13 = V.row(i1) - V.row(i3);
     Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v21 = V.row(i2) - V.row(i1);
-
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> n = v32.cross(v13);
     // area of parallelogram is twice area of triangle
     // area of parallelogram is || v1 x v2 ||
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> n  = v32.cross(v13);
-
     // This does correct l2 norm of rows, so that it contains #F list of twice
     // triangle areas
     double dblA = std::sqrt(n.dot(n));
+    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> u;
+    if (!uniform) {
+      // now normalize normals to get unit normals
+      u = n / dblA;
+    } else {
+      // Abstract equilateral triangle v1=(0,0), v2=(h,0), v3=(h/2, (sqrt(3)/2)*h)
+
+      // get h (by the area of the triangle)
+      double h = sqrt( (dblA)/sin(M_PI / 3.0)); // (h^2*sin(60))/2. = Area => h = sqrt(2*Area/sin_60)
+
+      Eigen::VectorXd v1,v2,v3;
+      v1 << 0,0,0;
+      v2 << h,0,0;
+      v3 << h/2.,(sqrt(3)/2.)*h,0;
 
-    // now normalize normals to get unit normals
-    Eigen::Matrix<typename DerivedV::Scalar, 1, 3> u = n / dblA;
+      // now fix v32,v13,v21 and the normal
+      v32 = v3-v2;
+      v13 = v1-v3;
+      v21 = v2-v1;
+      n = v32.cross(v13);
+    }
 
     // rotate each vector 90 degrees around normal
     double norm21 = std::sqrt(v21.dot(v21));
@@ -100,10 +221,21 @@ IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
   G.setFromTriplets(triplets.begin(), triplets.end());
 }
 
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
+                     const Eigen::PlainObjectBase<DerivedF>&F,
+                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                    bool uniform)
+{
+  assert(F.cols() == 3 || F.cols() == 4);
+  if (F.cols() == 3)
+    return grad_tri(V,F,G,uniform);
+  if (F.cols() == 4)
+    return grad_tet(V,F,G,uniform);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-// template void igl::grad<double, int>(Eigen::Matrix<double, -1, -1, 0, -1,-1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&,Eigen::SparseMatrix<double, 0, int>&);
-template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
-//template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&);
-template void igl::grad<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::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&);
+template void igl::grad<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::SparseMatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, 0, int>&, bool);
+template void igl::grad<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<Eigen::Matrix<double, -1, 3, 0, -1, 3>::Scalar, 0, int>&, bool);
 #endif

include/igl/grad.h

@@ -1,10 +1,10 @@
-// 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/.
+// 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_GRAD_MAT_H
 #define IGL_GRAD_MAT_H
 #include "igl_inline.h"
@@ -19,14 +19,15 @@ namespace igl {
   // Compute the numerical gradient operator
   //
   // Inputs:
-  //   V  #vertices by 3 list of mesh vertex positions
-  //   F  #faces by 3 list of mesh face indices
+  //   V          #vertices by 3 list of mesh vertex positions
+  //   F          #faces by 3 list of mesh face indices [or a #faces by 4 list of tetrahedral indices]
+  //   uniform    boolean (default false) - Use a uniform mesh instead of the vertices V
   // Outputs:
   //   G  #faces*dim by #V Gradient operator
   //
 
   // Gradient of a scalar function defined on piecewise linear elements (mesh)
-  // is constant on each triangle i,j,k:
+  // is constant on each triangle [tetrahedron] i,j,k:
   // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
   // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
   // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
@@ -35,10 +36,9 @@ namespace igl {
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void grad(const Eigen::PlainObjectBase<DerivedV>&V,
                      const Eigen::PlainObjectBase<DerivedF>&F,
-                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G);
-
+                    Eigen::SparseMatrix<typename DerivedV::Scalar> &G,
+                    bool uniform = false);
 }
-
 #ifndef IGL_STATIC_LIBRARY
 #  include "grad.cpp"
 #endif

include/igl/is_edge_manifold.cpp

@@ -52,11 +52,13 @@ IGL_INLINE bool igl::is_edge_manifold(
   return all;
 }
 
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedF>
 IGL_INLINE bool igl::is_edge_manifold(
-  const Eigen::PlainObjectBase<DerivedV>& /*V*/,
   const Eigen::PlainObjectBase<DerivedF>& F)
 {
+  // TODO: It's bothersome that this is not calling/reusing the code from the
+  // overload above. This could result in disagreement.
+
   // List of edges (i,j,f,c) where edge i<j is associated with corner i of face
   // f
   std::vector<std::vector<int> > TTT;
@@ -94,9 +96,9 @@ IGL_INLINE bool igl::is_edge_manifold(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template bool igl::is_edge_manifold<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
+template bool igl::is_edge_manifold<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
 // generated by autoexplicit.sh
 //template bool igl::is_edge_manifold<double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);
-template bool igl::is_edge_manifold<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&);
-template bool igl::is_edge_manifold<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
+template bool igl::is_edge_manifold<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
+template bool igl::is_edge_manifold<Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
 #endif

include/igl/is_edge_manifold.h

@@ -36,9 +36,8 @@ namespace igl
     Eigen::PlainObjectBase<DerivedE>& E,
     Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
     Eigen::PlainObjectBase<DerivedBE>& BE);
-  template <typename DerivedV, typename DerivedF>
+  template <typename DerivedF>
   IGL_INLINE bool is_edge_manifold(
-    const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F);
 }
 

include/igl/lexicographic_triangulation.cpp

@@ -0,0 +1,128 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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 "lexicographic_triangulation.h"
+#include "sortrows.h"
+
+#include <vector>
+#include <list>
+
+template<
+  typename DerivedP,
+  typename Orient2D,
+  typename DerivedF
+  >
+IGL_INLINE void igl::lexicographic_triangulation(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    Orient2D orient2D,
+    Eigen::PlainObjectBase<DerivedF>& F)
+{
+  typedef typename DerivedP::Scalar Scalar;
+  const size_t num_pts = P.rows();
+  if (num_pts < 3) {
+    throw "At least 3 points are required for triangulation!";
+  }
+
+  // Sort points in lexicographic order.
+  DerivedP ordered_P;
+  Eigen::VectorXi order;
+  igl::sortrows(P, true, ordered_P, order);
+
+  std::vector<Eigen::Vector3i> faces;
+  std::list<size_t> boundary;
+  const Scalar p0[] = {ordered_P(0, 0), ordered_P(0, 1)};
+  const Scalar p1[] = {ordered_P(1, 0), ordered_P(1, 1)};
+  for (size_t i=2; i<num_pts; i++) {
+    const Scalar curr_p[] = {ordered_P(i, 0), ordered_P(i, 1)};
+    if (faces.size() == 0) {
+      // All points processed so far are collinear.
+      // Check if the current point is collinear with every points before it.
+      auto orientation = orient2D(p0, p1, curr_p);
+      if (orientation != 0) {
+        // Add a fan of triangles eminating from curr_p.
+        if (orientation > 0) {
+          for (size_t j=0; j<=i-2; j++) {
+            faces.push_back({order[j], order[j+1], order[i]});
+          }
+        } else if (orientation < 0) {
+          for (size_t j=0; j<=i-2; j++) {
+            faces.push_back({order[j+1], order[j], order[i]});
+          }
+        }
+        // Initialize current boundary.
+        boundary.insert(boundary.end(), order.data(), order.data()+i+1);
+        if (orientation < 0) {
+          boundary.reverse();
+        }
+      }
+    } else {
+      const size_t bd_size = boundary.size();
+      assert(bd_size >= 3);
+      std::vector<short> orientations;
+      for (auto itr=boundary.begin(); itr!=boundary.end(); itr++) {
+        auto next_itr = std::next(itr, 1);
+        if (next_itr == boundary.end()) {
+          next_itr = boundary.begin();
+        }
+        const Scalar bd_p0[] = {P(*itr, 0), P(*itr, 1)};
+        const Scalar bd_p1[] = {P(*next_itr, 0), P(*next_itr, 1)};
+        auto orientation = orient2D(bd_p0, bd_p1, curr_p);
+        if (orientation < 0) {
+          faces.push_back({*next_itr, *itr, order[i]});
+        }
+        orientations.push_back(orientation);
+      }
+
+      auto left_itr = boundary.begin();
+      auto right_itr = boundary.begin();
+      auto curr_itr = boundary.begin();
+      for (size_t j=0; j<bd_size; j++, curr_itr++) {
+        size_t prev = (j+bd_size-1) % bd_size;
+        if (orientations[j] >= 0 && orientations[prev] < 0) {
+          right_itr = curr_itr;
+        } else if (orientations[j] < 0 && orientations[prev] >= 0) {
+          left_itr = curr_itr;
+        }
+      }
+      assert(left_itr != right_itr);
+
+      for (auto itr=left_itr; itr!=right_itr; itr++) {
+        if (itr == boundary.end()) itr = boundary.begin();
+        if (itr == right_itr) break;
+        if (itr == left_itr || itr == right_itr) continue;
+        itr = boundary.erase(itr);
+        if (itr == boundary.begin()) {
+            itr = boundary.end();
+        } else {
+            itr--;
+        }
+      }
+
+      if (right_itr == boundary.begin()) {
+        assert(std::next(left_itr, 1) == boundary.end());
+        boundary.insert(boundary.end(), order[i]);
+      } else {
+        assert(std::next(left_itr, 1) == right_itr);
+        boundary.insert(right_itr, order[i]);
+      }
+    }
+  }
+
+  const size_t num_faces = faces.size();
+  if (num_faces == 0) {
+      // All input points are collinear.
+      // Do nothing here.
+  } else {
+      F.resize(num_faces, 3);
+      for (size_t i=0; i<num_faces; i++) {
+          F.row(i) = faces[i];
+      }
+  }
+}
+

include/igl/lexicographic_triangulation.h

@@ -0,0 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//                    Qingan Zhou <qnzhou@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_LEXICOGRAPHIC_TRIANGULATION_H
+#define IGL_LEXICOGRAPHIC_TRIANGULATION_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Given a set of points in 2D, return a lexicographic triangulation of these
+  // points.
+  //
+  // Inputs:
+  //   P  #P by 2 list of vertex positions
+  //   orient2D  A functor such that orient2D(pa, pb, pc) returns
+  //               1 if pa,pb,pc forms a conterclockwise triangle.
+  //              -1 if pa,pb,pc forms a clockwise triangle.
+  //               0 if pa,pb,pc are collinear.
+  //              where the argument pa,pb,pc are of type Scalar[2].
+  //
+  // Outputs:
+  //   F  #F by 3 of faces in lexicographic triangulation.
+  template<
+    typename DerivedP,
+    typename Orient2D,
+    typename DerivedF
+    >
+  IGL_INLINE void lexicographic_triangulation(
+      const Eigen::PlainObjectBase<DerivedP>& P,
+      Orient2D orient2D,
+      Eigen::PlainObjectBase<DerivedF>& F);
+}
+
+
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lexicographic_triangulation.cpp"
+#endif
+#endif

include/igl/line_search.cpp

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 "line_search.h"
+
+IGL_INLINE double igl::line_search(
+  Eigen::MatrixXd& x,
+  const Eigen::MatrixXd& d,
+  double step_size,
+  std::function<double(Eigen::MatrixXd&)> energy,
+  double cur_energy)
+{
+  double old_energy;
+  if (cur_energy > 0) {
+    old_energy = cur_energy;
+  } else {
+    old_energy = energy(x); // no energy was given -> need to compute the current energy
+  }
+  double new_energy = old_energy;
+  int cur_iter = 0; int MAX_STEP_SIZE_ITER = 12;
+
+  while (new_energy >= old_energy && cur_iter < MAX_STEP_SIZE_ITER) {
+    Eigen::MatrixXd new_x = x + step_size * d;
+
+    double cur_e = energy(new_x);
+    if ( cur_e >= old_energy) {
+      step_size /= 2;
+    } else {
+      x = new_x;
+      new_energy = cur_e;
+    }
+    cur_iter++;
+  }
+  return new_energy;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/line_search.h

@@ -0,0 +1,42 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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_LINE_SEARCH_H
+#define IGL_LINE_SEARCH_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Implement a bisection linesearch to minimize a mesh-based energy on vertices given at 'x' at a search direction 'd',
+  // with initial step size. Stops when a point with lower energy is found, or after maximal iterations have been reached.
+  //
+  // Inputs:
+  //   x  						#X by dim list of variables
+  //   d  						#X by dim list of a given search direction
+  //   i_step_size  			initial step size
+  //   energy       			A function to compute the mesh-based energy (return an energy that is bigger than 0)
+  //   cur_energy(OPTIONAL)     The energy at the given point. Helps save redundant computations. 
+  //							This is optional. If not specified, the function will compute it.
+  // Outputs:
+  //		x  						#X by dim list of variables at the new location
+  // Returns the energy at the new point 'x'
+  IGL_INLINE double line_search(
+    Eigen::MatrixXd& x,
+    const Eigen::MatrixXd& d,
+    double i_step_size,
+    std::function<double(Eigen::MatrixXd&)> energy,
+    double cur_energy = -1);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "line_search.cpp"
+#endif
+
+#endif

include/igl/min_quad_with_fixed.cpp

@@ -90,7 +90,18 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   }
   // cache unknown followed by lagrange indices
   data.unknown_lagrange.resize(data.unknown.size()+data.lagrange.size());
-  data.unknown_lagrange << data.unknown, data.lagrange;
+  // Would like to do:
+  //data.unknown_lagrange << data.unknown, data.lagrange;
+  // but Eigen can't handle empty vectors in comma initialization
+  // https://forum.kde.org/viewtopic.php?f=74&t=107974&p=364947#p364947
+  if(data.unknown.size() > 0)
+  {
+    data.unknown_lagrange << data.unknown;
+  }
+  if(data.lagrange.size() > 0)
+  {
+    data.unknown_lagrange << data.lagrange;
+  }
 
   SparseMatrix<T> Auu;
   slice(A,data.unknown,data.unknown,Auu);
@@ -415,7 +426,19 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
     int neq = data.lagrange.size();
     // append lagrange multiplier rhs's
     VectorXT BBeq(B.size() + Beq.size());
-    BBeq << B, (Beq*-2.0);
+    // Would like to do:
+    // BBeq << B, (Beq*-2.0);
+    // but Eigen can't handle empty vectors in comma initialization
+    // https://forum.kde.org/viewtopic.php?f=74&t=107974&p=364947#p364947
+    if(B.size() > 0)
+    {
+      BBeq << B;
+    }
+    if(Beq.size() > 0)
+    {
+      BBeq << Beq*-2.;
+    }
+
     // Build right hand side
     VectorXT BBequl;
     igl::slice(BBeq,data.unknown_lagrange,BBequl);

include/igl/per_face_normals.cpp

@@ -121,4 +121,6 @@ template void igl::per_face_normals_stable<Eigen::Matrix<double, -1, -1, 0, -1,
 template void igl::per_face_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> >&);
 template void igl::per_face_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
 template void igl::per_face_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::per_face_normals<class Eigen::Matrix<double,-1,3,0,-1,3>,class Eigen::Matrix<int,-1,-1,0,-1,-1>,class Eigen::Matrix<double,-1,-1,0,-1,-1> >(class Eigen::PlainObjectBase<class Eigen::Matrix<double,-1,3,0,-1,3> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,-1,0,-1,-1> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<double,-1,-1,0,-1,-1> > &);
+
 #endif

include/igl/point_simplex_squared_distance.cpp

@@ -20,26 +20,29 @@ template <
   typename DerivedV,
   typename DerivedEle,
   typename Derivedsqr_d,
-  typename Derivedc>
+  typename Derivedc,
+  typename Derivedb>
 IGL_INLINE void igl::point_simplex_squared_distance(
   const Eigen::MatrixBase<Derivedp> & p,
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedEle> & Ele,
   const typename DerivedEle::Index primitive,
   Derivedsqr_d & sqr_d,
-  Eigen::PlainObjectBase<Derivedc> & c)
+  Eigen::PlainObjectBase<Derivedc> & c,
+  Eigen::PlainObjectBase<Derivedb> & bary)
 {
   typedef typename Derivedp::Scalar Scalar;
   typedef typename Eigen::Matrix<Scalar,1,DIM> Vector;
   typedef Vector Point;
+  typedef Eigen::PlainObjectBase<Derivedb> BaryPoint;
 
   const auto & Dot = [](const Point & a, const Point & b)->Scalar
   {
     return a.dot(b);
   };
   // Real-time collision detection, Ericson, Chapter 5
-  const auto & ClosestPtPointTriangle = 
-    [&Dot](Point p, Point a, Point b, Point c)->Point 
+  const auto & ClosestBaryPtPointTriangle = 
+    [&Dot](Point p, Point a, Point b, Point c, BaryPoint& bary_out )->Point 
   {
     // Check if P in vertex region outside A
     Vector ab = b - a;
@@ -47,42 +50,61 @@ IGL_INLINE void igl::point_simplex_squared_distance(
     Vector ap = p - a;
     Scalar d1 = Dot(ab, ap);
     Scalar d2 = Dot(ac, ap);
-    if (d1 <= 0.0 && d2 <= 0.0) return a; // barycentric coordinates (1,0,0)
+    if (d1 <= 0.0 && d2 <= 0.0) {
+      // barycentric coordinates (1,0,0)
+      bary_out << 1, 0, 0;
+      return a;
+    }
     // Check if P in vertex region outside B
     Vector bp = p - b;
     Scalar d3 = Dot(ab, bp);
     Scalar d4 = Dot(ac, bp);
-    if (d3 >= 0.0 && d4 <= d3) return b; // barycentric coordinates (0,1,0)
+    if (d3 >= 0.0 && d4 <= d3) {
+      // barycentric coordinates (0,1,0)
+      bary_out << 0, 1, 0;
+      return b;
+    }
     // Check if P in edge region of AB, if so return projection of P onto AB
     Scalar vc = d1*d4 - d3*d2;
     if( a != b)
     {
       if (vc <= 0.0 && d1 >= 0.0 && d3 <= 0.0) {
         Scalar v = d1 / (d1 - d3);
-        return a + v * ab; // barycentric coordinates (1-v,v,0)
+        // barycentric coordinates (1-v,v,0)
+        bary_out << 1-v, v, 0;
+        return a + v * ab;
       }
     }
     // Check if P in vertex region outside C
     Vector cp = p - c;
     Scalar d5 = Dot(ab, cp);
     Scalar d6 = Dot(ac, cp);
-    if (d6 >= 0.0 && d5 <= d6) return c; // barycentric coordinates (0,0,1)
+    if (d6 >= 0.0 && d5 <= d6) {
+      // barycentric coordinates (0,0,1)
+      bary_out << 0, 0, 1;
+      return c;
+    }
     // Check if P in edge region of AC, if so return projection of P onto AC
     Scalar vb = d5*d2 - d1*d6;
     if (vb <= 0.0 && d2 >= 0.0 && d6 <= 0.0) {
       Scalar w = d2 / (d2 - d6);
-      return a + w * ac; // barycentric coordinates (1-w,0,w)
+      // barycentric coordinates (1-w,0,w)
+      bary_out << 1-w, 0, w;
+      return a + w * ac;
     }
     // Check if P in edge region of BC, if so return projection of P onto BC
     Scalar va = d3*d6 - d5*d4;
     if (va <= 0.0 && (d4 - d3) >= 0.0 && (d5 - d6) >= 0.0) {
       Scalar w = (d4 - d3) / ((d4 - d3) + (d5 - d6));
-      return b + w * (c - b); // barycentric coordinates (0,1-w,w)
+      // barycentric coordinates (0,1-w,w)
+      bary_out << 0, 1-w, w;
+      return b + w * (c - b);
     }
     // P inside face region. Compute Q through its barycentric coordinates (u,v,w)
     Scalar denom = 1.0 / (va + vb + vc);
     Scalar v = vb * denom;
     Scalar w = vc * denom;
+    bary_out << 1.0-v-w, v, w;
     return a + ab * v + ac * w; // = u*a + v*b + w*c, u = va * denom = 1.0-v-w
   };
 
@@ -91,16 +113,41 @@ IGL_INLINE void igl::point_simplex_squared_distance(
   assert(Ele.cols() <= DIM+1);
   assert(Ele.cols() <= 3 && "Only simplices up to triangles are considered");
 
-  c = ClosestPtPointTriangle(
+  assert((Derivedb::RowsAtCompileTime == 1 || Derivedb::ColsAtCompileTime == 1) && "bary must be Eigen Vector or Eigen RowVector");
+  bary.resize( Derivedb::RowsAtCompileTime == 1 ? 1 : DIM, Derivedb::ColsAtCompileTime == 1 ? 1 : DIM );
+  c = ClosestBaryPtPointTriangle(
     p,
     V.row(Ele(primitive,0)),
     V.row(Ele(primitive,1%Ele.cols())),
-    V.row(Ele(primitive,2%Ele.cols())));
+    V.row(Ele(primitive,2%Ele.cols())),
+    bary);
   sqr_d = (p-c).squaredNorm();
 }
+template <
+  int DIM,
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedsqr_d,
+  typename Derivedc>
+IGL_INLINE void igl::point_simplex_squared_distance(
+  const Eigen::MatrixBase<Derivedp> & p,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedEle> & Ele,
+  const typename DerivedEle::Index primitive,
+  Derivedsqr_d & sqr_d,
+  Eigen::PlainObjectBase<Derivedc> & c)
+{
+    Eigen::PlainObjectBase<Derivedc> b;
+    point_simplex_squared_distance<DIM>( p, V, Ele, primitive, sqr_d, c, b );
+}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
+template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template void igl::point_simplex_squared_distance<3, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 1, 1, 3> >&);
+template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&);
+template void igl::point_simplex_squared_distance<2, Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, double&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 1, 1, 2> >&);
 #endif

include/igl/point_simplex_squared_distance.h

@@ -36,6 +36,35 @@ namespace igl
     const typename DerivedEle::Index i,
     Derivedsqr_d & sqr_d,
     Eigen::PlainObjectBase<Derivedc> & c);
+  // Determine squared distance from a point to linear simplex.
+  // Also return barycentric coordinate of closest point. 
+  //
+  // Inputs:
+  //   p  d-long query point
+  //   V  #V by d list of vertices
+  //   Ele  #Ele by ss<=d+1 list of simplex indices into V
+  //   i  index into Ele of simplex
+  // Outputs:
+  //   sqr_d  squared distance of Ele(i) to p
+  //   c  closest point on Ele(i) 
+  //   b  barycentric coordinates of closest point on Ele(i) 
+  //
+  template <
+    int DIM,
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedsqr_d,
+    typename Derivedc,
+    typename Derivedb>
+  IGL_INLINE void point_simplex_squared_distance(
+    const Eigen::MatrixBase<Derivedp> & p,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedEle> & Ele,
+    const typename DerivedEle::Index i,
+    Derivedsqr_d & sqr_d,
+    Eigen::PlainObjectBase<Derivedc> & c,
+    Eigen::PlainObjectBase<Derivedb> & b);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "point_simplex_squared_distance.cpp"

include/igl/readDMAT.cpp

@@ -213,6 +213,7 @@ IGL_INLINE bool igl::readDMAT(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
 template bool igl::readDMAT<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::string, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::readDMAT<double>(std::string, std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > >&);
 template bool igl::readDMAT<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::string, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/readSTL.h

@@ -23,9 +23,9 @@ namespace igl
   //   Scalar  type for positions and vectors (will be read as double and cast
   //     to Scalar)
   // Inputs:
-  //   filename path to .obj file
+  //   filename path to .stl file
   // Outputs:
-  //   V  double matrix of vertex positions  #F*3 by 3
+  //   V  double matrix of vertex positions  #V*3 by 3
   //   F  index matrix of triangle indices #F by 3
   //   N  double matrix of vertex positions  #F by 3
   // Returns true on success, false on errors

include/igl/remove_duplicates.cpp

@@ -84,4 +84,5 @@ IGL_INLINE void igl::remove_duplicates(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::remove_duplicates<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::Matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, -1, 1, 0, -1, 1>&, double);
+template void igl::remove_duplicates<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::Matrix<Eigen::Matrix<int, -1, 3, 1, -1, 3>::Scalar, -1, 1, 0, -1, 1>&, double);
 #endif

include/igl/seam_edges.cpp

@@ -11,23 +11,6 @@
 #include <unordered_set>
 #include <cassert>
 
-namespace {
-    // Computes the orientation of `c` relative to the line between `a` and `b`.
-    // Assumes 2D vector input.
-    // Based on: https://www.cs.cmu.edu/~quake/robust.html
-    template< typename scalar_t >
-    inline scalar_t orientation(
-        const Eigen::Matrix< scalar_t, 2, 1 >& a,
-        const Eigen::Matrix< scalar_t, 2, 1 >& b,
-        const Eigen::Matrix< scalar_t, 2, 1 >& c
-        ) {
-        typedef Eigen::Matrix< scalar_t, 2, 1 > Vector2S;
-        const Vector2S row0 = a - c;
-        const Vector2S row1 = b - c;
-        return row0(0)*row1(1) - row1(0)*row0(1);
-    }
-}
-
 // I have verified that this function produces the exact same output as
 // `find_seam_fast.py` for `cow_triangled.obj`.
 template <typename DerivedV, typename DerivedF, typename DerivedT>
@@ -49,6 +32,15 @@ IGL_INLINE void igl::seam_edges(
     // Assume 2D texture coordinates (foldovers tests).
     assert( TC.cols() == 2 );
     typedef Eigen::Matrix< typename DerivedT::Scalar, 2, 1 > Vector2S;
+    // Computes the orientation of `c` relative to the line between `a` and `b`.
+    // Assumes 2D vector input.
+    // Based on: https://www.cs.cmu.edu/~quake/robust.html
+    const auto& Orientation = []( const Vector2S& a, const Vector2S& b, const Vector2S& c ) -> typename DerivedT::Scalar
+    {
+        const Vector2S row0 = a - c;
+        const Vector2S row1 = b - c;
+        return row0(0)*row1(1) - row1(0)*row0(1);
+    };
     
     seams     .setZero( 3*F.rows(), 4 );
     boundaries.setZero( 3*F.rows(), 2 );
@@ -73,7 +65,7 @@ IGL_INLINE void igl::seam_edges(
     typedef std::pair< typename DerivedF::Scalar, typename DerivedF::Scalar > directed_edge;
 	const int numV = V.rows();
 	const int numF = F.rows();
-	auto edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
+	const auto& edge_hasher = [numV]( directed_edge const& e ) { return e.first*numV + e.second; };
 	// When we pass a hash function object, we also need to specify the number of buckets.
 	// The Euler characteristic says that the number of undirected edges is numV + numF -2*genus.
 	std::unordered_map< directed_edge, std::pair< int, int >, decltype( edge_hasher ) > directed_position_edge2face_position_index( 2*( numV + numF ), edge_hasher );
@@ -84,12 +76,13 @@ IGL_INLINE void igl::seam_edges(
         }
     }
     
-    // First find all undirected position edges (collect both orientations of
+    // First find all undirected position edges (collect a canonical orientation of
     // the directed edges).
     std::unordered_set< directed_edge, decltype( edge_hasher ) > undirected_position_edges( numV + numF, edge_hasher );
     for( const auto& el : directed_position_edge2face_position_index ) {
-        undirected_position_edges.insert( el.first );
-        undirected_position_edges.insert( std::make_pair( el.first.second, el.first.first ) );
+        // The canonical orientation is the one where the smaller of
+        // the two vertex indices is first.
+        undirected_position_edges.insert( std::make_pair( std::min( el.first.first, el.first.second ), std::max( el.first.first, el.first.second ) ) );
     }
     
     // Now we will iterate over all position edges.
@@ -97,6 +90,10 @@ IGL_INLINE void igl::seam_edges(
     // texcoord indices (or one doesn't exist at all in the case of a mesh
     // boundary).
     for( const auto& vp_edge : undirected_position_edges ) {
+        // We should only see canonical edges,
+        // where the first vertex index is smaller.
+        assert( vp_edge.first < vp_edge.second );
+        
         const auto vp_edge_reverse = std::make_pair( vp_edge.second, vp_edge.first );
         // If it and its opposite exist as directed edges, check if their
         // texture coordinate indices match.
@@ -108,12 +105,6 @@ IGL_INLINE void igl::seam_edges(
             // NOTE: They should never be equal.
             assert( forwards != backwards );
             
-            // NOTE: Non-matching seam edges will show up twice, once as
-            //       ( forwards, backwards ) and once as
-            //       ( backwards, forwards ). We don't need to examine both,
-            //       so continue only if forwards < backwards.
-            if( forwards < backwards ) continue;
-
             // If the texcoord indices match (are similarly flipped),
             // this edge is not a seam. It could be a foldover.
             if( std::make_pair( FTC( forwards.first, forwards.second ), FTC( forwards.first, ( forwards.second+1 ) % 3 ) )
@@ -129,8 +120,8 @@ IGL_INLINE void igl::seam_edges(
                 const Vector2S c_backwards = TC.row( FTC( backwards.first, (backwards.second+2) % 3 ) );
                 // If the opposite vertices' texture coordinates fall on the same side
                 // of the edge, we have a UV-space foldover.
-                const auto orientation_forwards = orientation( a, b, c_forwards );
-                const auto orientation_backwards = orientation( a, b, c_backwards );
+                const auto orientation_forwards = Orientation( a, b, c_forwards );
+                const auto orientation_backwards = Orientation( a, b, c_backwards );
                 if( ( orientation_forwards > 0 && orientation_backwards > 0 ) ||
                     ( orientation_forwards < 0 && orientation_backwards < 0 )
                     ) {

include/igl/slim.cpp

@@ -0,0 +1,755 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 "slim.h"
+
+#include <igl/boundary_loop.h>
+#include <igl/cotmatrix.h>
+#include <igl/edge_lengths.h>
+#include <igl/grad.h>
+#include <igl/local_basis.h>
+#include <igl/readOBJ.h>
+#include <igl/repdiag.h>
+#include <igl/vector_area_matrix.h>
+#include <iostream>
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <map>
+#include <set>
+#include <vector>
+
+#include "igl/arap.h"
+#include "igl/cat.h"
+#include "igl/doublearea.h"
+#include "igl/grad.h"
+#include "igl/local_basis.h"
+#include "igl/per_face_normals.h"
+#include "igl/slice_into.h"
+#include "igl/volume.h"
+#include "igl/polar_svd.h"
+
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/Sparse>
+#include <Eigen/SparseQR>
+#include <Eigen/SparseCholesky>
+#include <Eigen/IterativeLinearSolvers>
+
+#include "flip_avoiding_line_search.h"
+
+using namespace std;
+using namespace Eigen;
+
+///////// Helper functions to compute gradient matrices
+
+void compute_surface_gradient_matrix(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+                                     const Eigen::MatrixXd& F1, const Eigen::MatrixXd& F2,
+                                     Eigen::SparseMatrix<double>& D1, Eigen::SparseMatrix<double>& D2) {
+
+  Eigen::SparseMatrix<double> G;
+  igl::grad(V,F,G);
+  Eigen::SparseMatrix<double> Dx = G.block(0,0,F.rows(),V.rows());
+  Eigen::SparseMatrix<double> Dy = G.block(F.rows(),0,F.rows(),V.rows());
+  Eigen::SparseMatrix<double> Dz = G.block(2*F.rows(),0,F.rows(),V.rows());
+
+  D1 = F1.col(0).asDiagonal()*Dx + F1.col(1).asDiagonal()*Dy + F1.col(2).asDiagonal()*Dz;
+  D2 = F2.col(0).asDiagonal()*Dx + F2.col(1).asDiagonal()*Dy + F2.col(2).asDiagonal()*Dz;
+}
+
+// Computes the weights and solve the linear system for the quadratic proxy specified in the paper
+// The output of this is used to generate a search direction that will be fed to the Linesearch class
+class WeightedGlobalLocal {
+
+public:
+  WeightedGlobalLocal(igl::SLIMData& state);
+
+  // Compute necessary information before solving the proxy quadratic
+  void pre_calc();
+
+  // Solve the weighted proxy global step
+  // Output:
+  //    V_new #V by dim list of mesh positions (will be fed to a linesearch algorithm)
+  void solve_weighted_proxy(Eigen::MatrixXd& V_new);
+
+  // Compute the energy specified in the SLIMData structure + the soft constraint energy (in case there are soft constraints)
+  // Input:
+  //    V_new #V by dim list of mesh positions
+  virtual double compute_energy(Eigen::MatrixXd& V_new);
+
+//private:
+
+  void compute_jacobians(const Eigen::MatrixXd& V_o);
+  double compute_energy_with_jacobians(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& Ji, Eigen::MatrixXd& V_o, Eigen::VectorXd& areas);
+  double compute_soft_const_energy(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+                                             Eigen::MatrixXd& V_o);
+
+  void update_weights_and_closest_rotations(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXd& uv);
+  void solve_weighted_arap(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, Eigen::MatrixXd& uv, Eigen::VectorXi& b,
+      Eigen::MatrixXd& bc);
+
+  void build_linear_system(Eigen::SparseMatrix<double> &L);
+  void buildA(Eigen::SparseMatrix<double>& A);
+  void buildRhs(const Eigen::SparseMatrix<double>& At);
+
+  void add_soft_constraints(Eigen::SparseMatrix<double> &L);
+  void add_proximal_penalty();
+
+  igl::SLIMData& m_state;
+  Eigen::VectorXd M;
+  Eigen::VectorXd rhs;
+  Eigen::MatrixXd Ri,Ji;
+  Eigen::VectorXd W_11; Eigen::VectorXd W_12; Eigen::VectorXd W_13;
+  Eigen::VectorXd W_21; Eigen::VectorXd W_22; Eigen::VectorXd W_23;
+  Eigen::VectorXd W_31; Eigen::VectorXd W_32; Eigen::VectorXd W_33;
+  Eigen::SparseMatrix<double> Dx,Dy,Dz;
+
+  int f_n,v_n;
+
+  bool first_solve;
+  bool has_pre_calc = false;
+
+  int dim;
+};
+
+//// Implementation
+
+WeightedGlobalLocal::WeightedGlobalLocal(igl::SLIMData& state) :
+                                  m_state(state) {
+}
+
+void WeightedGlobalLocal::solve_weighted_proxy(Eigen::MatrixXd& V_new) {
+
+  update_weights_and_closest_rotations(m_state.V,m_state.F,V_new);
+  solve_weighted_arap(m_state.V,m_state.F,V_new,m_state.b,m_state.bc);
+}
+
+void WeightedGlobalLocal::compute_jacobians(const Eigen::MatrixXd& uv) {
+  if (m_state.F.cols() == 3){
+    // Ji=[D1*u,D2*u,D1*v,D2*v];
+    Ji.col(0) = Dx*uv.col(0); Ji.col(1) = Dy*uv.col(0);
+    Ji.col(2) = Dx*uv.col(1); Ji.col(3) = Dy*uv.col(1);
+  } else /*tet mesh*/{
+    // Ji=[D1*u,D2*u,D3*u, D1*v,D2*v, D3*v, D1*w,D2*w,D3*w];
+    Ji.col(0) = Dx*uv.col(0); Ji.col(1) = Dy*uv.col(0); Ji.col(2) = Dz*uv.col(0);
+    Ji.col(3) = Dx*uv.col(1); Ji.col(4) = Dy*uv.col(1); Ji.col(5) = Dz*uv.col(1);
+    Ji.col(6) = Dx*uv.col(2); Ji.col(7) = Dy*uv.col(2); Ji.col(8) = Dz*uv.col(2);
+  }
+}
+
+void WeightedGlobalLocal::update_weights_and_closest_rotations(const Eigen::MatrixXd& V,
+       const Eigen::MatrixXi& F, Eigen::MatrixXd& uv) {
+  compute_jacobians(uv);
+
+  const double eps = 1e-8;
+  double exp_f = m_state.exp_factor;
+
+  if (dim==2) {
+    for(int i=0; i <Ji.rows(); ++i ) {
+    typedef Eigen::Matrix<double,2,2> Mat2;
+    typedef Eigen::Matrix<double,2,1> Vec2;
+    Mat2 ji,ri,ti,ui,vi; Vec2 sing; Vec2 closest_sing_vec;Mat2 mat_W;
+    Vec2 m_sing_new;
+    double s1,s2;
+
+    ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1);
+    ji(1,0) = Ji(i,2); ji(1,1) = Ji(i,3);
+
+    igl::polar_svd(ji,ri,ti,ui,sing,vi);
+
+    s1 = sing(0); s2 = sing(1);
+
+    // Update Weights according to energy
+    switch(m_state.slim_energy) {
+    case igl::SLIMData::ARAP: {
+      m_sing_new << 1,1;
+      break;
+    } case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+        double s1_g = 2* (s1-pow(s1,-3));
+        double s2_g = 2 * (s2-pow(s2,-3));
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+    } case igl::SLIMData::LOG_ARAP: {
+        double s1_g = 2 * (log(s1)/s1);
+        double s2_g = 2 * (log(s2)/s2);
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+    } case igl::SLIMData::CONFORMAL: {
+        double s1_g = 1/(2*s2) - s2/(2*pow(s1,2));
+        double s2_g = 1/(2*s1) - s1/(2*pow(s2,2));
+
+        double geo_avg = sqrt(s1*s2);
+        double s1_min = geo_avg; double s2_min = geo_avg;
+
+        m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min)));
+
+        // change local step
+        closest_sing_vec << s1_min,s2_min;
+        ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
+        break;
+    } case igl::SLIMData::EXP_CONFORMAL: {
+        double s1_g = 2* (s1-pow(s1,-3));
+        double s2_g = 2 * (s2-pow(s2,-3));
+
+        double geo_avg = sqrt(s1*s2);
+        double s1_min = geo_avg; double s2_min = geo_avg;
+
+        double in_exp = exp_f*((pow(s1,2)+pow(s2,2))/(2*s1*s2));
+        double exp_thing = exp(in_exp);
+
+        s1_g *= exp_thing*exp_f;
+        s2_g *= exp_thing*exp_f;
+
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+    } case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+        double s1_g = 2* (s1-pow(s1,-3));
+        double s2_g = 2 * (s2-pow(s2,-3));
+
+        double in_exp = exp_f*(pow(s1,2)+pow(s1,-2)+pow(s2,2)+pow(s2,-2));
+        double exp_thing = exp(in_exp);
+
+        s1_g *= exp_thing*exp_f;
+        s2_g *= exp_thing*exp_f;
+
+        m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1)));
+        break;
+      }
+    }
+
+    if (abs(s1-1) < eps) m_sing_new(0) = 1; if (abs(s2-1) < eps) m_sing_new(1) = 1;
+    mat_W = ui*m_sing_new.asDiagonal()*ui.transpose();
+
+    W_11(i) = mat_W(0,0); W_12(i) = mat_W(0,1); W_21(i) = mat_W(1,0); W_22(i) = mat_W(1,1);
+
+    // 2) Update local step (doesn't have to be a rotation, for instance in case of conformal energy)
+    Ri(i,0) = ri(0,0); Ri(i,1) = ri(1,0); Ri(i,2) = ri(0,1); Ri(i,3) = ri(1,1);
+   }
+  } else {
+    typedef Eigen::Matrix<double,3,1> Vec3; typedef Eigen::Matrix<double,3,3> Mat3;
+    Mat3 ji; Vec3 m_sing_new; Vec3 closest_sing_vec;
+    const double sqrt_2 = sqrt(2);
+    for(int i=0; i <Ji.rows(); ++i ) {
+      ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1); ji(0,2) = Ji(i,2);
+      ji(1,0) = Ji(i,3); ji(1,1) = Ji(i,4); ji(1,2) = Ji(i,5);
+      ji(2,0) = Ji(i,6); ji(2,1) = Ji(i,7); ji(2,2) = Ji(i,8);
+
+      Mat3 ri,ti,ui,vi;
+      Vec3 sing;
+      igl::polar_svd(ji,ri,ti,ui,sing,vi);
+
+      double s1 = sing(0); double s2 = sing(1); double s3 = sing(2);
+
+      // 1) Update Weights
+      switch(m_state.slim_energy) {
+        case igl::SLIMData::ARAP: {
+          m_sing_new << 1,1,1;
+          break;
+        } case igl::SLIMData::LOG_ARAP: {
+            double s1_g = 2 * (log(s1)/s1);
+            double s2_g = 2 * (log(s2)/s2);
+            double s3_g = 2 * (log(s3)/s3);
+            m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+            break;
+          } case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+            double s1_g = 2* (s1-pow(s1,-3));
+            double s2_g = 2 * (s2-pow(s2,-3));
+            double s3_g = 2 * (s3-pow(s3,-3));
+            m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+            break;
+          } case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+           double s1_g = 2* (s1-pow(s1,-3));
+          double s2_g = 2 * (s2-pow(s2,-3));
+          double s3_g = 2 * (s3-pow(s3,-3));
+          m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+
+          double in_exp = exp_f*(pow(s1,2)+pow(s1,-2)+pow(s2,2)+pow(s2,-2)+pow(s3,2)+pow(s3,-2));
+          double exp_thing = exp(in_exp);
+
+          s1_g *= exp_thing*exp_f;
+          s2_g *= exp_thing*exp_f;
+          s3_g *= exp_thing*exp_f;
+
+          m_sing_new << sqrt(s1_g/(2*(s1-1))), sqrt(s2_g/(2*(s2-1))), sqrt(s3_g/(2*(s3-1)));
+
+          break;
+        }
+        case igl::SLIMData::CONFORMAL: {
+          double common_div = 9*(pow(s1*s2*s3,5./3.));
+
+          double s1_g = (-2*s2*s3*(pow(s2,2)+pow(s3,2)-2*pow(s1,2)) ) / common_div;
+          double s2_g = (-2*s1*s3*(pow(s1,2)+pow(s3,2)-2*pow(s2,2)) ) / common_div;
+          double s3_g = (-2*s1*s2*(pow(s1,2)+pow(s2,2)-2*pow(s3,2)) ) / common_div;
+
+          double closest_s = sqrt(pow(s1,2)+pow(s3,2)) / sqrt_2;
+          double s1_min = closest_s; double s2_min = closest_s; double s3_min = closest_s;
+
+          m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min))), sqrt(s3_g/(2*(s3-s3_min)));
+
+          // change local step
+          closest_sing_vec << s1_min,s2_min,s3_min;
+          ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
+          break;
+        }
+        case igl::SLIMData::EXP_CONFORMAL: {
+          // E_conf = (s1^2 + s2^2 + s3^2)/(3*(s1*s2*s3)^(2/3) )
+          // dE_conf/ds1 = (-2*(s2*s3)*(s2^2+s3^2 -2*s1^2) ) / (9*(s1*s2*s3)^(5/3))
+          // Argmin E_conf(s1): s1 = sqrt(s1^2+s2^2)/sqrt(2)
+          double common_div = 9*(pow(s1*s2*s3,5./3.));
+
+          double s1_g = (-2*s2*s3*(pow(s2,2)+pow(s3,2)-2*pow(s1,2)) ) / common_div;
+          double s2_g = (-2*s1*s3*(pow(s1,2)+pow(s3,2)-2*pow(s2,2)) ) / common_div;
+          double s3_g = (-2*s1*s2*(pow(s1,2)+pow(s2,2)-2*pow(s3,2)) ) / common_div;
+
+          double in_exp = exp_f*( (pow(s1,2)+pow(s2,2)+pow(s3,2))/ (3*pow((s1*s2*s3),2./3)) ); ;
+          double exp_thing = exp(in_exp);
+
+          double closest_s = sqrt(pow(s1,2)+pow(s3,2)) / sqrt_2;
+          double s1_min = closest_s; double s2_min = closest_s; double s3_min = closest_s;
+
+          s1_g *= exp_thing*exp_f;
+          s2_g *= exp_thing*exp_f;
+          s3_g *= exp_thing*exp_f;
+
+          m_sing_new << sqrt(s1_g/(2*(s1-s1_min))), sqrt(s2_g/(2*(s2-s2_min))), sqrt(s3_g/(2*(s3-s3_min)));
+
+          // change local step
+          closest_sing_vec << s1_min,s2_min,s3_min;
+          ri = ui*closest_sing_vec.asDiagonal()*vi.transpose();
+        }
+      }
+      if (abs(s1-1) < eps) m_sing_new(0) = 1; if (abs(s2-1) < eps) m_sing_new(1) = 1; if (abs(s3-1) < eps) m_sing_new(2) = 1;
+      Mat3 mat_W;
+      mat_W = ui*m_sing_new.asDiagonal()*ui.transpose();
+
+      W_11(i) = mat_W(0,0);
+      W_12(i) = mat_W(0,1);
+      W_13(i) = mat_W(0,2);
+      W_21(i) = mat_W(1,0);
+      W_22(i) = mat_W(1,1);
+      W_23(i) = mat_W(1,2);
+      W_31(i) = mat_W(2,0);
+      W_32(i) = mat_W(2,1);
+      W_33(i) = mat_W(2,2);
+
+      // 2) Update closest rotations (not rotations in case of conformal energy)
+      Ri(i,0) = ri(0,0); Ri(i,1) = ri(1,0); Ri(i,2) = ri(2,0);
+      Ri(i,3) = ri(0,1); Ri(i,4) = ri(1,1); Ri(i,5) = ri(2,1);
+      Ri(i,6) = ri(0,2); Ri(i,7) = ri(1,2); Ri(i,8) = ri(2,2);
+    } // for loop end
+
+  } // if dim end
+
+}
+
+void WeightedGlobalLocal::solve_weighted_arap(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F,
+        Eigen::MatrixXd& uv, Eigen::VectorXi& soft_b_p, Eigen::MatrixXd& soft_bc_p) {
+  using namespace Eigen;
+
+  Eigen::SparseMatrix<double> L;
+  build_linear_system(L);
+
+  // solve
+  Eigen::VectorXd Uc;
+  if (dim == 2) {
+    SimplicialLDLT<SparseMatrix<double> > solver;
+    Uc = solver.compute(L).solve(rhs);
+  } else { // seems like CG performs much worse for 2D and way better for 3D
+    Eigen::VectorXd guess(uv.rows()*dim);
+    for (int i = 0; i < dim; i++) for (int j = 0; j < dim; j++) guess(uv.rows()*i + j) = uv(i,j); // flatten vector
+    ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
+    solver.setTolerance(1e-8);
+    Uc = solver.compute(L).solveWithGuess(rhs,guess);
+  }
+
+  for (int i = 0; i < dim; i++)
+    uv.col(i) = Uc.block(i*v_n,0,v_n,1);
+}
+
+void WeightedGlobalLocal::pre_calc() {
+  if (!has_pre_calc) {
+    v_n = m_state.v_num; f_n = m_state.f_num;
+
+    if (m_state.F.cols() == 3) {
+      dim = 2;
+      Eigen::MatrixXd F1,F2,F3;
+      igl::local_basis(m_state.V,m_state.F,F1,F2,F3);
+      compute_surface_gradient_matrix(m_state.V,m_state.F,F1,F2,Dx,Dy);
+
+      W_11.resize(f_n); W_12.resize(f_n); W_21.resize(f_n); W_22.resize(f_n);
+    } else {
+      dim = 3;
+      Eigen::SparseMatrix<double> G;
+      igl::grad(m_state.V,m_state.F,G,m_state.mesh_improvement_3d /*use normal gradient, or one from a "regular" tet*/);
+      Dx = G.block(0,0,m_state.F.rows(),m_state.V.rows());
+      Dy = G.block(m_state.F.rows(),0,m_state.F.rows(),m_state.V.rows());
+      Dz = G.block(2*m_state.F.rows(),0,m_state.F.rows(),m_state.V.rows());
+
+
+      W_11.resize(f_n);W_12.resize(f_n);W_13.resize(f_n);
+      W_21.resize(f_n);W_22.resize(f_n);W_23.resize(f_n);
+      W_31.resize(f_n);W_32.resize(f_n);W_33.resize(f_n);
+    }
+
+    Dx.makeCompressed();Dy.makeCompressed(); Dz.makeCompressed();
+    Ri.resize(f_n, dim*dim); Ji.resize(f_n, dim*dim);
+    rhs.resize(dim*m_state.v_num);
+
+    // flattened weight matrix
+    M.resize(dim*dim*f_n);
+    for (int i = 0; i < dim*dim; i++)
+      for (int j = 0; j < f_n; j++)
+        M(i*f_n + j) = m_state.M(j);
+
+    first_solve = true;
+    has_pre_calc = true;
+  }
+}
+
+void WeightedGlobalLocal::build_linear_system(Eigen::SparseMatrix<double> &L) {
+  // formula (35) in paper
+  Eigen::SparseMatrix<double> A(dim*dim*f_n, dim*v_n);
+  buildA(A);
+
+  Eigen::SparseMatrix<double> At = A.transpose();
+  At.makeCompressed();
+
+  Eigen::SparseMatrix<double> id_m(At.rows(),At.rows()); id_m.setIdentity();
+
+  // add proximal penalty
+  L = At*M.asDiagonal()*A + m_state.proximal_p * id_m; //add also a proximal term
+  L.makeCompressed();
+
+  buildRhs(At);
+  Eigen::SparseMatrix<double> OldL = L;
+  add_soft_constraints(L);
+  L.makeCompressed();
+}
+
+void WeightedGlobalLocal::add_soft_constraints(Eigen::SparseMatrix<double> &L) {
+  int v_n = m_state.v_num;
+  for (int d = 0; d < dim; d++) {
+    for (int i = 0; i < m_state.b.rows(); i++) {
+      int v_idx = m_state.b(i);
+      rhs(d*v_n + v_idx) += m_state.soft_const_p * m_state.bc(i,d); // rhs
+      L.coeffRef(d*v_n + v_idx, d*v_n + v_idx) += m_state.soft_const_p; // diagonal of matrix
+    }
+  }
+}
+
+double WeightedGlobalLocal::compute_energy(Eigen::MatrixXd& V_new) {
+  compute_jacobians(V_new);
+  return compute_energy_with_jacobians(m_state.V,m_state.F, Ji, V_new,m_state.M) + compute_soft_const_energy(m_state.V,m_state.F,V_new);
+}
+
+double WeightedGlobalLocal::compute_soft_const_energy(const Eigen::MatrixXd& V,
+                                                       const Eigen::MatrixXi& F,
+                                                       Eigen::MatrixXd& V_o) {
+  double e = 0;
+  for (int i = 0; i < m_state.b.rows(); i++) {
+    e += m_state.soft_const_p*(m_state.bc.row(i)-V_o.row(m_state.b(i))).squaredNorm();
+  }
+  return e;
+}
+
+double WeightedGlobalLocal::compute_energy_with_jacobians(const Eigen::MatrixXd& V,
+       const Eigen::MatrixXi& F, const Eigen::MatrixXd& Ji, Eigen::MatrixXd& uv, Eigen::VectorXd& areas) {
+
+  double energy = 0;
+  if (dim == 2) {
+    Eigen::Matrix<double,2,2> ji;
+    for (int i = 0; i < f_n; i++) {
+      ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1);
+      ji(1,0) = Ji(i,2); ji(1,1) = Ji(i,3);
+
+      typedef Eigen::Matrix<double,2,2> Mat2;
+      typedef Eigen::Matrix<double,2,1> Vec2;
+      Mat2 ri,ti,ui,vi; Vec2 sing;
+      igl::polar_svd(ji,ri,ti,ui,sing,vi);
+      double s1 = sing(0); double s2 = sing(1);
+
+      switch(m_state.slim_energy) {
+        case igl::SLIMData::ARAP: {
+          energy+= areas(i) * (pow(s1-1,2) + pow(s2-1,2));
+          break;
+        }
+        case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * (pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2));
+          break;
+        }
+        case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * exp(m_state.exp_factor*(pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2)));
+          break;
+        }
+        case igl::SLIMData::LOG_ARAP: {
+          energy += areas(i) * (pow(log(s1),2) + pow(log(s2),2));
+          break;
+        }
+        case igl::SLIMData::CONFORMAL: {
+          energy += areas(i) * ( (pow(s1,2)+pow(s2,2))/(2*s1*s2) );
+          break;
+        }
+        case igl::SLIMData::EXP_CONFORMAL: {
+          energy += areas(i) * exp(m_state.exp_factor*((pow(s1,2)+pow(s2,2))/(2*s1*s2)));
+          break;
+        }
+
+      }
+
+    }
+  } else {
+    Eigen::Matrix<double,3,3> ji;
+    for (int i = 0; i < f_n; i++) {
+      ji(0,0) = Ji(i,0); ji(0,1) = Ji(i,1); ji(0,2) = Ji(i,2);
+      ji(1,0) = Ji(i,3); ji(1,1) = Ji(i,4); ji(1,2) = Ji(i,5);
+      ji(2,0) = Ji(i,6); ji(2,1) = Ji(i,7); ji(2,2) = Ji(i,8);
+
+      typedef Eigen::Matrix<double,3,3> Mat3;
+      typedef Eigen::Matrix<double,3,1> Vec3;
+      Mat3 ri,ti,ui,vi; Vec3 sing;
+      igl::polar_svd(ji,ri,ti,ui,sing,vi);
+      double s1 = sing(0); double s2 = sing(1); double s3 = sing(2);
+
+      switch(m_state.slim_energy) {
+        case igl::SLIMData::ARAP: {
+          energy+= areas(i) * (pow(s1-1,2) + pow(s2-1,2) + pow(s3-1,2));
+          break;
+        }
+        case igl::SLIMData::SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * (pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2) + pow(s3,2) + pow(s3,-2));
+          break;
+        }
+        case igl::SLIMData::EXP_SYMMETRIC_DIRICHLET: {
+          energy += areas(i) * exp(m_state.exp_factor*(pow(s1,2) +pow(s1,-2) + pow(s2,2) + pow(s2,-2) + pow(s3,2) + pow(s3,-2)));
+          break;
+        }
+        case igl::SLIMData::LOG_ARAP: {
+          energy += areas(i) * (pow(log(s1),2) + pow(log(abs(s2)),2) + pow(log(abs(s3)),2));
+          break;
+        }
+        case igl::SLIMData::CONFORMAL: {
+          energy += areas(i) * ( ( pow(s1,2)+pow(s2,2)+pow(s3,2) ) /(3*pow(s1*s2*s3,2./3.)) );
+          break;
+        }
+        case igl::SLIMData::EXP_CONFORMAL: {
+          energy += areas(i) * exp( ( pow(s1,2)+pow(s2,2)+pow(s3,2) ) /(3*pow(s1*s2*s3,2./3.)) );
+          break;
+        }
+      }
+    }
+  }
+
+  return energy;
+}
+
+void WeightedGlobalLocal::buildA(Eigen::SparseMatrix<double>& A) {
+  // formula (35) in paper
+  std::vector<Triplet<double> > IJV;
+  if (dim == 2) {
+    IJV.reserve(4*(Dx.outerSize()+ Dy.outerSize()));
+
+    /*A = [W11*Dx, W12*Dx;
+         W11*Dy, W12*Dy;
+         W21*Dx, W22*Dx;
+         W21*Dy, W22*Dy];*/
+    for (int k=0; k<Dx.outerSize(); ++k) {
+        for (SparseMatrix<double>::InnerIterator it(Dx,k); it; ++it) {
+          int dx_r = it.row();
+          int dx_c = it.col();
+          double val = it.value();
+
+          IJV.push_back(Triplet<double>(dx_r,dx_c, val*W_11(dx_r)));
+          IJV.push_back(Triplet<double>(dx_r,v_n + dx_c, val*W_12(dx_r)));
+
+          IJV.push_back(Triplet<double>(2*f_n+dx_r,dx_c, val*W_21(dx_r)));
+          IJV.push_back(Triplet<double>(2*f_n+dx_r,v_n + dx_c, val*W_22(dx_r)));
+      }
+    }
+
+    for (int k=0; k<Dy.outerSize(); ++k) {
+      for (SparseMatrix<double>::InnerIterator it(Dy,k); it; ++it) {
+        int dy_r = it.row();
+        int dy_c = it.col();
+        double val = it.value();
+
+        IJV.push_back(Triplet<double>(f_n+dy_r,dy_c, val*W_11(dy_r)));
+        IJV.push_back(Triplet<double>(f_n+dy_r,v_n + dy_c, val*W_12(dy_r)));
+
+        IJV.push_back(Triplet<double>(3*f_n+dy_r,dy_c, val*W_21(dy_r)));
+        IJV.push_back(Triplet<double>(3*f_n+dy_r,v_n + dy_c, val*W_22(dy_r)));
+      }
+    }
+  } else {
+
+    /*A = [W11*Dx, W12*Dx, W13*Dx;
+           W11*Dy, W12*Dy, W13*Dy;
+           W11*Dz, W12*Dz, W13*Dz;
+           W21*Dx, W22*Dx, W23*Dx;
+           W21*Dy, W22*Dy, W23*Dy;
+           W21*Dz, W22*Dz, W23*Dz;
+           W31*Dx, W32*Dx, W33*Dx;
+           W31*Dy, W32*Dy, W33*Dy;
+           W31*Dz, W32*Dz, W33*Dz;];*/
+    IJV.reserve(9*(Dx.outerSize()+ Dy.outerSize() + Dz.outerSize()));
+    for (int k = 0; k < Dx.outerSize(); k++) {
+      for (SparseMatrix<double>::InnerIterator it(Dx,k); it; ++it) {
+         int dx_r = it.row();
+         int dx_c = it.col();
+         double val = it.value();
+
+         IJV.push_back(Triplet<double>(dx_r,dx_c, val*W_11(dx_r)));
+         IJV.push_back(Triplet<double>(dx_r,v_n + dx_c, val*W_12(dx_r)));
+         IJV.push_back(Triplet<double>(dx_r,2*v_n + dx_c, val*W_13(dx_r)));
+
+         IJV.push_back(Triplet<double>(3*f_n+dx_r,dx_c, val*W_21(dx_r)));
+         IJV.push_back(Triplet<double>(3*f_n+dx_r,v_n + dx_c, val*W_22(dx_r)));
+         IJV.push_back(Triplet<double>(3*f_n+dx_r,2*v_n + dx_c, val*W_23(dx_r)));
+
+         IJV.push_back(Triplet<double>(6*f_n+dx_r,dx_c, val*W_31(dx_r)));
+         IJV.push_back(Triplet<double>(6*f_n+dx_r,v_n + dx_c, val*W_32(dx_r)));
+         IJV.push_back(Triplet<double>(6*f_n+dx_r,2*v_n + dx_c, val*W_33(dx_r)));
+      }
+    }
+
+    for (int k = 0; k < Dy.outerSize(); k++) {
+      for (SparseMatrix<double>::InnerIterator it(Dy,k); it; ++it) {
+         int dy_r = it.row();
+         int dy_c = it.col();
+         double val = it.value();
+
+         IJV.push_back(Triplet<double>(f_n+dy_r,dy_c, val*W_11(dy_r)));
+         IJV.push_back(Triplet<double>(f_n+dy_r,v_n + dy_c, val*W_12(dy_r)));
+         IJV.push_back(Triplet<double>(f_n+dy_r,2*v_n + dy_c, val*W_13(dy_r)));
+
+         IJV.push_back(Triplet<double>(4*f_n+dy_r,dy_c, val*W_21(dy_r)));
+         IJV.push_back(Triplet<double>(4*f_n+dy_r,v_n + dy_c, val*W_22(dy_r)));
+         IJV.push_back(Triplet<double>(4*f_n+dy_r,2*v_n + dy_c, val*W_23(dy_r)));
+
+         IJV.push_back(Triplet<double>(7*f_n+dy_r,dy_c, val*W_31(dy_r)));
+         IJV.push_back(Triplet<double>(7*f_n+dy_r,v_n + dy_c, val*W_32(dy_r)));
+         IJV.push_back(Triplet<double>(7*f_n+dy_r,2*v_n + dy_c, val*W_33(dy_r)));
+      }
+    }
+
+    for (int k = 0; k < Dz.outerSize(); k++) {
+      for (SparseMatrix<double>::InnerIterator it(Dz,k); it; ++it) {
+         int dz_r = it.row();
+         int dz_c = it.col();
+         double val = it.value();
+
+         IJV.push_back(Triplet<double>(2*f_n + dz_r,dz_c, val*W_11(dz_r)));
+         IJV.push_back(Triplet<double>(2*f_n + dz_r,v_n + dz_c, val*W_12(dz_r)));
+         IJV.push_back(Triplet<double>(2*f_n + dz_r,2*v_n + dz_c, val*W_13(dz_r)));
+
+         IJV.push_back(Triplet<double>(5*f_n+dz_r,dz_c, val*W_21(dz_r)));
+         IJV.push_back(Triplet<double>(5*f_n+dz_r,v_n + dz_c, val*W_22(dz_r)));
+         IJV.push_back(Triplet<double>(5*f_n+dz_r,2*v_n + dz_c, val*W_23(dz_r)));
+
+         IJV.push_back(Triplet<double>(8*f_n+dz_r,dz_c, val*W_31(dz_r)));
+         IJV.push_back(Triplet<double>(8*f_n+dz_r,v_n + dz_c, val*W_32(dz_r)));
+         IJV.push_back(Triplet<double>(8*f_n+dz_r,2*v_n + dz_c, val*W_33(dz_r)));
+      }
+    }
+  }
+  A.setFromTriplets(IJV.begin(),IJV.end());
+}
+
+void WeightedGlobalLocal::buildRhs(const Eigen::SparseMatrix<double>& At) {
+  VectorXd f_rhs(dim*dim*f_n); f_rhs.setZero();
+  if (dim==2) {
+    /*b = [W11*R11 + W12*R21; (formula (36))
+         W11*R12 + W12*R22;
+         W21*R11 + W22*R21;
+         W21*R12 + W22*R22];*/
+    for (int i = 0; i < f_n; i++) {
+      f_rhs(i+0*f_n) = W_11(i) * Ri(i,0) + W_12(i)*Ri(i,1);
+      f_rhs(i+1*f_n) = W_11(i) * Ri(i,2) + W_12(i)*Ri(i,3);
+      f_rhs(i+2*f_n) = W_21(i) * Ri(i,0) + W_22(i)*Ri(i,1);
+      f_rhs(i+3*f_n) = W_21(i) * Ri(i,2) + W_22(i)*Ri(i,3);
+    }
+  } else {
+    /*b = [W11*R11 + W12*R21 + W13*R31;
+         W11*R12 + W12*R22 + W13*R32;
+         W11*R13 + W12*R23 + W13*R33;
+         W21*R11 + W22*R21 + W23*R31;
+         W21*R12 + W22*R22 + W23*R32;
+         W21*R13 + W22*R23 + W23*R33;
+         W31*R11 + W32*R21 + W33*R31;
+         W31*R12 + W32*R22 + W33*R32;
+         W31*R13 + W32*R23 + W33*R33;];*/
+    for (int i = 0; i < f_n; i++) {
+      f_rhs(i+0*f_n) = W_11(i) * Ri(i,0) + W_12(i)*Ri(i,1) + W_13(i)*Ri(i,2);
+      f_rhs(i+1*f_n) = W_11(i) * Ri(i,3) + W_12(i)*Ri(i,4) + W_13(i)*Ri(i,5);
+      f_rhs(i+2*f_n) = W_11(i) * Ri(i,6) + W_12(i)*Ri(i,7) + W_13(i)*Ri(i,8);
+      f_rhs(i+3*f_n) = W_21(i) * Ri(i,0) + W_22(i)*Ri(i,1) + W_23(i)*Ri(i,2);
+      f_rhs(i+4*f_n) = W_21(i) * Ri(i,3) + W_22(i)*Ri(i,4) + W_23(i)*Ri(i,5);
+      f_rhs(i+5*f_n) = W_21(i) * Ri(i,6) + W_22(i)*Ri(i,7) + W_23(i)*Ri(i,8);
+      f_rhs(i+6*f_n) = W_31(i) * Ri(i,0) + W_32(i)*Ri(i,1) + W_33(i)*Ri(i,2);
+      f_rhs(i+7*f_n) = W_31(i) * Ri(i,3) + W_32(i)*Ri(i,4) + W_33(i)*Ri(i,5);
+      f_rhs(i+8*f_n) = W_31(i) * Ri(i,6) + W_32(i)*Ri(i,7) + W_33(i)*Ri(i,8);
+    }
+  }
+  VectorXd uv_flat(dim*v_n);
+  for (int i = 0; i < dim; i++)
+    for (int j = 0; j < v_n; j++)
+      uv_flat(v_n*i+j) = m_state.V_o(j,i);
+
+  rhs = (At*M.asDiagonal()*f_rhs + m_state.proximal_p * uv_flat);
+}
+
+
+
+// #define TwoPi  6.28318530717958648
+// const double eps=1e-14;
+
+/// Slim Implementation
+
+IGL_INLINE void igl::slim_precompute(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::MatrixXd& V_init, SLIMData& data,
+   SLIMData::SLIM_ENERGY slim_energy, Eigen::VectorXi& b, Eigen::MatrixXd& bc, double soft_p) {
+
+  data.V = V;
+  data.F = F;
+  data.V_o = V_init;
+
+  data.v_num = V.rows();
+  data.f_num = F.rows();
+
+  data.slim_energy = slim_energy;
+
+  data.b = b;
+  data.bc = bc;
+  data.soft_const_p = soft_p;
+
+  data.proximal_p = 0.0001;
+
+  igl::doublearea(V,F,data.M); data.M /= 2.;
+  data.mesh_area = data.M.sum();
+  data.mesh_improvement_3d = false; // whether to use a jacobian derived from a real mesh or an abstract regular mesh (used for mesh improvement)
+  data.exp_factor = 1.0; // param used only for exponential energies (e.g exponential symmetric dirichlet)
+
+  assert (F.cols() == 3 || F.cols() == 4);
+  data.wGlobalLocal = new WeightedGlobalLocal(data);
+
+  data.wGlobalLocal->pre_calc();
+  data.energy = data.wGlobalLocal->compute_energy(data.V_o)/data.mesh_area;
+}
+
+IGL_INLINE void igl::slim_solve(SLIMData& data, int iter_num) {
+
+  for (int i = 0; i < iter_num; i++) {
+    Eigen::MatrixXd dest_res;
+    dest_res = data.V_o;
+    data.wGlobalLocal->solve_weighted_proxy(dest_res);
+
+    double old_energy = data.energy;
+
+    std::function<double(Eigen::MatrixXd&)> compute_energy = [&](Eigen::MatrixXd& aaa) { return data.wGlobalLocal->compute_energy(aaa); };
+
+    data.energy = igl::flip_avoiding_line_search(data.F,data.V_o, dest_res, compute_energy,
+                                           data.energy*data.mesh_area)/data.mesh_area;
+  }
+}

include/igl/slim.h

@@ -0,0 +1,89 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Michael Rabinovich
+//
+// 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 SLIM_H
+#define SLIM_H
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include <string>
+
+#include <igl/jet.h>
+#include <igl/readOBJ.h>
+#include <igl/facet_components.h>
+#include <igl/slice.h>
+
+class WeightedGlobalLocal;
+
+namespace igl
+{
+
+// Compute a SLIM map as derived in "Scalable Locally Injective Maps" [Rabinovich et al. 2016].
+struct SLIMData {
+
+  // Input
+  Eigen::MatrixXd V; // #V by 3 list of mesh vertex positions
+  Eigen::MatrixXi F; // #F by 3/3 list of mesh faces (triangles/tets)
+  enum SLIM_ENERGY {
+    ARAP,
+    LOG_ARAP,
+    SYMMETRIC_DIRICHLET,
+    CONFORMAL,
+    EXP_CONFORMAL,
+    EXP_SYMMETRIC_DIRICHLET
+  };
+  SLIM_ENERGY slim_energy;
+
+  // Optional Input
+    // soft constraints
+    Eigen::VectorXi b;
+    Eigen::MatrixXd bc;
+    double soft_const_p;
+
+  double exp_factor; // used for exponential energies, ignored otherwise
+  bool mesh_improvement_3d; // only supported for 3d
+
+  // Output
+  Eigen::MatrixXd V_o; // #V by dim list of mesh vertex positions (dim = 2 for parametrization, 3 otherwise)
+  double energy; // objective value
+
+  // INTERNAL
+  Eigen::VectorXd M;
+  double mesh_area;
+  double avg_edge_length;
+  int v_num;
+  int f_num;
+  double proximal_p;
+
+  WeightedGlobalLocal* wGlobalLocal;
+};
+
+  // Compute necessary information to start using SLIM
+  // Inputs:
+  //		V           #V by 3 list of mesh vertex positions
+  //		F           #F by 3/3 list of mesh faces (triangles/tets)
+  //    b           list of boundary indices into V
+  //    bc          #b by dim list of boundary conditions
+  //    soft_p      Soft penalty factor (can be zero)
+  //    slim_energy Energy to minimize
+IGL_INLINE void slim_precompute(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::MatrixXd& V_init, SLIMData& data,
+   SLIMData::SLIM_ENERGY slim_energy, Eigen::VectorXi& b, Eigen::MatrixXd& bc, double soft_p);
+
+// Run iter_num iterations of SLIM
+// Outputs:
+//    V_o (in SLIMData): #V by dim list of mesh vertex positions
+IGL_INLINE void slim_solve(SLIMData& data, int iter_num);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slim.cpp"
+#endif
+
+
+#endif // SLIM_H

include/igl/triangle_triangle_adjacency.cpp

@@ -6,7 +6,6 @@
 // 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 "triangle_triangle_adjacency.h"
-#include "is_edge_manifold.h"
 #include "all_edges.h"
 #include "unique_simplices.h"
 #include "parallel_for.h"

include/igl/triangle_triangle_adjacency.h

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_TRIANGLE_TRIANGLE_ADJACENCY_H
 #define IGL_TRIANGLE_TRIANGLE_ADJACENCY_H
-#include <igl/igl_inline.h>
+#include "igl_inline.h"
 #include <Eigen/Core>
 #include <vector>
 

include/igl/volume.cpp

@@ -1,16 +1,16 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2014 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 "volume.h"
 #include "cross.h"
 #include <Eigen/Geometry>
 template <
-  typename DerivedV, 
-  typename DerivedT, 
+  typename DerivedV,
+  typename DerivedT,
   typename Derivedvol>
 IGL_INLINE void igl::volume(
   const Eigen::PlainObjectBase<DerivedV>& V,
@@ -68,7 +68,7 @@ IGL_INLINE typename VecA::Scalar igl::volume_single(
 
 
 template <
-  typename DerivedL, 
+  typename DerivedL,
   typename Derivedvol>
 IGL_INLINE void igl::volume(
   const Eigen::PlainObjectBase<DerivedL>& L,
@@ -92,10 +92,10 @@ IGL_INLINE void igl::volume(
     const ScalarS y = (V - w + u)*(w - u + V);
     const ScalarS Z = (v - W + u)*(W + u + v);
     const ScalarS z = (W - u + v)*(u - v + W);
-    const ScalarS a = sqrt(x*Y*Z); 
-    const ScalarS b = sqrt(y*Z*X); 
-    const ScalarS c = sqrt(z*X*Y); 
-    const ScalarS d = sqrt(x*y*z); 
+    const ScalarS a = sqrt(x*Y*Z);
+    const ScalarS b = sqrt(y*Z*X);
+    const ScalarS c = sqrt(z*X*Y);
+    const ScalarS d = sqrt(x*y*z);
     vol(t) = sqrt(
        (-a + b + c + d)*
        ( a - b + c + d)*
@@ -112,4 +112,6 @@ template Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar igl::volume_single<Eigen::
 template Eigen::Matrix<double, 3, 1, 0, 3, 1>::Scalar igl::volume_single<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&);
 template void igl::volume<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> >&);
 template void igl::volume<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::volume<class Eigen::Matrix<double,-1,3,0,-1,3>,class Eigen::Matrix<int,-1,3,0,-1,3>,class Eigen::Matrix<double,-1,1,0,-1,1> >(class Eigen::PlainObjectBase<class Eigen::Matrix<double,-1,3,0,-1,3> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<int,-1,3,0,-1,3> > const &,class Eigen::PlainObjectBase<class Eigen::Matrix<double,-1,1,0,-1,1> > &);
+
 #endif

include/igl/writeDMAT.cpp

@@ -18,7 +18,7 @@ IGL_INLINE bool igl::writeDMAT(
   const Mat & W,
   const bool ascii)
 {
-  FILE * fp = fopen(file_name.c_str(),"w");
+  FILE * fp = fopen(file_name.c_str(),"wb");
   if(fp == NULL)
   {
     fprintf(stderr,"IOError: writeDMAT() could not open %s...",file_name.c_str());

include/igl/writeOFF.cpp

@@ -43,4 +43,5 @@ template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matri
 template bool igl::writeOFF<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&);
 template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
 template bool igl::writeOFF<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
+template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&);
 #endif

index.html

@@ -59,6 +59,10 @@ and Windows with Visual Studio 2015 Community Edition.</p>
 <p>As of version 1.0, libigl includes an introductory
 <a href="http://libigl.github.io/libigl/tutorial/tutorial.html">tutorial</a> that covers many functionalities.</p>
 
+<h2 id="libiglexampleproject">libigl example project</h2>
+
+<p>We provide a <a href="https://github.com/libigl/libigl-example-project">blank project example</a> showing how to use libigl and cmake. Feel free and encouraged to copy or fork this project as a way of starting a new personal project using libigl.</p>
+
 <h2 id="installation">Installation</h2>
 
 <p>Libigl is a <strong>header-only</strong> library. You do <strong>not</strong> need to build anything to
@@ -213,8 +217,11 @@ Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:</p>
 
 <ul>
+<li><a href="http://www.activision.com">Activision</a></li>
 <li><a href="http://www.adobe.com/technology/">Adobe Research</a></li>
 <li><a href="http://www.ea.com">Electronic Arts, Inc</a></li>
+<li><a href="https://epicgames.com">Epic Games</a></li>
+<li><a href="https://research.google.com">Google Research</a></li>
 <li><a href="http://meshconsultants.ca/">Mesh</a>, consultants, Canada</li>
 <li><a href="http://graphics.pixar.com/research/">Pixar Research</a></li>
 <li><a href="http://esotericsoftware.com/">Spine by Esoteric Software</a> is an animation tool dedicated to 2D characters.</li>
@@ -225,6 +232,7 @@ few labs/companies/institutions using libigl:</p>
 <li>ETH Zurich, <a href="http://igl.ethz.ch/">Interactive Geometry Lab</a> and <a href="http://ait.inf.ethz.ch/">Advanced Technologies Lab</a>, Swizterland</li>
 <li>George Mason University, <a href="http://cs.gmu.edu/~ygingold/">CraGL</a>, USA</li>
 <li><a href="http://www.ust.hk/">Hong Kong University of Science and Technology</a>, Hong Kong</li>
+<li>[Inria](Université Grenoble Alpes), France</li>
 <li><a href="http://www.nii.ac.jp/en/">National Institute of Informatics</a>, Japan</li>
 <li>New York University, <a href="http://mrl.nyu.edu/">Media Research Lab</a>, USA</li>
 <li>NYUPoly, <a href="http://game.engineering.nyu.edu/">Game Innovation Lab</a>, USA</li>
@@ -232,13 +240,16 @@ few labs/companies/institutions using libigl:</p>
 <li><a href="http://www.tudelft.nl/en/">TU Delft</a>, Netherlands</li>
 <li><a href="https://www.tuwien.ac.at/en/tuwien_home/">TU Wien</a>, Austria</li>
 <li><a href="http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html">Telecom ParisTech</a>, Paris, France</li>
+<li><a href="http://www.staff.science.uu.nl/~vaxma001/">Utrecht University</a>, The Netherlands</li>
 <li><a href="http://mtm.ufsc.br/~leo/">Universidade Federal de Santa Catarina</a>, Brazil</li>
 <li><a href="http://vecg.cs.ucl.ac.uk/">University College London</a>, England</li>
 <li><a href="http://vis.berkeley.edu/">University of California Berkeley</a>, USA</li>
 <li><a href="http://www.cam.ac.uk/">University of Cambridge</a>, England</li>
 <li><a href="http://cg.cis.upenn.edu/">University of Pennsylvania</a>, USA</li>
 <li><a href="http://www.cs.utexas.edu/users/evouga/">University of Texas at Austin</a>, USA</li>
+<li><a href="http://dgp.toronto.edu">University of Toronto</a>, Canada</li>
 <li><a href="https://www.csc.uvic.ca/Research/graphics/">University of Victoria</a>, Canada</li>
+<li><a href="http://www.uwec.edu/computer-science/">University of Wisconsin-Eau Claire</a>, USA</li>
 <li><a href="http://www.usi.ch/en">Università della Svizzera Italiana</a>, Switzerland</li>
 <li><a href="http://www.univ-tlse3.fr/">Université Toulouse III Paul Sabatier</a>, France</li>
 <li><a href="http://www.math.zju.edu.cn/cagd/">Zhejiang University</a>, China</li>
@@ -264,7 +275,7 @@ page</a>.</p>
 <h2 id="copyright">Copyright</h2>
 
 <p>2016 Alec Jacobson, Daniele Panozzo, Christian Schüller, Olga Diamanti, Qingnan
-Zhou, Sebastian Koch, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
+Zhou, Sebastian Koch, Amir Vaxman, Nico Pietroni, Stefan Brugger, Kenshi Takayama, Wenzel Jakob, Nikolas De
 Giorgis, Luigi Rocca, Leonardo Sacht, Kevin Walliman, Olga Sorkine-Hornung, and others.</p>
 
 <p>Please see individual files for appropriate copyright notices.</p>

python/CMakeLists.txt

@@ -120,6 +120,11 @@ if (LIBIGL_WITH_COPYLEFT)
   list(APPEND SHARED_SOURCES "modules/copyleft/py_igl_copyleft.cpp")
 endif ()
 
+if (LIBIGL_WITH_BBW)
+  add_definitions(-DPY_BBW)
+  list(APPEND SHARED_SOURCES "modules/py_igl_bbw.cpp")
+endif ()
+
 if (LIBIGL_WITH_PNG)
   add_definitions(-DPY_PNG)
   list(APPEND SHARED_SOURCES "modules/py_igl_png.cpp")

python/README.md

@@ -128,9 +128,11 @@ We provide a few examples in the folder python/matlab.
 ## Documentation
 
 The python functions have exactly the same prototypes as their C++ counterpart.
-To get help for a certain function, please check the documentation in the
-corresponding .h file in libigl/include. We will add a proper docstring
-documentation in the future.
+Docstrings for all available python functions are extracted from the C++ header files and compiled into the python module. To get help for a certain function, you can run `help(pyigl.<function_name>)` in the python console.
+
+In the scripts folder there is the script `generate_docstrings.py` that automatically generates python docstrings for a new function. You can run it with `generate_docstrings.py <path_to_cpp_header_files> <path_to_python_files>`. 
+The script depends on additional libraries (joblib, mako, clang), make sure to install them (e.g. through pip. python-clang is included in external/nanogui/ext/pybind11/tools/clang).
+
 
 ## Known Issues
 

python/modules/py_igl_bbw.cpp

@@ -0,0 +1,18 @@
+//#include <Eigen/Geometry>
+//#include <Eigen/Dense>
+//#include <Eigen/Sparse>
+
+
+#include "../python_shared.h"
+
+#include <igl/bbw/bbw.h>
+
+
+void python_export_igl_bbw(py::module &me) {
+
+  py::module m = me.def_submodule(
+    "bbw", "Wrappers for libigl functions that use bbw");
+
+  #include "../py_igl/bbw/py_bbw.cpp"
+
+}

python/modules/py_igl_viewer.cpp

@@ -188,9 +188,15 @@ py::class_<igl::viewer::ViewerCore> viewercore_class(me, "ViewerCore");
     })
 
     .def_readwrite("lighting_factor",&igl::viewer::ViewerCore::lighting_factor)
-
     .def_readwrite("model_zoom",&igl::viewer::ViewerCore::model_zoom)
 
+    .def_property("trackball_angle",
+    [](const igl::viewer::ViewerCore& core) {return Eigen::Quaterniond(core.trackball_angle.cast<double>());},
+    [](igl::viewer::ViewerCore& core, const Eigen::Quaterniond& q)
+    {
+      core.trackball_angle = Eigen::Quaternionf(q.cast<float>());
+    })
+
     .def_property("model_translation",
     [](const igl::viewer::ViewerCore& core) {return Eigen::MatrixXd(core.model_translation.cast<double>());},
     [](igl::viewer::ViewerCore& core, const Eigen::MatrixXd& v)
@@ -316,13 +322,13 @@ py::class_<igl::viewer::ViewerCore> viewercore_class(me, "ViewerCore");
 // UI Enumerations
     py::class_<igl::viewer::Viewer> viewer_class(me, "Viewer");
 
-    #ifdef IGL_VIEWER_WITH_NANOGUI
-    py::object fh = (py::object) py::module::import("nanogui").attr("FormHelper");
-    py::class_<nanogui::FormHelper> form_helper_class(me, "FormHelper", fh);
+//    #ifdef IGL_VIEWER_WITH_NANOGUI
+//    py::object fh = (py::object) py::module::import("nanogui").attr("FormHelper");
+//    py::class_<nanogui::FormHelper> form_helper_class(me, "FormHelper", fh);
 
-    py::object screen = (py::object) py::module::import("nanogui").attr("Screen");
-    py::class_<nanogui::Screen> screen_class(me, "Screen", screen);
-    #endif
+//    py::object screen = (py::object) py::module::import("nanogui").attr("Screen");
+//    py::class_<nanogui::Screen> screen_class(me, "Screen", screen);
+//    #endif
 
     py::enum_<igl::viewer::Viewer::MouseButton>(viewer_class, "MouseButton")
         .value("Left", igl::viewer::Viewer::MouseButton::Left)

python/modules/py_typedefs.cpp

@@ -0,0 +1,13 @@
+py::class_<RotationList>(m, "RotationList")
+    .def(py::init<>())
+    .def(py::init<size_t>())
+    .def("pop_back", &RotationList::pop_back)
+    /* There are multiple versions of push_back(), etc. Select the right ones. */
+    .def("append", (void (RotationList::*)(const Eigen::Quaterniond &)) &RotationList::push_back)
+    .def("back", (Eigen::Quaterniond &(RotationList::*)()) &RotationList::back)
+    .def("__len__", [](const RotationList &v) { return v.size(); })
+    .def("__getitem__", [](const RotationList &v, int b) { return v.at(b); })
+    .def("__setitem__", [](RotationList &v, int b, Eigen::Quaterniond &c) { return v.at(b) = c; })
+    .def("__iter__", [](RotationList &v) {
+       return py::make_iterator(v.begin(), v.end());
+}, py::keep_alive<0, 1>());

python/modules/py_typedefs.h

@@ -0,0 +1,5 @@
+typedef std::vector<Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> > RotationList;
+PYBIND11_MAKE_OPAQUE(RotationList);
+
+//typedef std::vector<Eigen::Vector3d> TranslationList;
+//PYBIND11_MAKE_OPAQUE(TranslationList);

python/modules/py_vector.cpp

@@ -68,7 +68,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         })
         .def("__init__", [](Type &m, py::buffer b) {
             py::buffer_info info = b.request();
-            if (info.format != py::format_descriptor<Scalar>::value())
+            if (info.format != py::format_descriptor<Scalar>::value)
                 throw std::runtime_error("Incompatible buffer format!");
             if (info.ndim == 1) {
                 new (&m) Type(info.shape[0], 1);
@@ -319,7 +319,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
                 m.data(),                /* Pointer to buffer */
                 sizeof(Scalar),          /* Size of one scalar */
                 /* Python struct-style format descriptor */
-                py::format_descriptor<Scalar>::value(),
+                py::format_descriptor<Scalar>::value,
                 2,                       /* Number of dimensions */
                 { (size_t) m.rows(),     /* Buffer dimensions */
                   (size_t) m.cols() },
@@ -692,10 +692,12 @@ void python_export_vector(py::module &m) {
     .def("solve",[](const Eigen::SimplicialLLT<Eigen::SparseMatrix<double > >& s, const Eigen::MatrixXd& rhs) { return Eigen::MatrixXd(s.solve(rhs)); })
     ;
 
+    // Bindings for Affine3d
     py::class_<Eigen::Affine3d > affine3d(me, "Affine3d");
 
     affine3d
     .def(py::init<>())
+    .def_static("Identity", []() { return Eigen::Affine3d::Identity(); })
     .def("setIdentity",[](Eigen::Affine3d& a){
         return a.setIdentity();
     })
@@ -703,6 +705,9 @@ void python_export_vector(py::module &m) {
         assert_is_Vector3("axis", axis);
         return a.rotate(Eigen::AngleAxisd(angle, Eigen::Vector3d(axis)));
     })
+    .def("rotate",[](Eigen::Affine3d& a, Eigen::Quaterniond quat) {
+        return a.rotate(quat);
+    })
     .def("translate",[](Eigen::Affine3d& a, Eigen::MatrixXd offset) {
         assert_is_Vector3("offset", offset);
         return a.translate(Eigen::Vector3d(offset));
@@ -711,13 +716,51 @@ void python_export_vector(py::module &m) {
         return Eigen::MatrixXd(a.matrix());
     })
     ;
-    /* Bindings for Quaterniond*/
-    //py::class_<Eigen::Quaterniond > quaterniond(me, "Quaterniond");
-    //
-    // quaterniond
-    // .def(py::init<>())
-    // .def(py::init<double, double, double, double>())
-    // ;
+    // Bindings for Quaterniond
+    py::class_<Eigen::Quaterniond > quaterniond(me, "Quaterniond");
+    
+    quaterniond
+    .def(py::init<>())
+    .def(py::init<double, double, double, double>())
+    .def("__init__", [](Eigen::Quaterniond &q, double angle, Eigen::MatrixXd axis) {
+        assert_is_Vector3("axis", axis);
+        new (&q) Eigen::Quaterniond(Eigen::AngleAxisd(angle, Eigen::Vector3d(axis)));
+    })
+    .def_static("Identity", []() { return Eigen::Quaterniond::Identity(); })
+    .def("__repr__", [](const Eigen::Quaterniond &v) {
+        std::ostringstream oss;
+        oss << "(" << v.w() << ", " << v.x() << ", " << v.y() << ", " << v.z() << ")";
+        return oss.str();
+    })
+    .def("conjugate",[](Eigen::Quaterniond& q) {
+        return q.conjugate();
+    })
+    .def("normalize",[](Eigen::Quaterniond& q) {
+        return q.normalize();
+    })
+    .def("slerp",[](Eigen::Quaterniond& q, double & t, Eigen::Quaterniond other) {
+        return q.slerp(t, other);
+    })
+//    .def_cast(-py::self)
+//    .def_cast(py::self + py::self)
+//    .def_cast(py::self - py::self)
+    .def_cast(py::self * py::self)
+    // .def_cast(py::self - Scalar())
+    // .def_cast(py::self * Scalar())
+    // .def_cast(py::self / Scalar())
+
+//    .def("__mul__", []
+//    (const Type &a, const Scalar& b)
+//    {
+//      return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(a * b);
+//    })
+//    .def("__rmul__", [](const Type& a, const Scalar& b)
+//    {
+//      return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(b * a);
+//    })
+    ;
+
+    
 
 
 

python/py_doc.cpp

@@ -90,6 +90,45 @@ const char *__doc_igl_barycentric_to_global = R"igl_Qu8mg5v7(// Converts barycen
   // Output:
   // #X: #Xx3 3D coordinates of all points in bc
   //   )igl_Qu8mg5v7";
+const char *__doc_igl_bbw_bbw = R"igl_Qu8mg5v7(// Compute Bounded Biharmonic Weights on a given domain (V,Ele) with a given
+    // set of boundary conditions
+    //
+    // Templates
+    //   DerivedV  derived type of eigen matrix for V (e.g. MatrixXd)
+    //   DerivedF  derived type of eigen matrix for F (e.g. MatrixXi)
+    //   Derivedb  derived type of eigen matrix for b (e.g. VectorXi)
+    //   Derivedbc  derived type of eigen matrix for bc (e.g. MatrixXd)
+    //   DerivedW  derived type of eigen matrix for W (e.g. MatrixXd)
+    // Inputs:
+    //   V  #V by dim vertex positions
+    //   Ele  #Elements by simplex-size list of element indices
+    //   b  #b boundary indices into V
+    //   bc #b by #W list of boundary values
+    //   data  object containing options, intial guess --> solution and results
+    // Outputs:
+    //   W  #V by #W list of *unnormalized* weights to normalize use
+    //    igl::normalize_row_sums(W,W);
+    // Returns true on success, false on failure)igl_Qu8mg5v7";
+const char *__doc_igl_boundary_conditions = R"igl_Qu8mg5v7(// Compute boundary conditions for automatic weights computation. This
+  // function expects that the given mesh (V,Ele) has sufficient samples
+  // (vertices) exactly at point handle locations and exactly along bone and
+  // cage edges.
+  //
+  // Inputs:
+  //   V  #V by dim list of domain vertices
+  //   Ele  #Ele by simplex-size list of simplex indices
+  //   C  #C by dim list of handle positions
+  //   P  #P by 1 list of point handle indices into C
+  //   BE  #BE by 2 list of bone edge indices into C
+  //   CE  #CE by 2 list of cage edge indices into *P*
+  // Outputs:
+  //   b  #b list of boundary indices (indices into V of vertices which have
+  //     known, fixed values)
+  //   bc #b by #weights list of known/fixed values for boundary vertices
+  //     (notice the #b != #weights in general because #b will include all the
+  //     intermediary samples along each bone, etc.. The ordering of the
+  //     weights corresponds to [P;BE]
+  // Returns true if boundary conditions make sense)igl_Qu8mg5v7";
 const char *__doc_igl_boundary_facets = R"igl_Qu8mg5v7(// BOUNDARY_FACETS Determine boundary faces (edges) of tetrahedra (triangles)
   // stored in T (analogous to qptoolbox's `outline` and `boundary_faces`).
   //
@@ -146,6 +185,13 @@ const char *__doc_igl_colon = R"igl_Qu8mg5v7(// Colon operator like matlab's col
   //     than hi, vice versa if hi<low
   // Output:
   //   I  list of values from low to hi with step size step)igl_Qu8mg5v7";
+const char *__doc_igl_column_to_quats = R"igl_Qu8mg5v7(// "Columnize" a list of quaternions (q1x,q1y,q1z,q1w,q2x,q2y,q2z,q2w,...)
+  //
+  // Inputs:
+  //   Q  n*4-long list of coefficients
+  // Outputs:
+  //   vQ  n-long list of quaternions
+  // Returns false if n%4!=0)igl_Qu8mg5v7";
 const char *__doc_igl_comb_cross_field = R"igl_Qu8mg5v7(// Inputs:
   //   V          #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F          #F by 4 eigen Matrix of face (quad) indices
@@ -229,7 +275,12 @@ const char *__doc_igl_copyleft_cgal_remesh_self_intersections = R"igl_Qu8mg5v7(/
       //     // remove any vertices now unreferenced after duplicate mapping.
       //     igl::remove_unreferenced(VV,FF,SV,SF,UIM);
       //     // Now (SV,SF) is ready to extract outer hull
+<<<<<<< HEAD
       //     igl::copyleft::cgal::outer_hull(SV,SF,G,J,flip);)igl_Qu8mg5v7";
+=======
+      //     igl::copyleft::cgal::outer_hull(SV,SF,G,J,flip);
+      //)igl_Qu8mg5v7";
+>>>>>>> upstream/master
 const char *__doc_igl_copyleft_comiso_miq = R"igl_Qu8mg5v7(// Inputs:
     //   V              #V by 3 list of mesh vertex 3D positions
     //   F              #F by 3 list of faces indices in V
@@ -383,6 +434,32 @@ const char *__doc_igl_cut_mesh_from_singularities = R"igl_Qu8mg5v7(// Given a me
   //   seams  #F by 3 list of per corner booleans that denotes if an edge is a
   //     seam or not
   //)igl_Qu8mg5v7";
+const char *__doc_igl_deform_skeleton = R"igl_Qu8mg5v7(// Deform a skeleton.
+  //
+  // Inputs:
+  //   C  #C by 3 list of joint positions
+  //   BE  #BE by 2 list of bone edge indices
+  //   vA  #BE list of bone transformations
+  // Outputs
+  //   CT  #BE*2 by 3 list of deformed joint positions
+  //   BET  #BE by 2 list of bone edge indices (maintains order)
+  //)igl_Qu8mg5v7";
+const char *__doc_igl_directed_edge_orientations = R"igl_Qu8mg5v7(// Determine rotations that take each edge from the x-axis to its given rest
+  // orientation.
+  //
+  // Inputs:
+  //   C  #C by 3 list of edge vertex positions
+  //   E  #E by 2 list of directed edges
+  // Outputs:
+  //   Q  #E list of quaternions 
+  //)igl_Qu8mg5v7";
+const char *__doc_igl_directed_edge_parents = R"igl_Qu8mg5v7(// Recover "parents" (preceeding edges) in a tree given just directed edges.
+  //
+  // Inputs:
+  //   E  #E by 2 list of directed edges
+  // Outputs:
+  //   P  #E list of parent indices into E (-1) means root
+  //)igl_Qu8mg5v7";
 const char *__doc_igl_doublearea = R"igl_Qu8mg5v7(// DOUBLEAREA computes twice the area for each input triangle[quad]
   //
   // Templates:
@@ -411,6 +488,15 @@ const char *__doc_igl_doublearea_quad = R"igl_Qu8mg5v7(// DOUBLEAREA_QUAD comput
   // Outputs:
   //   dblA  #F list of quadrilateral double areas
   //)igl_Qu8mg5v7";
+const char *__doc_igl_dqs = R"igl_Qu8mg5v7(// Dual quaternion skinning
+  //
+  // Inputs:
+  //   V  #V by 3 list of rest positions
+  //   W  #W by #C list of weights
+  //   vQ  #C list of rotation quaternions 
+  //   vT  #C list of translation vectors
+  // Outputs:
+  //   U  #V by 3 list of new positions)igl_Qu8mg5v7";
 const char *__doc_igl_edge_lengths = R"igl_Qu8mg5v7(// Constructs a list of lengths of edges opposite each index in a face
   // (triangle/tet) list
   //
@@ -426,7 +512,7 @@ const char *__doc_igl_edge_lengths = R"igl_Qu8mg5v7(// Constructs a list of leng
   //    or
   //   T  #T by 4 list of mesh elements (must be tets)
   // Outputs:
-  //   L  #F by {1|3|6} list of edge lengths
+  //   L  #F by {1|3|6} list of edge lengths 
   //     for edges, column of lengths
   //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
   //     for tets, columns correspond to edges
@@ -438,7 +524,8 @@ const char *__doc_igl_edge_topology = R"igl_Qu8mg5v7(// Initialize Edges and the
   // EV  : #Ex2, Stores the edge description as pair of indices to vertices
   // FE : #Fx3, Stores the Triangle-Edge relation
   // EF : #Ex2: Stores the Edge-Triangle relation
-  //)igl_Qu8mg5v7";
+  //
+  // TODO: This seems to be a duplicate of edge_flaps.h)igl_Qu8mg5v7";
 const char *__doc_igl_eigs = R"igl_Qu8mg5v7(See eigs for the documentation.)igl_Qu8mg5v7";
 const char *__doc_igl_embree_ambient_occlusion = R"igl_Qu8mg5v7(// Compute ambient occlusion per given point
     //
@@ -522,6 +609,18 @@ const char *__doc_igl_floor = R"igl_Qu8mg5v7(// Floor a given matrix to nearest
   //   X  m by n matrix of scalars
   // Outputs:
   //   Y  m by n matrix of floored integers)igl_Qu8mg5v7";
+const char *__doc_igl_forward_kinematics = R"igl_Qu8mg5v7(// Given a skeleton and a set of relative bone rotations compute absolute
+  // rigid transformations for each bone.
+  //
+  // Inputs:
+  //   C  #C by dim list of joint positions
+  //   BE  #BE by 2 list of bone edge indices
+  //   P  #BE list of parent indices into BE
+  //   dQ  #BE list of relative rotations
+  //   dT  #BE list of relative translations
+  // Outputs:
+  //   vQ  #BE list of absolute rotations
+  //   vT  #BE list of absolute translations)igl_Qu8mg5v7";
 const char *__doc_igl_gaussian_curvature = R"igl_Qu8mg5v7(// Compute discrete local integral gaussian curvature (angle deficit, without
   // averaging by local area).
   //
@@ -547,12 +646,12 @@ const char *__doc_igl_get_seconds = R"igl_Qu8mg5v7(// Return the current time in
   //    ... // part 2
   //    cout<<"part 2: "<<tictoc()<<endl;
   //    ... // etc)igl_Qu8mg5v7";
-const char *__doc_igl_grad = R"igl_Qu8mg5v7(// Gradient of a scalar function defined on piecewise linear elements (mesh)
-  // is constant on each triangle i,j,k:
-  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
-  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
-  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
-  // 90 degrees
+const char *__doc_igl_grad = R"igl_Qu8mg5v7(// Gradient of a scalar function defined on piecewise linear elements (mesh)
+  // is constant on each triangle i,j,k:
+  // grad(Xijk) = (Xj-Xi) * (Vi - Vk)^R90 / 2A + (Xk-Xi) * (Vj - Vi)^R90 / 2A
+  // where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
+  // i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
+  // 90 degrees
   //)igl_Qu8mg5v7";
 const char *__doc_igl_harmonic = R"igl_Qu8mg5v7(// Compute k-harmonic weight functions "coordinates".
   //
@@ -620,6 +719,39 @@ const char *__doc_igl_jet = R"igl_Qu8mg5v7(// JET like MATLAB's jet
   //   r  red value
   //   g  green value
   //   b  blue value)igl_Qu8mg5v7";
+const char *__doc_igl_lbs_matrix = R"igl_Qu8mg5v7(// LBS_MATRIX Linear blend skinning can be expressed by V' = M * T where V' is
+  // a #V by dim matrix of deformed vertex positions (one vertex per row), M is a
+  // #V by (dim+1)*#T (composed of weights and rest positions) and T is a
+  // #T*(dim+1) by dim matrix of #T stacked transposed transformation matrices.
+  // See equations (1) and (2) in "Fast Automatic Skinning Transformations"
+  // [Jacobson et al 2012]
+  //
+  // Inputs:
+  //   V  #V by dim list of rest positions
+  //   W  #V+ by #T  list of weights
+  // Outputs:
+  //   M  #V by #T*(dim+1)
+  //
+  // In MATLAB:
+  //   kron(ones(1,size(W,2)),[V ones(size(V,1),1)]).*kron(W,ones(1,size(V,2)+1)))igl_Qu8mg5v7";
+const char *__doc_igl_lbs_matrix_column = R"igl_Qu8mg5v7(// LBS_MATRIX  construct a matrix that when multiplied against a column of
+  // affine transformation entries computes new coordinates of the vertices
+  //
+  // I'm not sure it makes since that the result is stored as a sparse matrix.
+  // The number of non-zeros per row *is* dependent on the number of mesh
+  // vertices and handles.
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex rest positions
+  //   W  #V by #handles list of correspondence weights
+  // Output:
+  //   M  #V * dim by #handles * dim * (dim+1) matrix such that
+  //     new_V(:) = LBS(V,W,A) = reshape(M * A,size(V)), where A is a column
+  //     vectors formed by the entries in each handle's dim by dim+1 
+  //     transformation matrix. Specifcally, A =
+  //       reshape(permute(Astack,[3 1 2]),n*dim*(dim+1),1)
+  //     or A = [Lxx;Lyx;Lxy;Lyy;tx;ty], and likewise for other dim
+  //     if Astack(:,:,i) is the dim by (dim+1) transformation at handle i)igl_Qu8mg5v7";
 const char *__doc_igl_local_basis = R"igl_Qu8mg5v7(// Compute a local orthogonal reference system for each triangle in the given mesh
   // Templates:
   //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
@@ -731,6 +863,22 @@ const char *__doc_igl_n_polyvector = R"igl_Qu8mg5v7(// Inputs:
   // Output:
   //                  3 by 3 rotation matrix that takes v0 to v1
   //)igl_Qu8mg5v7";
+const char *__doc_igl_normalize_row_lengths = R"igl_Qu8mg5v7(// Obsolete: just use A.rowwise().normalize() or B=A.rowwise().normalized();
+  //
+  // Normalize the rows in A so that their lengths are each 1 and place the new
+  // entries in B
+  // Inputs:
+  //   A  #rows by k input matrix
+  // Outputs:
+  //   B  #rows by k input matrix, can be the same as A)igl_Qu8mg5v7";
+const char *__doc_igl_normalize_row_sums = R"igl_Qu8mg5v7(// Normalize the rows in A so that their sums are each 1 and place the new
+  // entries in B
+  // Inputs:
+  //   A  #rows by k input matrix
+  // Outputs:
+  //   B  #rows by k input matrix, can be the same as A
+  //
+  // Note: This is just calling an Eigen one-liner.)igl_Qu8mg5v7";
 const char *__doc_igl_parula = R"igl_Qu8mg5v7(// PARULA like MATLAB's parula
   //
   // Inputs:
@@ -923,6 +1071,23 @@ const char *__doc_igl_readOFF = R"igl_Qu8mg5v7(// Read a mesh from an ascii obj
   //   N  double matrix of corner normals #N by 3
   //   FN  #F list of face indices into vertex normals
   // Returns true on success, false on errors)igl_Qu8mg5v7";
+const char *__doc_igl_readTGF = R"igl_Qu8mg5v7(// READTGF
+  //
+  // [V,E,P,BE,CE,PE] = readTGF(filename)
+  //
+  // Read a graph from a .tgf file
+  //
+  // Input:
+  //  filename  .tgf file name
+  // Ouput:
+  //  V  # vertices by 3 list of vertex positions
+  //  E  # edges by 2 list of edge indices
+  //  P  # point-handles list of point handle indices
+  //  BE # bone-edges by 2 list of bone-edge indices
+  //  CE # cage-edges by 2 list of cage-edge indices
+  //  PE # pseudo-edges by 2 list of pseudo-edge indices
+  // 
+  // Assumes that graph vertices are 3 dimensional)igl_Qu8mg5v7";
 const char *__doc_igl_read_triangle_mesh = R"igl_Qu8mg5v7(// read mesh from an ascii file with automatic detection of file format.
   // supported: obj, off, stl, wrl, ply, mesh)
   // 
@@ -936,7 +1101,7 @@ const char *__doc_igl_read_triangle_mesh = R"igl_Qu8mg5v7(// read mesh from an a
   //   V  eigen double matrix #V by 3
   //   F  eigen int matrix #F by 3
   // Returns true iff success)igl_Qu8mg5v7";
-const char *__doc_igl_remove_duplicate_vertices = R"igl_Qu8mg5v7(  // REMOVE_DUPLICATE_VERTICES Remove duplicate vertices upto a uniqueness
+const char *__doc_igl_remove_duplicate_vertices = R"igl_Qu8mg5v7(// REMOVE_DUPLICATE_VERTICES Remove duplicate vertices upto a uniqueness
   // tolerance (epsilon)
   //
   // Inputs:
@@ -945,7 +1110,7 @@ const char *__doc_igl_remove_duplicate_vertices = R"igl_Qu8mg5v7(  // REMOVE_DUP
   //     this as a tolerance on L1 distance
   // Outputs:
   //   SV  #SV by dim new list of vertex positions
-  //   SVI #V by 1 list of indices so SV = V(SVI,:)
+  //   SVI #V by 1 list of indices so SV = V(SVI,:) 
   //   SVJ #SV by 1 list of indices so V = SV(SVJ,:)
   //
   // Example:
@@ -953,7 +1118,6 @@ const char *__doc_igl_remove_duplicate_vertices = R"igl_Qu8mg5v7(  // REMOVE_DUP
   //   [SV,SVI,SVJ] = remove_duplicate_vertices(V,1e-7);
   //   % remap faces
   //   SF = SVJ(F);
-  //
   //)igl_Qu8mg5v7";
 const char *__doc_igl_rotate_vectors = R"igl_Qu8mg5v7(// Rotate the vectors V by A radiants on the tangent plane spanned by B1 and
   // B2
@@ -1072,7 +1236,7 @@ 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'
+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
@@ -1084,12 +1248,12 @@ const char *__doc_igl_streamlines_init = R"igl_Qu8mg5v7(    // Given a mesh and
     //   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
+    //   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";
+    //   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).
   //
@@ -1106,6 +1270,9 @@ const char *__doc_igl_triangle_triangle_adjacency = R"igl_Qu8mg5v7(// Constructs
   // NOTE: the first edge of a triangle is [0,1] the second [1,2] and the third [2,3].
   //       this convention is DIFFERENT from cotmatrix_entries.h
   // Known bug: this should not need to take V as input.)igl_Qu8mg5v7";
+const char *__doc_igl_triangle_triangle_adjacency_preprocess = R"igl_Qu8mg5v7(// Preprocessing)igl_Qu8mg5v7";
+const char *__doc_igl_triangle_triangle_adjacency_extractTT = R"igl_Qu8mg5v7(// Extract the face adjacencies)igl_Qu8mg5v7";
+const char *__doc_igl_triangle_triangle_adjacency_extractTTi = R"igl_Qu8mg5v7(// Extract the face adjacencies indices (needed for fast traversal))igl_Qu8mg5v7";
 const char *__doc_igl_triangle_triangulate = R"igl_Qu8mg5v7(// Triangulate the interior of a polygon using the triangle library.
     //
     // Inputs:

python/py_doc.h

@@ -5,11 +5,14 @@ extern const char *__doc_igl_avg_edge_length;
 extern const char *__doc_igl_barycenter;
 extern const char *__doc_igl_barycentric_coordinates;
 extern const char *__doc_igl_barycentric_to_global;
+extern const char *__doc_igl_bbw_bbw;
+extern const char *__doc_igl_boundary_conditions;
 extern const char *__doc_igl_boundary_facets;
 extern const char *__doc_igl_boundary_loop;
 extern const char *__doc_igl_cat;
 extern const char *__doc_igl_collapse_edge;
 extern const char *__doc_igl_colon;
+extern const char *__doc_igl_column_to_quats;
 extern const char *__doc_igl_comb_cross_field;
 extern const char *__doc_igl_comb_frame_field;
 extern const char *__doc_igl_compute_frame_field_bisectors;
@@ -24,9 +27,13 @@ extern const char *__doc_igl_cotmatrix;
 extern const char *__doc_igl_covariance_scatter_matrix;
 extern const char *__doc_igl_cross_field_missmatch;
 extern const char *__doc_igl_cut_mesh_from_singularities;
+extern const char *__doc_igl_deform_skeleton;
+extern const char *__doc_igl_directed_edge_orientations;
+extern const char *__doc_igl_directed_edge_parents;
 extern const char *__doc_igl_doublearea;
 extern const char *__doc_igl_doublearea_single;
 extern const char *__doc_igl_doublearea_quad;
+extern const char *__doc_igl_dqs;
 extern const char *__doc_igl_edge_lengths;
 extern const char *__doc_igl_edge_topology;
 extern const char *__doc_igl_eigs;
@@ -38,6 +45,7 @@ extern const char *__doc_igl_fit_rotations;
 extern const char *__doc_igl_fit_rotations_planar;
 extern const char *__doc_igl_fit_rotations_SSE;
 extern const char *__doc_igl_floor;
+extern const char *__doc_igl_forward_kinematics;
 extern const char *__doc_igl_gaussian_curvature;
 extern const char *__doc_igl_get_seconds;
 extern const char *__doc_igl_grad;
@@ -47,6 +55,8 @@ extern const char *__doc_igl_internal_angles;
 extern const char *__doc_igl_invert_diag;
 extern const char *__doc_igl_is_irregular_vertex;
 extern const char *__doc_igl_jet;
+extern const char *__doc_igl_lbs_matrix;
+extern const char *__doc_igl_lbs_matrix_column;
 extern const char *__doc_igl_local_basis;
 extern const char *__doc_igl_lscm;
 extern const char *__doc_igl_map_vertices_to_circle;
@@ -55,6 +65,8 @@ extern const char *__doc_igl_min_quad_with_fixed_precompute;
 extern const char *__doc_igl_min_quad_with_fixed_solve;
 extern const char *__doc_igl_min_quad_with_fixed;
 extern const char *__doc_igl_n_polyvector;
+extern const char *__doc_igl_normalize_row_lengths;
+extern const char *__doc_igl_normalize_row_sums;
 extern const char *__doc_igl_parula;
 extern const char *__doc_igl_per_corner_normals;
 extern const char *__doc_igl_per_edge_normals;
@@ -73,6 +85,7 @@ extern const char *__doc_igl_readDMAT;
 extern const char *__doc_igl_readMESH;
 extern const char *__doc_igl_readOBJ;
 extern const char *__doc_igl_readOFF;
+extern const char *__doc_igl_readTGF;
 extern const char *__doc_igl_read_triangle_mesh;
 extern const char *__doc_igl_remove_duplicate_vertices;
 extern const char *__doc_igl_rotate_vectors;
@@ -88,6 +101,9 @@ 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_triangle_adjacency_preprocess;
+extern const char *__doc_igl_triangle_triangle_adjacency_extractTT;
+extern const char *__doc_igl_triangle_triangle_adjacency_extractTTi;
 extern const char *__doc_igl_triangle_triangulate;
 extern const char *__doc_igl_unique;
 extern const char *__doc_igl_unique_rows;

python/py_igl.cpp

@@ -1,6 +1,7 @@
 #include <Eigen/Dense>
 
 #include "python_shared.h"
+#include "modules/py_typedefs.h"
 
 #include <igl/AABB.h>
 #include <igl/ARAPEnergyType.h>
@@ -11,12 +12,17 @@
 #include <igl/avg_edge_length.h>
 #include <igl/barycenter.h>
 #include <igl/barycentric_coordinates.h>
+<<<<<<< HEAD
 #include <igl/barycentric_to_global.h>
+=======
+#include <igl/boundary_conditions.h>
+>>>>>>> upstream/master
 #include <igl/boundary_facets.h>
 #include <igl/boundary_loop.h>
 #include <igl/cat.h>
 #include <igl/collapse_edge.h>
 #include <igl/colon.h>
+#include <igl/column_to_quats.h>
 #include <igl/comb_cross_field.h>
 #include <igl/comb_frame_field.h>
 #include <igl/compute_frame_field_bisectors.h>
@@ -24,13 +30,18 @@
 #include <igl/covariance_scatter_matrix.h>
 #include <igl/cross_field_missmatch.h>
 #include <igl/cut_mesh_from_singularities.h>
+#include <igl/deform_skeleton.h>
+#include <igl/directed_edge_orientations.h>
+#include <igl/directed_edge_parents.h>
 #include <igl/doublearea.h>
+#include <igl/dqs.h>
 #include <igl/edge_lengths.h>
 #include <igl/edge_topology.h>
 #include <igl/eigs.h>
 #include <igl/find_cross_field_singularities.h>
 #include <igl/fit_rotations.h>
 #include <igl/floor.h>
+#include <igl/forward_kinematics.h>
 #include <igl/gaussian_curvature.h>
 #include <igl/get_seconds.h>
 #include <igl/grad.h>
@@ -40,12 +51,15 @@
 #include <igl/invert_diag.h>
 #include <igl/is_irregular_vertex.h>
 #include <igl/jet.h>
+#include <igl/lbs_matrix.h>
 #include <igl/local_basis.h>
 #include <igl/lscm.h>
 #include <igl/map_vertices_to_circle.h>
 #include <igl/massmatrix.h>
 #include <igl/min_quad_with_fixed.h>
 #include <igl/n_polyvector.h>
+#include <igl/normalize_row_lengths.h>
+#include <igl/normalize_row_sums.h>
 #include <igl/parula.h>
 #include <igl/per_corner_normals.h>
 #include <igl/per_edge_normals.h>
@@ -61,6 +75,7 @@
 #include <igl/readMESH.h>
 #include <igl/readOBJ.h>
 #include <igl/readOFF.h>
+#include <igl/readTGF.h>
 #include <igl/read_triangle_mesh.h>
 #include <igl/remove_duplicate_vertices.h>
 #include <igl/rotate_vectors.h>
@@ -83,6 +98,8 @@
 
 void python_export_igl(py::module &m)
 {
+#include "modules/py_typedefs.cpp"
+
 #include "py_igl/py_AABB.cpp"
 #include "py_igl/py_ARAPEnergyType.cpp"
 #include "py_igl/py_MeshBooleanType.cpp"
@@ -93,11 +110,13 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_barycenter.cpp"
 #include "py_igl/py_barycentric_coordinates.cpp"
 #include "py_igl/py_barycentric_to_global.cpp"
+#include "py_igl/py_boundary_conditions.cpp"
 #include "py_igl/py_boundary_facets.cpp"
 #include "py_igl/py_boundary_loop.cpp"
 #include "py_igl/py_cat.cpp"
 #include "py_igl/py_collapse_edge.cpp"
 #include "py_igl/py_colon.cpp"
+#include "py_igl/py_column_to_quats.cpp"
 #include "py_igl/py_comb_cross_field.cpp"
 #include "py_igl/py_comb_frame_field.cpp"
 #include "py_igl/py_compute_frame_field_bisectors.cpp"
@@ -105,13 +124,18 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_covariance_scatter_matrix.cpp"
 #include "py_igl/py_cross_field_missmatch.cpp"
 #include "py_igl/py_cut_mesh_from_singularities.cpp"
+#include "py_igl/py_deform_skeleton.cpp"
+#include "py_igl/py_directed_edge_orientations.cpp"
+#include "py_igl/py_directed_edge_parents.cpp"
 #include "py_igl/py_doublearea.cpp"
+#include "py_igl/py_dqs.cpp"
 #include "py_igl/py_edge_lengths.cpp"
 #include "py_igl/py_edge_topology.cpp"
 #include "py_igl/py_eigs.cpp"
 #include "py_igl/py_find_cross_field_singularities.cpp"
 #include "py_igl/py_fit_rotations.cpp"
 #include "py_igl/py_floor.cpp"
+#include "py_igl/py_forward_kinematics.cpp"
 #include "py_igl/py_gaussian_curvature.cpp"
 #include "py_igl/py_get_seconds.cpp"
 #include "py_igl/py_grad.cpp"
@@ -121,12 +145,15 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_invert_diag.cpp"
 #include "py_igl/py_is_irregular_vertex.cpp"
 #include "py_igl/py_jet.cpp"
+#include "py_igl/py_lbs_matrix.cpp"
 #include "py_igl/py_local_basis.cpp"
 #include "py_igl/py_lscm.cpp"
 #include "py_igl/py_map_vertices_to_circle.cpp"
 #include "py_igl/py_massmatrix.cpp"
 #include "py_igl/py_min_quad_with_fixed.cpp"
 #include "py_igl/py_n_polyvector.cpp"
+#include "py_igl/py_normalize_row_lengths.cpp"
+#include "py_igl/py_normalize_row_sums.cpp"
 #include "py_igl/py_parula.cpp"
 #include "py_igl/py_per_corner_normals.cpp"
 #include "py_igl/py_per_edge_normals.cpp"
@@ -142,6 +169,7 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_readMESH.cpp"
 #include "py_igl/py_readOBJ.cpp"
 #include "py_igl/py_readOFF.cpp"
+#include "py_igl/py_readTGF.cpp"
 #include "py_igl/py_read_triangle_mesh.cpp"
 #include "py_igl/py_remove_duplicate_vertices.cpp"
 #include "py_igl/py_rotate_vectors.cpp"

python/py_igl/bbw/py_bbw.cpp

@@ -0,0 +1,42 @@
+py::enum_<igl::bbw::QPSolver>(m, "QPSolver")
+    .value("QP_SOLVER_IGL_ACTIVE_SET", igl::bbw::QP_SOLVER_IGL_ACTIVE_SET)
+    .value("QP_SOLVER_MOSEK", igl::bbw::QP_SOLVER_MOSEK)
+    .value("NUM_QP_SOLVERS", igl::bbw::NUM_QP_SOLVERS)
+    .export_values();
+
+// Wrap the BBWData class
+py::class_<igl::bbw::BBWData > BBWData(m, "BBWData");
+
+BBWData
+.def(py::init<>())
+.def_readwrite("partition_unity", &igl::bbw::BBWData::partition_unity)
+.def_readwrite("W0", &igl::bbw::BBWData::W0)
+.def_readwrite("active_set_params", &igl::bbw::BBWData::active_set_params)
+.def_readwrite("qp_solver", &igl::bbw::BBWData::qp_solver)
+.def_readwrite("verbosity", &igl::bbw::BBWData::verbosity)
+#ifndef IGL_NO_MOSEK
+.def_readwrite("mosek_data", &igl::bbw::BBWData::mosek_data)
+#endif
+.def("print", [](igl::bbw::BBWData& data)
+{
+    return data.print();
+})
+;
+
+m.def("bbw", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& Ele,
+  const Eigen::MatrixXi& b,
+  const Eigen::MatrixXd& bc,
+  igl::bbw::BBWData& data,
+  Eigen::MatrixXd& W
+)
+{
+  assert_is_VectorX("b",b);
+  Eigen::VectorXi bv;
+  if (b.size() != 0)
+    bv = b;
+  return igl::bbw::bbw(V, Ele, bv, bc, data, W);
+}, __doc_igl_bbw_bbw,
+py::arg("V"), py::arg("Ele"), py::arg("b"), py::arg("bc"), py::arg("data"), py::arg("W"));

python/py_igl/py_boundary_conditions.cpp

@@ -0,0 +1,24 @@
+
+
+m.def("boundary_conditions", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& Ele,
+  const Eigen::MatrixXd& C,
+  const Eigen::MatrixXi& P,
+  const Eigen::MatrixXi& BE,
+  const Eigen::MatrixXi& CE,
+  Eigen::MatrixXi& b,
+  Eigen::MatrixXd& bc
+)
+{
+  assert_is_VectorX("P", P);
+  Eigen::VectorXi Pv;
+  if (P.size() != 0)
+    Pv = P;
+  Eigen::VectorXi bv;
+  igl::boundary_conditions(V, Ele, C, Pv, BE, CE, bv, bc);
+  b = bv;
+}, __doc_igl_boundary_conditions,
+py::arg("V"), py::arg("Ele"), py::arg("C"), py::arg("P"), py::arg("BE"), py::arg("CE"), py::arg("b"), py::arg("bc"));
+

python/py_igl/py_column_to_quats.cpp

@@ -0,0 +1,10 @@
+m.def("column_to_quats", []
+(
+  const Eigen::MatrixXd& Q,
+  RotationList& vQ
+)
+{
+  return igl::column_to_quats(Q, vQ);
+}, __doc_igl_column_to_quats,
+py::arg("Q"), py::arg("vQ"));
+

python/py_igl/py_deform_skeleton.cpp

@@ -0,0 +1,36 @@
+// COMPLETE BINDINGS ========================
+
+
+m.def("deform_skeleton", []
+(
+  const Eigen::MatrixXd& C,
+  const Eigen::MatrixXi& BE,
+  const Eigen::MatrixXd& T,
+  Eigen::MatrixXd& CT,
+  Eigen::MatrixXi& BET
+)
+{
+  return igl::deform_skeleton(C, BE, T, CT, BET);
+}, __doc_igl_deform_skeleton,
+py::arg("C"), py::arg("BE"), py::arg("T"), py::arg("CT"), py::arg("BET"));
+
+
+
+
+
+// INCOMPLETE BINDINGS ========================
+
+
+//m.def("deform_skeleton", []
+//(
+//  const Eigen::MatrixXd& C,
+//  const Eigen::MatrixXi& BE,
+//  std::vector<Eigen::Affine3d, Eigen::aligned_allocator<Eigen::Affine3d> > & vA,
+//  Eigen::MatrixXd& CT,
+//  Eigen::MatrixXi& BET
+//)
+//{
+//  return igl::deform_skeleton(C, BE, vA, CT, BET);
+//}, __doc_igl_deform_skeleton,
+//py::arg("C"), py::arg("BE"), py::arg("vA"), py::arg("CT"), py::arg("BET"));
+

python/py_igl/py_directed_edge_orientations.cpp

@@ -0,0 +1,12 @@
+
+m.def("directed_edge_orientations", []
+(
+  const Eigen::MatrixXd& C,
+  const Eigen::MatrixXi& E,
+  RotationList& Q
+)
+{
+  return igl::directed_edge_orientations(C, E, Q);
+}, __doc_igl_directed_edge_orientations,
+py::arg("C"), py::arg("E"), py::arg("Q"));
+

python/py_igl/py_directed_edge_parents.cpp

@@ -0,0 +1,14 @@
+
+
+m.def("directed_edge_parents", []
+(
+  const Eigen::MatrixXi& E,
+  Eigen::MatrixXi& P
+)
+{
+  Eigen::VectorXi Pv;
+  igl::directed_edge_parents(E, Pv);
+  P = Pv;
+}, __doc_igl_directed_edge_parents,
+py::arg("E"), py::arg("P"));
+

python/py_igl/py_dqs.cpp

@@ -0,0 +1,21 @@
+
+
+m.def("dqs", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& W,
+  const RotationList& vQ,
+  const std::vector<Eigen::MatrixXd> & vT,
+  Eigen::MatrixXd& U
+)
+{
+  std::vector<Eigen::Vector3d> vTv;
+  for (auto item : vT) {
+    assert_is_Vector3("item", item);
+    Eigen::Vector3d obj = Eigen::Vector3d(item);
+    vTv.push_back(obj);
+  }
+  return igl::dqs(V, W, vQ, vTv, U);
+}, __doc_igl_dqs,
+py::arg("V"), py::arg("W"), py::arg("vQ"), py::arg("vT"), py::arg("U"));
+

python/py_igl/py_forward_kinematics.cpp

@@ -0,0 +1,62 @@
+
+//m.def("forward_kinematics", []
+//(
+//  const Eigen::MatrixXd& C,
+//  const Eigen::MatrixXi& BE,
+//  const Eigen::MatrixXi& P,
+//  std::vector<Eigen::Quaterniond, Eigen::aligned_allocator<Eigen::Quaterniond> > & dQ,
+//  std::vector<Eigen::Vector3d> & dT,
+//  std::vector<Eigen::Quaterniond, Eigen::aligned_allocator<Eigen::Quaterniond> > & vQ,
+//  std::vector<Eigen::Vector3d> & vT
+//)
+//{
+//  return igl::forward_kinematics(C, BE, P, dQ, dT, vQ, vT);
+//}, __doc_igl_forward_kinematics,
+//py::arg("C"), py::arg("BE"), py::arg("P"), py::arg("dQ"), py::arg("dT"), py::arg("vQ"), py::arg("vT"));
+
+m.def("forward_kinematics", []
+(
+  const Eigen::MatrixXd& C,
+  const Eigen::MatrixXi& BE,
+  const Eigen::MatrixXi& P,
+  const RotationList& dQ,
+  RotationList& vQ,
+  py::list vT
+)
+{
+  std::vector<Eigen::Vector3d> vTl;
+  igl::forward_kinematics(C, BE, P, dQ, vQ, vTl);
+  for (auto item : vTl) {
+    py::object obj = py::cast(Eigen::MatrixXd(item));
+    vT.append(obj);
+  }
+}, __doc_igl_forward_kinematics,
+py::arg("C"), py::arg("BE"), py::arg("P"), py::arg("dQ"), py::arg("vQ"), py::arg("vT"));
+
+//m.def("forward_kinematics", []
+//(
+//  const Eigen::MatrixXd& C,
+//  const Eigen::MatrixXi& BE,
+//  const Eigen::MatrixXi& P,
+//  std::vector<Eigen::Quaterniond, Eigen::aligned_allocator<Eigen::Quaterniond> > & dQ,
+//  std::vector<Eigen::Vector3d> & dT,
+//  Eigen::MatrixXd& T
+//)
+//{
+//  return igl::forward_kinematics(C, BE, P, dQ, dT, T);
+//}, __doc_igl_forward_kinematics,
+//py::arg("C"), py::arg("BE"), py::arg("P"), py::arg("dQ"), py::arg("dT"), py::arg("T"));
+
+//m.def("forward_kinematics", []
+//(
+//  const Eigen::MatrixXd& C,
+//  const Eigen::MatrixXi& BE,
+//  const Eigen::MatrixXi& P,
+//  std::vector<Eigen::Quaterniond, Eigen::aligned_allocator<Eigen::Quaterniond> > & dQ,
+//  Eigen::MatrixXd& T
+//)
+//{
+//  return igl::forward_kinematics(C, BE, P, dQ, T);
+//}, __doc_igl_forward_kinematics,
+//py::arg("C"), py::arg("BE"), py::arg("P"), py::arg("dQ"), py::arg("T"));
+

python/py_igl/py_lbs_matrix.cpp

@@ -0,0 +1,67 @@
+// COMPLETE BINDINGS ========================
+
+
+m.def("lbs_matrix", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& W,
+  Eigen::MatrixXd& M
+)
+{
+  return igl::lbs_matrix(V, W, M);
+}, __doc_igl_lbs_matrix,
+py::arg("V"), py::arg("W"), py::arg("M"));
+
+m.def("lbs_matrix_column", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& W,
+  Eigen::MatrixXd& M
+)
+{
+  return igl::lbs_matrix_column(V, W, M);
+}, __doc_igl_lbs_matrix_column,
+py::arg("V"), py::arg("W"), py::arg("M"));
+
+m.def("lbs_matrix_column", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& W,
+  const Eigen::MatrixXi& WI,
+  Eigen::MatrixXd& M
+)
+{
+  return igl::lbs_matrix_column(V, W, WI, M);
+}, __doc_igl_lbs_matrix_column,
+py::arg("V"), py::arg("W"), py::arg("WI"), py::arg("M"));
+
+
+
+
+
+// INCOMPLETE BINDINGS ========================
+
+
+m.def("lbs_matrix_column", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& W,
+  Eigen::SparseMatrix<double>& M
+)
+{
+  return igl::lbs_matrix_column(V, W, M);
+}, __doc_igl_lbs_matrix_column,
+py::arg("V"), py::arg("W"), py::arg("M"));
+
+m.def("lbs_matrix_column", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXd& W,
+  const Eigen::MatrixXi& WI,
+  Eigen::SparseMatrix<double>& M
+)
+{
+  return igl::lbs_matrix_column(V, W, WI, M);
+}, __doc_igl_lbs_matrix_column,
+py::arg("V"), py::arg("W"), py::arg("WI"), py::arg("M"));
+

python/py_igl/py_normalize_row_lengths.cpp

@@ -0,0 +1,12 @@
+
+
+m.def("normalize_row_lengths", []
+(
+  const Eigen::MatrixXd& A,
+  Eigen::MatrixXd& B
+)
+{
+  return igl::normalize_row_lengths(A, B);
+}, __doc_igl_normalize_row_lengths,
+py::arg("A"), py::arg("B"));
+

python/py_igl/py_normalize_row_sums.cpp

@@ -0,0 +1,12 @@
+
+
+m.def("normalize_row_sums", []
+(
+  const Eigen::MatrixXd& A,
+  Eigen::MatrixXd& B
+)
+{
+  return igl::normalize_row_sums(A, B);
+}, __doc_igl_normalize_row_sums,
+py::arg("A"), py::arg("B"));
+

python/py_igl/py_readTGF.cpp

@@ -0,0 +1,54 @@
+// COMPLETE BINDINGS ========================
+
+
+m.def("readTGF", []
+(
+  const std::string tgf_filename,
+  Eigen::MatrixXd& C,
+  Eigen::MatrixXi& E,
+  Eigen::MatrixXi& P,
+  Eigen::MatrixXi& BE,
+  Eigen::MatrixXi& CE,
+  Eigen::MatrixXi& PE
+)
+{
+  Eigen::VectorXi Pv;
+  bool ret = igl::readTGF(tgf_filename, C, E, Pv, BE, CE, PE);
+  P = Pv;
+  return ret;
+}, __doc_igl_readTGF,
+py::arg("tgf_filename"), py::arg("C"), py::arg("E"), py::arg("P"), py::arg("BE"), py::arg("CE"), py::arg("PE"));
+
+m.def("readTGF", []
+(
+  const std::string tgf_filename,
+  Eigen::MatrixXd& C,
+  Eigen::MatrixXi& E
+)
+{
+  return igl::readTGF(tgf_filename, C, E);
+}, __doc_igl_readTGF,
+py::arg("tgf_filename"), py::arg("C"), py::arg("E"));
+
+
+
+
+
+// INCOMPLETE BINDINGS ========================
+
+
+//m.def("readTGF", []
+//(
+//  const std::string tgf_filename,
+//  std::vector<std::vector<double> > & C,
+//  std::vector<std::vector<int> > & E,
+//  std::vector<int> & P,
+//  std::vector<std::vector<int> > & BE,
+//  std::vector<std::vector<int> > & CE,
+//  std::vector<std::vector<int> > & PE
+//)
+//{
+//  return igl::readTGF(tgf_filename, C, E, P, BE, CE, PE);
+//}, __doc_igl_readTGF,
+//py::arg("tgf_filename"), py::arg("C"), py::arg("E"), py::arg("P"), py::arg("BE"), py::arg("CE"), py::arg("PE"));
+

python/python_shared.cpp

@@ -30,12 +30,16 @@ extern void python_export_igl_triangle(py::module &);
 extern void python_export_igl_cgal(py::module &);
 #endif
 
+#ifdef PY_COPYLEFT
+extern void python_export_igl_copyleft(py::module &);
+#endif
+
 #ifdef PY_PNG
 extern void python_export_igl_png(py::module &);
 #endif
 
-#ifdef PY_COPYLEFT
-extern void python_export_igl_copyleft(py::module &);
+#ifdef PY_BBW
+extern void python_export_igl_bbw(py::module &);
 #endif
 
 PYBIND11_PLUGIN(pyigl) {
@@ -58,17 +62,20 @@ PYBIND11_PLUGIN(pyigl) {
            barycenter
            barycentric_coordinates
            barycentric_to_global
+           bbw_bbw
+           boundary_conditions
            boundary_facets
            boundary_loop
            cat
            collapse_edge
            colon
+           column_to_quats
            comb_cross_field
            comb_frame_field
            compute_frame_field_bisectors
+           copyleft_cgal_RemeshSelfIntersectionsParam
            copyleft_cgal_mesh_boolean
            copyleft_cgal_remesh_self_intersections
-           copyleft_cgal_RemeshSelfIntersectionsParam
            copyleft_comiso_miq
            copyleft_comiso_nrosy
            copyleft_marching_cubes
@@ -78,7 +85,11 @@ PYBIND11_PLUGIN(pyigl) {
            covariance_scatter_matrix
            cross_field_missmatch
            cut_mesh_from_singularities
+           deform_skeleton
+           directed_edge_orientations
+           directed_edge_parents
            doublearea
+           dqs
            edge_lengths
            edge_topology
            eigs
@@ -87,6 +98,7 @@ PYBIND11_PLUGIN(pyigl) {
            find_cross_field_singularities
            fit_rotations
            floor
+           forward_kinematics
            gaussian_curvature
            get_seconds
            grad
@@ -96,12 +108,15 @@ PYBIND11_PLUGIN(pyigl) {
            invert_diag
            is_irregular_vertex
            jet
+           lbs_matrix
            local_basis
            lscm
            map_vertices_to_circle
            massmatrix
            min_quad_with_fixed
            n_polyvector
+           normalize_row_lengths
+           normalize_row_sums
            parula
            per_corner_normals
            per_edge_normals
@@ -119,6 +134,7 @@ PYBIND11_PLUGIN(pyigl) {
            readMESH
            readOBJ
            readOFF
+           readTGF
            read_triangle_mesh
            remove_duplicate_vertices
            rotate_vectors
@@ -129,8 +145,7 @@ PYBIND11_PLUGIN(pyigl) {
            slice_mask
            slice_tets
            sortrows
-           streamlines_init
-           streamlines_next
+           streamlines
            triangle_triangle_adjacency
            triangle_triangulate
            unique
@@ -170,12 +185,16 @@ PYBIND11_PLUGIN(pyigl) {
     python_export_igl_cgal(m);
     #endif
 
+    #ifdef PY_COPYLEFT
+    python_export_igl_copyleft(m);
+    #endif
+
     #ifdef PY_PNG
     python_export_igl_png(m);
     #endif
 
-    #ifdef PY_COPYLEFT
-    python_export_igl_copyleft(m);
+    #ifdef PY_BBW
+    python_export_igl_bbw(m);
     #endif
 
     return m.ptr();

python/scripts/py_igl.mako

@@ -1,6 +1,7 @@
 #include <Eigen/Dense>
 
 #include "python_shared.h"
+#include "modules/py_typedefs.h"
 
 % for f in functions:
 #include <igl/${f}.h>
@@ -9,6 +10,8 @@
 
 void python_export_igl(py::module &m)
 {
+#include "modules/py_typedefs.cpp"
+
 % for f in functions:
 #include "py_igl/py_${f}.cpp"
 % endfor

python/scripts/python_shared.mako

@@ -30,12 +30,16 @@ extern void python_export_igl_triangle(py::module &);
 extern void python_export_igl_cgal(py::module &);
 #endif
 
+#ifdef PY_COPYLEFT
+extern void python_export_igl_copyleft(py::module &);
+#endif
+
 #ifdef PY_PNG
 extern void python_export_igl_png(py::module &);
 #endif
 
-#ifdef PY_COPYLEFT
-extern void python_export_igl_copyleft(py::module &);
+#ifdef PY_BBW
+extern void python_export_igl_bbw(py::module &);
 #endif
 
 PYBIND11_PLUGIN(pyigl) {
@@ -82,12 +86,16 @@ PYBIND11_PLUGIN(pyigl) {
     python_export_igl_cgal(m);
     #endif
 
+    #ifdef PY_COPYLEFT
+    python_export_igl_copyleft(m);
+    #endif
+
     #ifdef PY_PNG
     python_export_igl_png(m);
     #endif
 
-    #ifdef PY_COPYLEFT
-    python_export_igl_copyleft(m);
+    #ifdef PY_BBW
+    python_export_igl_bbw(m);
     #endif
 
     return m.ptr();

python/tutorial/402_PolyharmonicDeformation.py

@@ -93,8 +93,8 @@ viewer = igl.viewer.Viewer()
 viewer.data.set_mesh(U, F)
 viewer.core.show_lines = False
 viewer.data.set_colors(C)
-# viewer.core.trackball_angle = igl.eigen.Quaterniond(0.81,-0.58,-0.03,-0.03)
-# viewer.core.trackball_angle.normalize()
+viewer.core.trackball_angle = igl.eigen.Quaterniond(0.81,-0.58,-0.03,-0.03)
+viewer.core.trackball_angle.normalize()
 viewer.callback_pre_draw = pre_draw
 viewer.callback_key_down = key_down
 viewer.core.is_animating = True

python/tutorial/403_BoundedBiharmonicWeights.py

@@ -0,0 +1,145 @@
+import sys, os
+
+# 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, print_usage
+
+dependencies = ["viewer", "bbw"]
+check_dependencies(dependencies)
+
+
+def pre_draw(viewer):
+    global pose, anim_t, C, BE, P, U, M, anim_t_dir
+
+    if viewer.core.is_animating:
+        # Interpolate pose and identity
+        anim_pose = igl.RotationList(len(pose))
+
+        for e in range(len(pose)):
+            anim_pose[e] = pose[e].slerp(anim_t, igl.eigen.Quaterniond.Identity())
+
+        # Propogate relative rotations via FK to retrieve absolute transformations
+        vQ = igl.RotationList()
+        vT = []
+        igl.forward_kinematics(C, BE, P, anim_pose, vQ, vT)
+        dim = C.cols()
+        T = igl.eigen.MatrixXd(BE.rows() * (dim + 1), dim)
+        for e in range(BE.rows()):
+            a = igl.eigen.Affine3d.Identity()
+            a.translate(vT[e])
+            a.rotate(vQ[e])
+            T.setBlock(e * (dim + 1), 0, dim + 1, dim, a.matrix().transpose().block(0, 0, dim + 1, dim))
+
+        # Compute deformation via LBS as matrix multiplication
+        U = M * T
+
+        # Also deform skeleton edges
+        CT = igl.eigen.MatrixXd()
+        BET = igl.eigen.MatrixXi()
+        igl.deform_skeleton(C, BE, T, CT, BET)
+
+        viewer.data.set_vertices(U)
+        viewer.data.set_edges(CT, BET, sea_green)
+        viewer.data.compute_normals()
+        anim_t += anim_t_dir
+        anim_t_dir *= -1.0 if (0.0 >= anim_t or anim_t >= 1.0) else 1.0
+
+    return False
+
+
+def key_down(viewer, key, mods):
+    global selected, W
+    if key == ord('.'):
+        selected += 1
+        selected = min(max(selected, 0), W.cols()-1)
+        set_color(viewer)
+    elif key == ord(','):
+        selected -= 1
+        selected = min(max(selected, 0), W.cols()-1)
+        set_color(viewer)
+    elif key == ord(' '):
+        viewer.core.is_animating = not viewer.core.is_animating
+
+    return True
+
+
+def set_color(viewer):
+    global selected, W
+    C = igl.eigen.MatrixXd()
+    igl.jet(W.col(selected), True, C)
+    viewer.data.set_colors(C)
+
+
+if __name__ == "__main__":
+    keys = {".": "show next weight function",
+            ",": "show previous weight function",
+            "space": "toggle animation"}
+
+    print_usage(keys)
+
+    V = igl.eigen.MatrixXd()
+    W = igl.eigen.MatrixXd()
+    U = igl.eigen.MatrixXd()
+    C = igl.eigen.MatrixXd()
+    M = igl.eigen.MatrixXd()
+    Q = igl.eigen.MatrixXd()
+    T = igl.eigen.MatrixXi()
+    F = igl.eigen.MatrixXi()
+    BE = igl.eigen.MatrixXi()
+    P = igl.eigen.MatrixXi()
+
+    sea_green = igl.eigen.MatrixXd([[70. / 255., 252. / 255., 167. / 255.]])
+
+    selected = 0
+    pose = igl.RotationList()
+    anim_t = 1.0
+    anim_t_dir = -0.03
+
+    igl.readMESH(TUTORIAL_SHARED_PATH + "hand.mesh", V, T, F)
+    U = igl.eigen.MatrixXd(V)
+    igl.readTGF(TUTORIAL_SHARED_PATH + "hand.tgf", C, BE)
+
+    # Retrieve parents for forward kinematics
+    igl.directed_edge_parents(BE, P)
+
+    # Read pose as matrix of quaternions per row
+    igl.readDMAT(TUTORIAL_SHARED_PATH + "hand-pose.dmat", Q)
+    igl.column_to_quats(Q, pose)
+    assert (len(pose) == BE.rows())
+
+    # List of boundary indices (aka fixed value indices into VV)
+    b = igl.eigen.MatrixXi()
+    # List of boundary conditions of each weight function
+    bc = igl.eigen.MatrixXd()
+
+    igl.boundary_conditions(V, T, C, igl.eigen.MatrixXi(), BE, igl.eigen.MatrixXi(), b, bc)
+
+    # compute BBW weights matrix
+    bbw_data = igl.bbw.BBWData()
+    # only a few iterations for sake of demo
+    bbw_data.active_set_params.max_iter = 8
+    bbw_data.verbosity = 2
+    if not igl.bbw.bbw(V, T, b, bc, bbw_data, W):
+        exit(-1)
+
+    # Normalize weights to sum to one
+    igl.normalize_row_sums(W, W)
+    # precompute linear blend skinning matrix
+    igl.lbs_matrix(V, W, M)
+
+    # Plot the mesh with pseudocolors
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(U, F)
+    set_color(viewer)
+    viewer.data.set_edges(C, BE, sea_green)
+    viewer.core.show_lines = False
+    viewer.core.show_overlay_depth = False
+    viewer.core.line_width = 1
+    viewer.core.trackball_angle.normalize()
+    viewer.callback_pre_draw = pre_draw
+    viewer.callback_key_down = key_down
+    viewer.core.is_animating = False
+    viewer.core.animation_max_fps = 30.0
+    viewer.launch()

python/tutorial/404_DualQuaternionSkinning.py

@@ -0,0 +1,140 @@
+import sys, os
+from math import sin, cos, pi
+
+# Add the igl library to the modules search path
+import math
+
+sys.path.insert(0, os.getcwd() + "/../")
+import pyigl as igl
+
+from shared import TUTORIAL_SHARED_PATH, check_dependencies, print_usage
+
+dependencies = ["viewer"]
+check_dependencies(dependencies)
+
+
+def pre_draw(viewer):
+    global recompute, anim_t, poses, C, BE, P, U, M, anim_t_dir
+
+    if recompute:
+        # Find pose interval
+        begin = int(math.floor(anim_t)) % len(poses)
+        end = int(math.floor(anim_t) + 1) % len(poses)
+        t = anim_t - math.floor(anim_t)
+
+        # Interpolate pose and identity
+        anim_pose = igl.RotationList()
+        for e in range(len(poses[begin])):
+            anim_pose.append(poses[begin][e].slerp(t, poses[end][e]))
+
+        # Propogate relative rotations via FK to retrieve absolute transformations
+        vQ = igl.RotationList()
+        vT = []
+        igl.forward_kinematics(C, BE, P, anim_pose, vQ, vT)
+        dim = C.cols()
+        T = igl.eigen.MatrixXd(BE.rows() * (dim + 1), dim)
+        for e in range(BE.rows()):
+            a = igl.eigen.Affine3d.Identity()
+            a.translate(vT[e])
+            a.rotate(vQ[e])
+            T.setBlock(e * (dim + 1), 0, dim + 1, dim, a.matrix().transpose().block(0, 0, dim + 1, dim))
+
+        # Compute deformation via LBS as matrix multiplication
+        if use_dqs:
+            igl.dqs(V, W, vQ, vT, U)
+        else:
+            U = M * T
+
+        # Also deform skeleton edges
+        CT = igl.eigen.MatrixXd()
+        BET = igl.eigen.MatrixXi()
+        igl.deform_skeleton(C, BE, T, CT, BET)
+
+        viewer.data.set_vertices(U)
+        viewer.data.set_edges(CT, BET, sea_green)
+        viewer.data.compute_normals()
+        if viewer.core.is_animating:
+            anim_t += anim_t_dir
+        else:
+            recompute = False
+
+    return False
+
+
+def key_down(viewer, key, mods):
+    global recompute, use_dqs, animation
+    recompute = True
+    if key == ord('D') or key == ord('d'):
+        use_dqs = not use_dqs
+        viewer.core.is_animating = False
+        animation = False
+        if use_dqs:
+            print("Switched to Dual Quaternion Skinning")
+        else:
+            print("Switched to Linear Blend Skinning")
+    elif key == ord(' '):
+        if animation:
+            viewer.core.is_animating = False
+            animation = False
+        else:
+            viewer.core.is_animating = True
+            animation = True
+    return False
+
+
+if __name__ == "__main__":
+    keys = {"d": "toggle between LBS and DQS",
+            "space": "toggle animation"}
+
+    print_usage(keys)
+
+    V = igl.eigen.MatrixXd()
+    F = igl.eigen.MatrixXi()
+    C = igl.eigen.MatrixXd()
+    BE = igl.eigen.MatrixXi()
+    P = igl.eigen.MatrixXi()
+    W = igl.eigen.MatrixXd()
+    M = igl.eigen.MatrixXd()
+
+    sea_green = igl.eigen.MatrixXd([[70. / 255., 252. / 255., 167. / 255.]])
+
+    anim_t = 0.0
+    anim_t_dir = 0.015
+    use_dqs = False
+    recompute = True
+    animation = False  # Flag needed as there is some synchronization problem with viewer.core.is_animating
+
+    poses = [[]]
+
+    igl.readOBJ(TUTORIAL_SHARED_PATH + "arm.obj", V, F)
+    U = igl.eigen.MatrixXd(V)
+    igl.readTGF(TUTORIAL_SHARED_PATH + "arm.tgf", C, BE)
+
+    # retrieve parents for forward kinematics
+    igl.directed_edge_parents(BE, P)
+    rest_pose = igl.RotationList()
+    igl.directed_edge_orientations(C, BE, rest_pose)
+    poses = [[igl.eigen.Quaterniond.Identity() for i in range(4)] for j in range(4)]
+
+    twist = igl.eigen.Quaterniond(pi, igl.eigen.MatrixXd([1, 0, 0]))
+    poses[1][2] = rest_pose[2] * twist * rest_pose[2].conjugate()
+    bend = igl.eigen.Quaterniond(-pi * 0.7, igl.eigen.MatrixXd([0, 0, 1]))
+    poses[3][2] = rest_pose[2] * bend * rest_pose[2].conjugate()
+
+    igl.readDMAT(TUTORIAL_SHARED_PATH + "arm-weights.dmat", W)
+    igl.lbs_matrix(V, W, M)
+
+    # Plot the mesh with pseudocolors
+    viewer = igl.viewer.Viewer()
+    viewer.data.set_mesh(U, F)
+    viewer.data.set_edges(C, BE, sea_green)
+    viewer.core.show_lines = False
+    viewer.core.show_overlay_depth = False
+    viewer.core.line_width = 1
+    viewer.core.trackball_angle.normalize()
+    viewer.callback_pre_draw = pre_draw
+    viewer.callback_key_down = key_down
+    viewer.core.is_animating = False
+    viewer.core.camera_zoom = 2.5
+    viewer.core.animation_max_fps = 30.0
+    viewer.launch()

shared/cmake/CMakeLists.txt

@@ -387,6 +387,7 @@ endif()
 
 ### Compile the viewer ###
 if(LIBIGL_WITH_VIEWER)
+  find_package(OpenGL REQUIRED)
   set(NANOGUI_DIR "${LIBIGL_EXTERNAL}/nanogui")
 
   if(LIBIGL_WITH_NANOGUI)
@@ -415,7 +416,7 @@ if(LIBIGL_WITH_VIEWER)
     add_subdirectory("${NANOGUI_DIR}/ext/glfw" "glfw")
 
     set(VIEWER_INCLUDE_DIRS "${NANOGUI_DIR}/ext/glfw/include")
-    set(LIBIGL_VIEWER_EXTRA_LIBRARIES "glfw" ${GLFW_LIBRARIES})
+    set(LIBIGL_VIEWER_EXTRA_LIBRARIES "glfw" ${OPENGL_gl_LIBRARY})
   endif()
 
   ### GLEW for linux and windows

tutorial/611_SLIM/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 2.8.12)
+project(611_SLIM)
+
+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})

tutorial/611_SLIM/main.cpp

@@ -0,0 +1,394 @@
+#include <iostream>
+
+#include "igl/slim.h"
+
+#include "igl/components.h"
+#include "igl/readOBJ.h"
+#include "igl/writeOBJ.h"
+#include "igl/Timer.h"
+
+#include "igl/boundary_loop.h"
+#include "igl/map_vertices_to_circle.h"
+#include "igl/harmonic.h"
+#include <igl/serialize.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/flipped_triangles.h>
+#include <igl/euler_characteristic.h>
+#include <igl/barycenter.h>
+#include <igl/adjacency_matrix.h>
+#include <igl/is_edge_manifold.h>
+#include <igl/doublearea.h>
+#include <igl/cat.h>
+
+#include <stdlib.h>
+
+#include <string>
+#include <vector>
+
+using namespace std;
+using namespace Eigen;
+
+void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F);
+void param_2d_demo_iter(igl::viewer::Viewer& viewer);
+void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
+void soft_const_demo_iter(igl::viewer::Viewer& viewer);
+void deform_3d_demo_iter(igl::viewer::Viewer& viewer);
+void get_cube_corner_constraints(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc);
+void display_3d_mesh(igl::viewer::Viewer& viewer);
+void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec);
+
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+bool first_iter = true;
+igl::SLIMData sData;
+igl::Timer timer;
+
+double uv_scale_param;
+
+enum DEMO_TYPE {
+  PARAM_2D,
+  SOFT_CONST,
+  DEFORM_3D
+};
+DEMO_TYPE demo_type;
+
+bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier){
+  if (key == ' ') {
+    switch (demo_type) {
+      case PARAM_2D: {
+        param_2d_demo_iter(viewer);
+        break;
+      }
+      case SOFT_CONST: {
+        soft_const_demo_iter(viewer);
+        break;
+      }
+      case DEFORM_3D: {
+        deform_3d_demo_iter(viewer);
+        break;
+      }
+      default:
+        break;
+    }
+  }
+
+  return false;
+}
+
+void param_2d_demo_iter(igl::viewer::Viewer& viewer) {
+  if (first_iter) {
+    timer.start();
+    igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/face.obj", V, F);
+    check_mesh_for_issues(V,F);
+    cout << "\tMesh is valid!" << endl;
+
+    Eigen::MatrixXd uv_init;
+    Eigen::VectorXi bnd; Eigen::MatrixXd bnd_uv;
+    igl::boundary_loop(F,bnd);
+    igl::map_vertices_to_circle(V,bnd,bnd_uv);
+
+    igl::harmonic(V,F,bnd,bnd_uv,1,uv_init);
+    if (igl::flipped_triangles(uv_init,F).size() != 0) {
+      igl::harmonic(F,bnd,bnd_uv,1,uv_init); // use uniform laplacian
+    }
+
+    cout << "initialized parametrization" << endl;
+
+    sData.slim_energy = igl::SLIMData::SYMMETRIC_DIRICHLET;
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    slim_precompute(V,F,uv_init,sData, igl::SLIMData::SYMMETRIC_DIRICHLET, b,bc,0);
+
+    uv_scale_param = 15 * (1./sqrt(sData.mesh_area));
+    viewer.data.set_mesh(V, F);
+    viewer.core.align_camera_center(V,F);
+    viewer.data.set_uv(sData.V_o*uv_scale_param);
+    viewer.data.compute_normals();
+    viewer.core.show_texture = true;
+
+    first_iter = false;
+  } else {
+    slim_solve(sData,1); // 1 iter
+    viewer.data.set_uv(sData.V_o*uv_scale_param);
+    cout << "time = " << timer.getElapsedTime() << endl;
+  }
+}
+
+void soft_const_demo_iter(igl::viewer::Viewer& viewer) {
+  if (first_iter) {
+
+    igl::read_triangle_mesh(TUTORIAL_SHARED_PATH "/circle.obj", V, F);
+
+    check_mesh_for_issues(V,F);
+    cout << "\tMesh is valid!" << endl;
+    Eigen::MatrixXd V_0 = V.block(0,0,V.rows(),2);
+
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    get_soft_constraint_for_circle(V_0,F,b,bc);
+    double soft_const_p = 1e5;
+    slim_precompute(V,F,V_0,sData,igl::SLIMData::SYMMETRIC_DIRICHLET,b,bc,soft_const_p);
+
+    viewer.data.set_mesh(V, F);
+    viewer.core.align_camera_center(V,F);
+    viewer.data.compute_normals();
+    viewer.core.show_lines = true;
+
+    first_iter = false;
+
+  } else {
+    slim_solve(sData,1); // 1 iter
+    viewer.data.set_mesh(sData.V_o, F);
+  }
+}
+
+void deform_3d_demo_iter(igl::viewer::Viewer& viewer) {
+  if (first_iter) {
+    igl::readOBJ(TUTORIAL_SHARED_PATH "/cube_40k.obj", V, F);
+
+    Eigen::MatrixXd V_0 = V;
+    Eigen::VectorXi b; Eigen::MatrixXd bc;
+    get_cube_corner_constraints(V_0,F,b,bc);
+
+    double soft_const_p = 1e5;
+    sData.exp_factor = 5.0;
+    slim_precompute(V,F,V_0,sData,igl::SLIMData::EXP_CONFORMAL,b,bc,soft_const_p);
+    cout << "precomputed" << endl;
+
+    first_iter = false;
+    display_3d_mesh(viewer);
+
+  } else {
+    slim_solve(sData,1); // 1 iter
+    display_3d_mesh(viewer);
+  }
+}
+
+void display_3d_mesh(igl::viewer::Viewer& viewer) {
+  MatrixXd V_temp; MatrixXi F_temp;
+  Eigen::MatrixXd Barycenters;
+
+  igl::barycenter(sData.V,sData.F,Barycenters);
+  //cout << "Barycenters.rows() = " << Barycenters.rows() << endl;
+  //double t = double((key - '1')+1) / 9.0;
+  double view_depth = 10.;
+  double t = view_depth/9.;
+
+  VectorXd v = Barycenters.col(2).array() - Barycenters.col(2).minCoeff();
+  v /= v.col(0).maxCoeff();
+
+  vector<int> s;
+
+  for (unsigned i=0; i<v.size();++i)
+    if (v(i) < t)
+      s.push_back(i);
+
+  V_temp.resize(s.size()*4,3);
+  F_temp.resize(s.size()*4,3);
+
+  for (unsigned i=0; i<s.size();++i){
+    V_temp.row(i*4+0) = sData.V_o.row(sData.F(s[i],0));
+    V_temp.row(i*4+1) = sData.V_o.row(sData.F(s[i],1));
+    V_temp.row(i*4+2) = sData.V_o.row(sData.F(s[i],2));
+    V_temp.row(i*4+3) = sData.V_o.row(sData.F(s[i],3));
+    F_temp.row(i*4+0) << (i*4)+0, (i*4)+1, (i*4)+3;
+    F_temp.row(i*4+1) << (i*4)+0, (i*4)+2, (i*4)+1;
+    F_temp.row(i*4+2) << (i*4)+3, (i*4)+2, (i*4)+0;
+    F_temp.row(i*4+3) << (i*4)+1, (i*4)+2, (i*4)+3;
+  }
+  viewer.data.set_mesh(V_temp,F_temp);
+  viewer.core.align_camera_center(V_temp,F_temp);
+  viewer.data.set_face_based(true);
+  viewer.core.show_lines = true;
+}
+
+int main(int argc, char *argv[]) {
+
+  cerr << "Press space for running an iteration." << std::endl;
+  cerr << "Syntax: " << argv[0] << " demo_number (1 to 3)" << std::endl;
+  cerr << "1. 2D unconstrained parametrization" << std::endl;
+  cerr << "2. 2D deformation with positional constraints" << std::endl;
+  cerr << "3. 3D mesh deformation with positional constraints" << std::endl;
+
+  demo_type = PARAM_2D;
+
+   if (argc == 2) {
+     switch (std::atoi(argv[1])) {
+       case 1: {
+         demo_type = PARAM_2D;
+         break;
+       } case 2: {
+         demo_type = SOFT_CONST;
+         break;
+       } case 3: {
+         demo_type = DEFORM_3D;
+         break;
+       }
+       default: {
+         cerr << "Wrong demo number - Please choose one between 1-3" << std:: endl;
+         exit(1);
+       }
+     }
+   }
+
+
+  // Launch the viewer
+  igl::viewer::Viewer viewer;
+  viewer.callback_key_down = &key_down;
+
+  // Disable wireframe
+  viewer.core.show_lines = false;
+
+  // Draw checkerboard texture
+  viewer.core.show_texture = false;
+
+  // First iteration
+  key_down(viewer, ' ', 0);
+
+  viewer.launch();
+
+  return 0;
+}
+
+void check_mesh_for_issues(Eigen::MatrixXd& V, Eigen::MatrixXi& F) {
+
+  Eigen::SparseMatrix<double> A;
+  igl::adjacency_matrix(F,A);
+
+  Eigen::MatrixXi C, Ci;
+  igl::components(A, C, Ci);
+
+  int connected_components = Ci.rows();
+  if (connected_components!=1) {
+    cout << "Error! Input has multiple connected components" << endl; exit(1);
+  }
+  int euler_char = igl::euler_characteristic(V, F);
+  if (!euler_char) {
+    cout << "Error! Input does not have a disk topology, it's euler char is " << euler_char << endl; exit(1);
+  }
+  bool is_edge_manifold = igl::is_edge_manifold(F);
+  if (!is_edge_manifold) {
+    cout << "Error! Input is not an edge manifold" << endl; exit(1);
+  }
+
+  Eigen::VectorXd areas; igl::doublearea(V,F,areas);
+  const double eps = 1e-14;
+  for (int i = 0; i < areas.rows(); i++) {
+    if (areas(i) < eps) {
+      cout << "Error! Input has zero area faces" << endl; exit(1);
+    }
+  }
+}
+
+void get_soft_constraint_for_circle(Eigen::MatrixXd& V_o, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
+
+    Eigen::VectorXi bnd;
+    igl::boundary_loop(F,bnd);
+    const int B_STEPS = 22; // constraint every B_STEPS vertices of the boundary
+
+    b.resize(bnd.rows()/B_STEPS);
+    bc.resize(b.rows(),2);
+
+    int c_idx = 0;
+    for (int i = B_STEPS; i < bnd.rows(); i+=B_STEPS) {
+        b(c_idx) = bnd(i);
+        c_idx++;
+    }
+
+    bc.row(0) = V_o.row(b(0)); // keep it there for now
+    bc.row(1) = V_o.row(b(2));
+    bc.row(2) = V_o.row(b(3));
+    bc.row(3) = V_o.row(b(4));
+    bc.row(4) = V_o.row(b(5));
+
+
+    bc.row(0) << V_o(b(0),0), 0;
+    bc.row(4) << V_o(b(4),0), 0;
+    bc.row(2) << V_o(b(2),0), 0.1;
+    bc.row(3) << V_o(b(3),0), 0.05;
+    bc.row(1) << V_o(b(1),0), -0.15;
+    bc.row(5) << V_o(b(5),0), +0.15;
+}
+
+void get_cube_corner_constraints(Eigen::MatrixXd& V, Eigen::MatrixXi& F, Eigen::VectorXi& b, Eigen::MatrixXd& bc) {
+  double min_x,max_x,min_y,max_y,min_z,max_z;
+  min_x = V.col(0).minCoeff(); max_x = V.col(0).maxCoeff();
+  min_y = V.col(1).minCoeff(); max_y = V.col(1).maxCoeff();
+  min_z = V.col(2).minCoeff(); max_z = V.col(2).maxCoeff();
+
+
+  // get all cube corners
+  b.resize(8,1); bc.resize(8,3);
+  int x;
+  for (int i = 0; i < V.rows(); i++) {
+    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,min_z)) b(0) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,min_y,max_z)) b(1) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,min_z)) b(2) = i;
+    if (V.row(i) == Eigen::RowVector3d(min_x,max_y,max_z)) b(3) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,min_z)) b(4) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,min_z)) b(5) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,min_y,max_z)) b(6) = i;
+    if (V.row(i) == Eigen::RowVector3d(max_x,max_y,max_z)) b(7) = i;
+  }
+
+  // get all cube edges
+  std::set<int> cube_edge1; Eigen::VectorXi cube_edge1_vec;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == min_x && V(i,1) == min_y)) {
+      cube_edge1.insert(i);
+    }
+  }
+  Eigen::VectorXi edge1;
+  int_set_to_eigen_vector(cube_edge1, edge1);
+
+  std::set<int> cube_edge2; Eigen::VectorXi edge2;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == max_x && V(i,1) == max_y)) {
+      cube_edge2.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge2, edge2);
+  b = igl::cat(1,edge1,edge2);
+
+  std::set<int> cube_edge3; Eigen::VectorXi edge3;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == max_x && V(i,1) == min_y)) {
+      cube_edge3.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge3, edge3);
+  b = igl::cat(1,b,edge3);
+
+  std::set<int> cube_edge4; Eigen::VectorXi edge4;
+  for (int i = 0; i < V.rows(); i++) {
+    if ((V(i,0) == min_x && V(i,1) == max_y)) {
+      cube_edge4.insert(i);
+    }
+  }
+  int_set_to_eigen_vector(cube_edge4, edge4);
+  b = igl::cat(1,b,edge4);
+
+  bc.resize(b.rows(),3);
+  Eigen::Matrix3d m; m = Eigen::AngleAxisd(0.3 * M_PI, Eigen::Vector3d(1./sqrt(2.),1./sqrt(2.),0.)/*Eigen::Vector3d::UnitX()*/);
+  int i = 0;
+  for (; i < cube_edge1.size(); i++) {
+    Eigen::RowVector3d edge_rot_center(min_x,min_y,(min_z+max_z)/2.);
+    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m + edge_rot_center;
+  }
+  for (; i < cube_edge1.size() + cube_edge2.size(); i++) {
+    Eigen::RowVector3d edge_rot_center(max_x,max_y,(min_z+max_z)/2.);
+    bc.row(i) = (V.row(b(i)) - edge_rot_center) * m.transpose() + edge_rot_center;
+  }
+  for (; i < cube_edge1.size() + cube_edge2.size() + cube_edge3.size(); i++) {
+    bc.row(i) = 0.75*V.row(b(i));
+  }
+  for (; i < b.rows(); i++) {
+    bc.row(i) = 0.75*V.row(b(i));
+  }
+}
+
+void int_set_to_eigen_vector(const std::set<int>& int_set, Eigen::VectorXi& vec) {
+  vec.resize(int_set.size()); int idx = 0;
+  for(auto f : int_set) {
+      vec(idx) = f; idx++;
+    }
+}

tutorial/CMakeLists.txt

@@ -167,6 +167,7 @@ if(TUTORIALS_CHAPTER6)
     add_subdirectory("609_Boolean")
     add_subdirectory("610_CSGTree")
   endif()
+  add_subdirectory("611_SLIM")
 endif()
 
 # Chapter 7

tutorial/images/611_SLIM.png.REMOVED.git-id

@@ -0,0 +1 @@
+9fcca7386e32522dc9425a7722ab7d7adb8c3f56

tutorial/shared/circle.obj.REMOVED.git-id

@@ -0,0 +1 @@
+8765e0b494718013bdd1305490a1d60d637f1b61

tutorial/shared/cube_40k.obj.REMOVED.git-id

@@ -0,0 +1 @@
+15badea46771f5da1ce678b5321d01b84c1ac6da

tutorial/shared/face.obj.REMOVED.git-id

@@ -0,0 +1 @@
+b07bc4c5241330e8bc15cef5d50c166e4346240c

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-b4d291ea3ce115a6b9af9ffe3559ffd5f157bea6
+cd11b42bda4357cfaae8f62d5036f840df5b627d

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-0e9dc3a10233852a045d53b248481da495b37e4c
+91d360c45d833747534b3dd524f5d2b39a1e73b4