Merge remote-tracking branch 'upstream/master'

Former-commit-id: 118b679308fb5321f32df8b5fb66276cf51bfbe7

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

include/igl/AABB.cpp

@@ -391,7 +391,7 @@ igl::AABB<DerivedV,DIM>::squared_distance(
       int i_left;
       RowVectorDIMS c_left = c;
       Scalar sqr_d_left = m_left->squared_distance(V,Ele,p,sqr_d,i_left,c_left);
-      set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
+      this->set_min(p,sqr_d_left,i_left,c_left,sqr_d,i,c);
       looked_left = true;
     };
     const auto & look_right = [&]()
@@ -399,7 +399,7 @@ igl::AABB<DerivedV,DIM>::squared_distance(
       int i_right;
       RowVectorDIMS c_right = c;
       Scalar sqr_d_right = m_right->squared_distance(V,Ele,p,sqr_d,i_right,c_right);
-      set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
+      this->set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
       looked_right = true;
     };
 

include/igl/bbw.cpp

@@ -0,0 +1,143 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "bbw.h"
+#include "min_quad_with_fixed.h"
+#include "harmonic.h"
+#include "parallel_for.h"
+#include <Eigen/Sparse>
+#include <iostream>
+#include <mutex>
+#include <cstdio>
+
+igl::BBWData::BBWData():
+  partition_unity(false),
+  W0(),
+  active_set_params(),
+  verbosity(0)
+{
+  // We know that the Bilaplacian is positive semi-definite
+  active_set_params.Auu_pd = true;
+}
+
+void igl::BBWData::print()
+{
+  using namespace std;
+  cout<<"partition_unity: "<<partition_unity<<endl;
+  cout<<"W0=["<<endl<<W0<<endl<<"];"<<endl;
+}
+
+
+template <
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedW>
+IGL_INLINE bool igl::bbw(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedEle> & Ele,
+  const Eigen::PlainObjectBase<Derivedb> & b,
+  const Eigen::PlainObjectBase<Derivedbc> & bc,
+  igl::BBWData & data,
+  Eigen::PlainObjectBase<DerivedW> & W
+  )
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(!data.partition_unity && "partition_unity not implemented yet");
+  // number of domain vertices
+  int n = V.rows();
+  // number of handles
+  int m = bc.cols();
+  // Build biharmonic operator
+  Eigen::SparseMatrix<typename DerivedV::Scalar> Q;
+  harmonic(V,Ele,2,Q);
+  W.derived().resize(n,m);
+  // No linear terms
+  VectorXd c = VectorXd::Zero(n);
+  // No linear constraints
+  SparseMatrix<typename DerivedW::Scalar> A(0,n),Aeq(0,n),Aieq(0,n);
+  VectorXd Beq(0,1),Bieq(0,1);
+  // Upper and lower box constraints (Constant bounds)
+  VectorXd ux = VectorXd::Ones(n);
+  VectorXd lx = VectorXd::Zero(n);
+  active_set_params eff_params = data.active_set_params;
+  if(data.verbosity >= 1)
+  {
+    cout<<"BBW: max_iter: "<<data.active_set_params.max_iter<<endl;
+    cout<<"BBW: eff_max_iter: "<<eff_params.max_iter<<endl;
+  }
+  if(data.verbosity >= 1)
+  {
+    cout<<"BBW: Computing initial weights for "<<m<<" handle"<<
+      (m!=1?"s":"")<<"."<<endl;
+  }
+  min_quad_with_fixed_data<typename DerivedW::Scalar > mqwf;
+  min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
+  min_quad_with_fixed_solve(mqwf,c,bc,Beq,W);
+  // decrement
+  eff_params.max_iter--;
+  bool error = false;
+  // Loop over handles
+  std::mutex critical;
+  const auto & optimize_weight = [&](const int i)
+  {
+    // Quicker exit for paralle_for
+    if(error)
+    {
+      return;
+    }
+    if(data.verbosity >= 1)
+    {
+      std::lock_guard<std::mutex> lock(critical);
+      cout<<"BBW: Computing weight for handle "<<i+1<<" out of "<<m<<
+        "."<<endl;
+    }
+    VectorXd bci = bc.col(i);
+    VectorXd Wi;
+    // use initial guess
+    Wi = W.col(i);
+    SolverStatus ret = active_set(
+        Q,c,b,bci,Aeq,Beq,Aieq,Bieq,lx,ux,eff_params,Wi);
+    switch(ret)
+    {
+      case SOLVER_STATUS_CONVERGED:
+        break;
+      case SOLVER_STATUS_MAX_ITER:
+        cerr<<"active_set: max iter without convergence."<<endl;
+        break;
+      case SOLVER_STATUS_ERROR:
+      default:
+        cerr<<"active_set error."<<endl;
+        error = true;
+    }
+    W.col(i) = Wi;
+  };
+  parallel_for(m,optimize_weight,2);
+  if(error)
+  {
+    return false;
+  }
+
+#ifndef NDEBUG
+  const double min_rowsum = W.rowwise().sum().array().abs().minCoeff();
+  if(min_rowsum < 0.1)
+  {
+    cerr<<"bbw.cpp: Warning, minimum row sum is very low. Consider more "
+      "active set iterations or enforcing partition of unity."<<endl;
+  }
+#endif
+
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::bbw<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::BBWData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif
+

include/igl/bbw.h

@@ -0,0 +1,77 @@
+// 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_BBW_H
+#define IGL_BBW_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <igl/active_set.h>
+
+namespace igl
+{
+  // Container for BBW computation related data and flags
+  class BBWData
+  {
+    public:
+      // Enforce partition of unity during optimization (optimize all weight
+      // simultaneously)
+      bool partition_unity;
+      // Initial guess
+      Eigen::MatrixXd W0;
+      igl::active_set_params active_set_params;
+      // Verbosity level
+      // 0: quiet
+      // 1: loud
+      // 2: louder
+      int verbosity;
+    public:
+      IGL_INLINE BBWData();
+      // Print current state of object
+      IGL_INLINE void print();
+  };
+
+  // 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
+  template <
+    typename DerivedV,
+    typename DerivedEle,
+    typename Derivedb,
+    typename Derivedbc,
+    typename DerivedW>
+  IGL_INLINE bool bbw(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedEle> & Ele,
+    const Eigen::PlainObjectBase<Derivedb> & b,
+    const Eigen::PlainObjectBase<Derivedbc> & bc,
+    BBWData & data,
+    Eigen::PlainObjectBase<DerivedW> & W);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bbw.cpp"
+#endif
+
+#endif
+

include/igl/bbw/bbw.cpp

@@ -64,29 +64,17 @@ IGL_INLINE bool igl::bbw::bbw(
   int n = V.rows();
   // number of handles
   int m = bc.cols();
-
+  // Build biharmonic operator
   SparseMatrix<typename DerivedW::Scalar> L;
   cotmatrix(V,Ele,L);
-  MassMatrixType mmtype = MASSMATRIX_TYPE_VORONOI;
-  if(Ele.cols() == 4)
-  {
-    mmtype = MASSMATRIX_TYPE_BARYCENTRIC;
-  }
   SparseMatrix<typename DerivedW::Scalar> M;
   SparseMatrix<typename DerivedW::Scalar> Mi;
-  massmatrix(V,Ele,mmtype,M);
-
+  massmatrix(V,Ele,MASSMATRIX_TYPE_DEFAULT,M);
   invert_diag(M,Mi);
-
-  // Biharmonic operator
   SparseMatrix<typename DerivedW::Scalar> Q = L.transpose() * Mi * L;
+  assert(!data.partition_unity && "partition_unity not implemented yet");
 
   W.derived().resize(n,m);
-  if(data.partition_unity)
-  {
-    // Not yet implemented
-    assert(false);
-  }else
   {
     // No linear terms
     VectorXd c = VectorXd::Zero(n);

include/igl/boundary_loop.cpp

@@ -114,6 +114,13 @@ IGL_INLINE void igl::boundary_loop(
     }
   }
 
+  //Check for meshes without boundary
+  if (idxMax == -1)
+  {
+      L.clear();
+      return;
+  }
+
   L.resize(Lall[idxMax].size());
   for (size_t i = 0; i < Lall[idxMax].size(); ++i)
   {

include/igl/comb_frame_field.cpp

@@ -49,7 +49,7 @@ IGL_INLINE void igl::comb_frame_field(const Eigen::PlainObjectBase<DerivedV> &V,
     {
       a[j] = atan2(B2.row(i).dot(DIRs.row(j)),B1.row(i).dot(DIRs.row(j))) - a_combed;
       //make it positive by adding some multiple of 2pi
-      a[j] += ceil (std::max(0., -a[j]) / (M_PI*2.)) * (M_PI*2.);
+      a[j] += std::ceil (std::max(0., -a[j]) / (M_PI*2.)) * (M_PI*2.);
       //take modulo 2pi
       a[j] = fmod(a[j], (M_PI*2.));
     }

include/igl/compute_frame_field_bisectors.cpp

@@ -16,14 +16,14 @@
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::compute_frame_field_bisectors(
-                                                   const Eigen::PlainObjectBase<DerivedV>& V,
-                                                   const Eigen::PlainObjectBase<DerivedF>& F,
-                                                   const Eigen::PlainObjectBase<DerivedV>& B1,
-                                                   const Eigen::PlainObjectBase<DerivedV>& B2,
-                                                   const Eigen::PlainObjectBase<DerivedV>& PD1,
-                                                   const Eigen::PlainObjectBase<DerivedV>& PD2,
-                                                   Eigen::PlainObjectBase<DerivedV>& BIS1,
-                                                   Eigen::PlainObjectBase<DerivedV>& BIS2)
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedV>& B1,
+  const Eigen::PlainObjectBase<DerivedV>& B2,
+  const Eigen::PlainObjectBase<DerivedV>& PD1,
+  const Eigen::PlainObjectBase<DerivedV>& PD2,
+  Eigen::PlainObjectBase<DerivedV>& BIS1,
+  Eigen::PlainObjectBase<DerivedV>& BIS2)
 {
   BIS1.resize(PD1.rows(),3);
   BIS2.resize(PD1.rows(),3);
@@ -34,24 +34,24 @@ IGL_INLINE void igl::compute_frame_field_bisectors(
     // Convert to angle
     double a1 = atan2(B2.row(i).dot(PD1.row(i)),B1.row(i).dot(PD1.row(i)));
     //make it positive by adding some multiple of 2pi
-    a1 += ceil (std::max(0., -a1) / (M_PI*2.)) * (M_PI*2.);
+    a1 += std::ceil (std::max(0., -a1) / (M_PI*2.)) * (M_PI*2.);
     //take modulo 2pi
     a1 = fmod(a1, (M_PI*2.));
     double a2 = atan2(B2.row(i).dot(PD2.row(i)),B1.row(i).dot(PD2.row(i)));
     //make it positive by adding some multiple of 2pi
-    a2 += ceil (std::max(0., -a2) / (M_PI*2.)) * (M_PI*2.);
+    a2 += std::ceil (std::max(0., -a2) / (M_PI*2.)) * (M_PI*2.);
     //take modulo 2pi
     a2 = fmod(a2, (M_PI*2.));
 
     double b1 = (a1+a2)/2.0;
     //make it positive by adding some multiple of 2pi
-    b1 += ceil (std::max(0., -b1) / (M_PI*2.)) * (M_PI*2.);
+    b1 += std::ceil (std::max(0., -b1) / (M_PI*2.)) * (M_PI*2.);
     //take modulo 2pi
     b1 = fmod(b1, (M_PI*2.));
 
     double b2 = b1+(M_PI/2.);
     //make it positive by adding some multiple of 2pi
-    b2 += ceil (std::max(0., -b2) / (M_PI*2.)) * (M_PI*2.);
+    b2 += std::ceil (std::max(0., -b2) / (M_PI*2.)) * (M_PI*2.);
     //take modulo 2pi
     b2 = fmod(b2, (M_PI*2.));
 

include/igl/copyleft/cgal/mesh_to_cgal_triangle_list.cpp

@@ -55,6 +55,10 @@ IGL_INLINE void igl::copyleft::cgal::mesh_to_cgal_triangle_list(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::copyleft::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);
+// generated by autoexplicit.sh
 template void igl::copyleft::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
 // generated by autoexplicit.sh
 template void igl::copyleft::cgal::mesh_to_cgal_triangle_list<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);

include/igl/copyleft/cgal/order_facets_around_edge.cpp

@@ -1,4 +1,5 @@
 #include "order_facets_around_edge.h"
+#include <Eigen/Geometry>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 #include <stdexcept>

include/igl/copyleft/cgal/order_facets_around_edges.cpp

@@ -8,6 +8,7 @@
 #include "order_facets_around_edges.h"
 #include "order_facets_around_edge.h"
 #include "../../sort_angles.h"
+#include <Eigen/Geometry>
 #include <type_traits>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
@@ -320,7 +321,12 @@ IGL_INLINE void igl::copyleft::cgal::order_facets_around_edges(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, int, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
+// generated by autoexplicit.sh
+template std::enable_if<!(std::is_same<Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, CGAL::Lazy_exact_nt<CGAL::Gmpq> >::value), void>::type igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, int, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
 template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(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, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
+template void igl::copyleft::cgal::order_facets_around_edges<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, long, long, bool>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > >&);
 #endif

include/igl/copyleft/cgal/outer_element.cpp

@@ -214,4 +214,6 @@ template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<
 template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 1, -1, -1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 1, -1, -1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/outer_facet.cpp

@@ -149,6 +149,8 @@ IGL_INLINE void igl::copyleft::cgal::outer_facet(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, unsigned long&, bool&);
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
@@ -160,4 +162,6 @@ template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, 3, 0, -
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, unsigned long&, bool&);
 template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, int&, bool&);
+template void igl::copyleft::cgal::outer_facet<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int&, bool&);
 #endif

include/igl/copyleft/cgal/outer_hull.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "outer_hull.h"
 #include "extract_cells.h"
@@ -35,12 +35,12 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull(
   // Exact types
   typedef CGAL::Epeck Kernel;
   typedef Kernel::FT ExactScalar;
-  typedef 
+  typedef
     Eigen::Matrix<
     ExactScalar,
     Eigen::Dynamic,
     Eigen::Dynamic,
-    DerivedHV::IsRowMajor> 
+    DerivedHV::IsRowMajor>
       MatrixXES;
   // Remesh self-intersections
   MatrixXES Vr;
@@ -535,4 +535,11 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull_legacy(
 // Explicit template specialization
 template void igl::copyleft::cgal::outer_hull<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::outer_hull_legacy< Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::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> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > &, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > &);
+template void igl::copyleft::cgal::outer_hull_legacy< Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::copyleft::cgal::outer_hull_legacy<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::outer_hull_legacy<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#include "../../barycenter.cpp"
+template void igl::barycenter<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
+#include "../../remove_unreferenced.cpp"
+template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -115,5 +115,9 @@ IGL_INLINE size_t igl::copyleft::cgal::peel_outer_hull_layers(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+template unsigned long igl::copyleft::cgal::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template size_t igl::copyleft::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/points_inside_component.cpp

@@ -344,4 +344,5 @@ template void igl::copyleft::cgal::points_inside_component< Eigen::Matrix<double
 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::points_inside_component<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -513,6 +513,10 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -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&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/copyleft/cgal/remesh_self_intersections.cpp

@@ -96,4 +96,5 @@ template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<CGAL:
 template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/progressive_hulls.cpp

@@ -18,6 +18,7 @@ IGL_INLINE bool igl::copyleft::progressive_hulls(
   Eigen::VectorXi & J)
 {
   int m = F.rows();
+  Eigen::VectorXi I;
   return decimate(
     V,
     F,
@@ -25,5 +26,6 @@ IGL_INLINE bool igl::copyleft::progressive_hulls(
     max_faces_stopping_condition(m,max_m),
     U,
     G,
-    J);
+    J,
+    I);
 }

include/igl/cotmatrix.h

@@ -43,8 +43,6 @@ namespace igl
   // therefore in general negative and the matrix is **negative** semi-definite
   // (immediately, -L is **positive** semi-definite)
   //
-  // Known bugs: off by 1e-16 on regular grid. I think its a problem of
-  // arithmetic order in cotmatrix_entries.h: C(i,e) = (arithmetic)/dblA/4
   template <typename DerivedV, typename DerivedF, typename Scalar>
   IGL_INLINE void cotmatrix(
     const Eigen::PlainObjectBase<DerivedV> & V, 

include/igl/cotmatrix_entries.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cotmatrix_entries.h"
 #include "doublearea.h"
+#include "squared_edge_lengths.h"
 #include "edge_lengths.h"
 #include "face_areas.h"
 #include "volume.h"
@@ -34,11 +35,13 @@ IGL_INLINE void igl::cotmatrix_entries(
     case 3:
     {
       // Triangles
-      //Matrix<typename DerivedC::Scalar,Dynamic,3> l;
-      //edge_lengths(V,F,l);
-      // edge lengths numbered same as opposite vertices
+      //Compute Squared Edge lenghts 
+      Matrix<typename DerivedC::Scalar,Dynamic,3> l2;
+      igl::squared_edge_lengths(V,F,l2);
+      //Compute Edge lenghts 
       Matrix<typename DerivedC::Scalar,Dynamic,3> l;
-      igl::edge_lengths(V,F,l);
+      l = l2.array().sqrt();
+      
       // double area
       Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
       doublearea(l,dblA);
@@ -47,9 +50,9 @@ IGL_INLINE void igl::cotmatrix_entries(
       C.resize(m,3);
       for(int i = 0;i<m;i++)
       {
-        C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
-        C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
-        C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
+        C(i,0) = (l2(i,1) + l2(i,2) - l2(i,0))/dblA(i)/4.0;
+        C(i,1) = (l2(i,2) + l2(i,0) - l2(i,1))/dblA(i)/4.0;
+        C(i,2) = (l2(i,0) + l2(i,1) - l2(i,2))/dblA(i)/4.0;
       }
       break;
     }

include/igl/decimate.cpp

@@ -10,6 +10,7 @@
 #include "edge_flaps.h"
 #include "remove_unreferenced.h"
 #include "slice_mask.h"
+#include "slice.h"
 #include "connect_boundary_to_infinity.h"
 #include <iostream>
 
@@ -19,8 +20,10 @@ IGL_INLINE bool igl::decimate(
   const size_t max_m,
   Eigen::MatrixXd & U,
   Eigen::MatrixXi & G,
-  Eigen::VectorXi & J)
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & I)
 {
+
   // Original number of faces
   const int orig_m = F.rows();
   // Tracking number of faces
@@ -73,15 +76,29 @@ IGL_INLINE bool igl::decimate(
     max_non_infinite_faces_stopping_condition,
     U,
     G,
-    J);
+    J,
+    I);
   const Eigen::Array<bool,Eigen::Dynamic,1> keep = (J.array()<orig_m);
   igl::slice_mask(Eigen::MatrixXi(G),keep,1,G);
   igl::slice_mask(Eigen::VectorXi(J),keep,1,J);
-  Eigen::VectorXi _1;
-  igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1);
+  Eigen::VectorXi _1,I2;
+  igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1,I2);
+  igl::slice(Eigen::VectorXi(I),I2,1,I);
   return ret;
 }
 
+IGL_INLINE bool igl::decimate(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const size_t max_m,
+  Eigen::MatrixXd & U,
+  Eigen::MatrixXi & G,
+  Eigen::VectorXi & J)
+{
+  Eigen::VectorXi I;
+  return igl::decimate(V,F,max_m,U,G,J,I);
+}
+
 IGL_INLINE bool igl::decimate(
   const Eigen::MatrixXd & OV,
   const Eigen::MatrixXi & OF,
@@ -112,7 +129,9 @@ IGL_INLINE bool igl::decimate(
       const int)> & stopping_condition,
   Eigen::MatrixXd & U,
   Eigen::MatrixXi & G,
-  Eigen::VectorXi & J)
+  Eigen::VectorXi & J,
+  Eigen::VectorXi & I
+  )
 {
   using namespace Eigen;
   using namespace std;
@@ -189,6 +208,6 @@ IGL_INLINE bool igl::decimate(
   F2.conservativeResize(m,F2.cols());
   J.conservativeResize(m);
   VectorXi _1;
-  remove_unreferenced(V,F2,U,G,_1);
+  remove_unreferenced(V,F2,U,G,_1,I);
   return clean_finish;
 }

include/igl/decimate.h

@@ -25,6 +25,24 @@ namespace igl
   //   U  #U by dim list of output vertex posistions (can be same ref as V)
   //   G  #G by 3 list of output face indices into U (can be same ref as G)
   //   J  #G list of indices into F of birth face
+  //   I  #U list of indices into V of birth vertices
+  // Returns true if m was reached (otherwise #G > m)
+  IGL_INLINE bool decimate(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const size_t max_m,
+    Eigen::MatrixXd & U,
+    Eigen::MatrixXi & G,
+    Eigen::VectorXi & J,
+    Eigen::VectorXi & I);
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   F  #F by 3 list of face indices into V.
+  //   max_m  desired number of output faces
+  // Outputs:
+  //   U  #U by dim list of output vertex posistions (can be same ref as V)
+  //   G  #G by 3 list of output face indices into U (can be same ref as G)
+  //   J  #G list of indices into F of birth face
   // Returns true if m was reached (otherwise #G > m)
   IGL_INLINE bool decimate(
     const Eigen::MatrixXd & V,
@@ -79,7 +97,8 @@ namespace igl
       )> & stopping_condition,
     Eigen::MatrixXd & U,
     Eigen::MatrixXi & G,
-    Eigen::VectorXi & J);
+    Eigen::VectorXi & J,
+    Eigen::VectorXi & I);
 
 }
 

include/igl/doublearea.cpp

@@ -146,10 +146,10 @@ IGL_INLINE void igl::doublearea(
   //
   // "Miscalculating Area and Angles of a Needle-like Triangle"
   // https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf
-  sort(ul,2,false,l,_);
+  igl::sort(ul,2,false,l,_);
   // semiperimeters
-  Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-  assert((Index)s.rows() == m);
+  //Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
+  //assert((Index)s.rows() == m);
   // resize output
   dblA.resize(l.rows(),1);
   parallel_for(

include/igl/edge_flaps.h

@@ -26,6 +26,12 @@ namespace igl
   //   EI  #E by 2 list of edge flap corners (see above).
   //
   // TODO: This seems to be a duplicate of edge_topology.h
+  // igl::edge_topology(V,F,etEV,etFE,etEF);
+  // igl::edge_flaps(F,efE,efEMAP,efEF,efEI);
+  // [~,I] = sort(efE,2)
+  // all( efE(sub2ind(size(efE),repmat(1:size(efE,1),2,1)',I)) == etEV )
+  // all( efEF(sub2ind(size(efE),repmat(1:size(efE,1),2,1)',I)) == etEF )
+  // all(efEMAP(sub2ind(size(F),repmat(1:size(F,1),3,1)',repmat([1 2 3],size(F,1),1))) == etFE(:,[2 3 1]))
   IGL_INLINE void edge_flaps(
     const Eigen::MatrixXi & F,
     const Eigen::MatrixXi & E,

include/igl/edge_lengths.cpp

@@ -6,69 +6,18 @@
 // 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_lengths.h"
-#include "parallel_for.h"
-#include <iostream>
+#include "squared_edge_lengths.h"
 
 template <typename DerivedV, typename DerivedF, typename DerivedL>
 IGL_INLINE void igl::edge_lengths(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedL>& L)
-{
-  using namespace std;
-  const int m = F.rows();
-  switch(F.cols())
   {
-    case 2:
-    {
-      L.resize(F.rows(),1);
-      for(int i = 0;i<F.rows();i++)
-      {
-        L(i,0) = (V.row(F(i,1))-V.row(F(i,0))).norm();
-      }
-      break;
-    }
-    case 3:
-    {
-      L.resize(m,3);
-      // loop over faces
-      parallel_for(
-        m,
-        [&V,&F,&L](const int i)
-        {
-          L(i,0) = (V.row(F(i,1))-V.row(F(i,2))).norm();
-          L(i,1) = (V.row(F(i,2))-V.row(F(i,0))).norm();
-          L(i,2) = (V.row(F(i,0))-V.row(F(i,1))).norm();
-        },
-        1000);
-      break;
-    }
-    case 4:
-    {
-      L.resize(m,6);
-      // loop over faces
-      parallel_for(
-        m,
-        [&V,&F,&L](const int i)
-        {
-          L(i,0) = (V.row(F(i,3))-V.row(F(i,0))).norm();
-          L(i,1) = (V.row(F(i,3))-V.row(F(i,1))).norm();
-          L(i,2) = (V.row(F(i,3))-V.row(F(i,2))).norm();
-          L(i,3) = (V.row(F(i,1))-V.row(F(i,2))).norm();
-          L(i,4) = (V.row(F(i,2))-V.row(F(i,0))).norm();
-          L(i,5) = (V.row(F(i,0))-V.row(F(i,1))).norm();
-        },
-        1000);
-      break;
-    }
-    default:
-    {
-      cerr<< "edge_lengths.h: Error: Simplex size ("<<F.cols()<<
-        ") not supported"<<endl;
-      assert(false);
-    }
+      igl::squared_edge_lengths(V,F,L);
+      L=L.array().sqrt().eval();
   }
-}
+  
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization

include/igl/find.cpp

@@ -127,4 +127,5 @@ template void igl::find<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::M
 // generated by autoexplicit.sh
 template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::find<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::find<double, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/flip_avoiding_line_search.cpp

@@ -6,252 +6,296 @@
 // 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 "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;
+namespace igl
+{
+  namespace flip_avoiding
+  {
+    //---------------------------------------------------------------------------
+    // 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
+    IGL_INLINE int SolveP3(std::vector<double>& x,double a,double b,double c)
+    { // solve cubic equation x^3 + a*x^2 + b*x + c
+      using namespace std;
+      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+(2*M_PI))/3)-a;
+          x[2]=q*cos((t-(2*M_PI))/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 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;
+    IGL_INLINE double get_smallest_pos_quad_zero(double a,double b, double c)
+    {
+      using namespace std;
+      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;
+      }
     }
-  } 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;
+    IGL_INLINE double get_min_pos_root_2D(const Eigen::MatrixXd& uv,
+                                          const Eigen::MatrixXi& F,
+                                          Eigen::MatrixXd& d,
+                                          int f)
+    {
+      using namespace std;
+    /*
+          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)
+      const double& U11 = uv(v1,0);
+      const double& U12 = uv(v1,1);
+      const double& U21 = uv(v2,0);
+      const double& U22 = uv(v2,1);
+      const double& U31 = uv(v3,0);
+      const double& U32 = uv(v3,1);
+
+      const double& V11 = d(v1,0);
+      const double& V12 = d(v1,1);
+      const double& V21 = d(v2,0);
+      const double& V22 = d(v2,1);
+      const double& V31 = d(v3,0);
+      const double& 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);
     }
-    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;
+    IGL_INLINE double get_min_pos_root_3D(const Eigen::MatrixXd& uv,
+                                          const Eigen::MatrixXi& F,
+                                          Eigen::MatrixXd& direc,
+                                          int f)
+    {
+      using namespace std;
+      /*
+          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);
+      const double& a_x = uv(v1,0);
+      const double& a_y = uv(v1,1);
+      const double& a_z = uv(v1,2);
+      const double& b_x = uv(v2,0);
+      const double& b_y = uv(v2,1);
+      const double& b_z = uv(v2,2);
+      const double& c_x = uv(v3,0);
+      const double& c_y = uv(v3,1);
+      const double& c_z = uv(v3,2);
+      const double& d_x = uv(v4,0);
+      const double& d_y = uv(v4,1);
+      const double& d_z = uv(v4,2);
+
+      const double& a_dx = direc(v1,0);
+      const double& a_dy = direc(v1,1);
+      const double& a_dz = direc(v1,2);
+      const double& b_dx = direc(v2,0);
+      const double& b_dy = direc(v2,1);
+      const double& b_dz = direc(v2,2);
+      const double& c_dx = direc(v3,0);
+      const double& c_dy = direc(v3,1);
+      const double& c_dz = direc(v3,2);
+      const double& d_dx = direc(v4,0);
+      const double& d_dy = direc(v4,1);
+      const double& 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;
+        }
+      }
+    }
 
-    // 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);
+    IGL_INLINE double compute_max_step_from_singularities(const Eigen::MatrixXd& uv,
+                                                          const Eigen::MatrixXi& F,
+                                                          Eigen::MatrixXd& d)
+    {
+      using namespace std;
+      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);
+      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;
     }
-    return max_step;
- }
+  }
+}
 
 IGL_INLINE double igl::flip_avoiding_line_search(
   const Eigen::MatrixXi F,
@@ -260,12 +304,13 @@ IGL_INLINE double igl::flip_avoiding_line_search(
   std::function<double(Eigen::MatrixXd&)> energy,
   double cur_energy)
 {
-    Eigen::MatrixXd d = dst_v - cur_v;
+  using namespace std;
+  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);
+  double min_step_to_singularity = igl::flip_avoiding::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);
+  return igl::line_search(cur_v,d,max_step_size, energy, cur_energy);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/flip_avoiding_line_search.h

@@ -7,13 +7,13 @@
 // 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 "igl_inline.h"
 
 #include <Eigen/Dense>
 
 namespace igl
 {
-  // A bisection line search for a mesh based energy that avoids triangle flips as suggested in 
+  // 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).
@@ -27,7 +27,7 @@ namespace igl
   //   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. 
+  //   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

include/igl/harmonic.cpp

@@ -35,16 +35,7 @@ IGL_INLINE bool igl::harmonic(
   typedef typename DerivedV::Scalar Scalar;
   SparseMatrix<Scalar> L,M;
   cotmatrix(V,F,L);
-  switch(F.cols())
-  {
-    case 3:
-      massmatrix(V,F,MASSMATRIX_TYPE_VORONOI,M);
-      break;
-    case 4:
-    default:
-      massmatrix(V,F,MASSMATRIX_TYPE_BARYCENTRIC,M);
-    break;
-  }
+  massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
   return harmonic(L,M,b,bc,k,W);
 }
 
@@ -96,14 +87,9 @@ IGL_INLINE bool igl::harmonic(
   assert(n == M.rows() && "M must be square");
   assert(igl::isdiag(M) && "Mass matrix should be diagonal");
 
-  Eigen::SparseMatrix<DerivedM> Mi;
-  invert_diag(M,Mi);
   Eigen::SparseMatrix<DerivedL> Q;
-  Q = -L;
-  for(int p = 1;p<k;p++)
-  {
-    Q = (Q*Mi*-L).eval();
-  }
+  igl::harmonic(L,M,k,Q);
+
   typedef DerivedL Scalar;
   min_quad_with_fixed_data<Scalar> data;
   min_quad_with_fixed_precompute(Q,b,Eigen::SparseMatrix<Scalar>(),true,data);
@@ -123,6 +109,43 @@ IGL_INLINE bool igl::harmonic(
   return true;
 }
 
+template <
+  typename DerivedL,
+  typename DerivedM,
+  typename DerivedQ>
+IGL_INLINE void igl::harmonic(
+  const Eigen::SparseMatrix<DerivedL> & L,
+  const Eigen::SparseMatrix<DerivedM> & M,
+  const int k,
+  Eigen::SparseMatrix<DerivedQ> & Q)
+{
+  assert(L.rows() == L.cols()&&"L should be square");
+  assert(M.rows() == M.cols()&&"M should be square");
+  assert(L.rows() == M.rows()&&"L should match M's dimensions");
+  Eigen::SparseMatrix<DerivedM> Mi;
+  invert_diag(M,Mi);
+  Q = -L;
+  for(int p = 1;p<k;p++)
+  {
+    Q = (Q*Mi*-L).eval();
+  }
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedQ>
+IGL_INLINE void igl::harmonic(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const int k,
+  Eigen::SparseMatrix<DerivedQ> & Q)
+{
+  Eigen::SparseMatrix<DerivedQ> L,M;
+  cotmatrix(V,F,L);
+  massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
+  return harmonic(L,M,k,Q);
+}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
@@ -133,4 +156,5 @@ template bool igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Mat
 template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template bool igl::harmonic<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::harmonic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/harmonic.h

@@ -81,6 +81,39 @@ namespace igl
     const Eigen::PlainObjectBase<Derivedbc> & bc,
     const int k,
     Eigen::PlainObjectBase<DerivedW> & W);
+  // Build the discrete k-harmonic operator (computing integrated quantities).
+  // That is, if the k-harmonic PDE is Q x = 0, then this minimizes x' Q x
+  //
+  // Inputs:
+  //   L  #V by #V discrete (integrated) Laplacian  
+  //   M  #V by #V mass matrix
+  //   k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
+  // Outputs:
+  //   Q  #V by #V discrete (integrated) k-Laplacian  
+  template <
+    typename DerivedL,
+    typename DerivedM,
+    typename DerivedQ>
+  IGL_INLINE void harmonic(
+    const Eigen::SparseMatrix<DerivedL> & L,
+    const Eigen::SparseMatrix<DerivedM> & M,
+    const int k,
+    Eigen::SparseMatrix<DerivedQ> & Q);
+  // Inputs:
+  //   V  #V by dim vertex positions
+  //   F  #F by simplex-size list of element indices
+  //   k  power of harmonic operation (1: harmonic, 2: biharmonic, etc)
+  // Outputs:
+  //   Q  #V by #V discrete (integrated) k-Laplacian  
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedQ>
+  IGL_INLINE void harmonic(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const int k,
+    Eigen::SparseMatrix<DerivedQ> & Q);
 };
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/internal_angles.cpp

@@ -7,7 +7,7 @@
 // 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 "internal_angles.h"
-#include "edge_lengths.h"
+#include "squared_edge_lengths.h"
 #include "parallel_for.h"
 #include "get_seconds.h"
 
@@ -26,11 +26,11 @@ IGL_INLINE void igl::internal_angles(
     Matrix<
       Scalar,
       DerivedF::RowsAtCompileTime,
-      DerivedF::ColsAtCompileTime> L;
-    edge_lengths(V,F,L);
+      DerivedF::ColsAtCompileTime> L_sq;
+      igl::squared_edge_lengths(V,F,L_sq);
 
     assert(F.cols() == 3 && "F should contain triangles");
-    internal_angles(L,K);
+    igl::internal_angles_using_squared_edge_lengths(L_sq,K);
   }else
   {
     assert(V.cols() == 3 && "If F contains non-triangle facets, V must be 3D");
@@ -63,10 +63,37 @@ IGL_INLINE void igl::internal_angles(
 }
 
 template <typename DerivedL, typename DerivedK>
-IGL_INLINE void igl::internal_angles(
+IGL_INLINE void igl::internal_angles_using_squared_edge_lengths(
+  const Eigen::PlainObjectBase<DerivedL>& L_sq,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  typedef typename DerivedL::Index Index;
+  assert(L_sq.cols() == 3 && "Edge-lengths should come from triangles");
+  const Index m = L_sq.rows();
+  K.resize(m,3);
+  parallel_for(
+    m,
+    [&L_sq,&K](const Index f)
+    {
+      for(size_t d = 0;d<3;d++)
+      {
+        const auto & s1 = L_sq(f,d);
+        const auto & s2 = L_sq(f,(d+1)%3);
+        const auto & s3 = L_sq(f,(d+2)%3);
+        K(f,d) = acos((s3 + s2 - s1)/(2.*sqrt(s3*s2)));
+      }
+    },
+    1000l);
+}
+
+template <typename DerivedL, typename DerivedK>
+IGL_INLINE void igl::internal_angles_using_edge_lengths(
   const Eigen::PlainObjectBase<DerivedL>& L,
   Eigen::PlainObjectBase<DerivedK> & K)
 {
+  // Note:
+  //   Usage of internal_angles_using_squared_edge_lengths() is preferred to internal_angles_using_squared_edge_lengths()
+  //   This function is deprecated and probably will be removed in future versions
   typedef typename DerivedL::Index Index;
   assert(L.cols() == 3 && "Edge-lengths should come from triangles");
   const Index m = L.rows();

include/igl/internal_angles.h

@@ -28,12 +28,30 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedK> & K);
   // Inputs:
+  //   L_sq  #F by 3 list of squared edge lengths
+  // Output:
+  //   K  #F by poly-size eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //
+  // Note:
+  //   Usage of internal_angles_using_squared_edge_lengths is preferred to internal_angles_using_squared_edge_lengths
+  template <typename DerivedL, typename DerivedK>
+  IGL_INLINE void internal_angles_using_squared_edge_lengths(
+    const Eigen::PlainObjectBase<DerivedL>& L_sq,
+    Eigen::PlainObjectBase<DerivedK> & K);
+  // Inputs:
   //   L  #F by 3 list of edge lengths
+  // Output:
+  //   K  #F by poly-size eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //
+  // Note:
+  //   Usage of internal_angles_using_squared_edge_lengths is preferred to internal_angles_using_squared_edge_lengths
+  //   This function is deprecated and probably will be removed in future versions
   template <typename DerivedL, typename DerivedK>
-  IGL_INLINE void internal_angles(
+  IGL_INLINE void internal_angles_using_edge_lengths(
     const Eigen::PlainObjectBase<DerivedL>& L,
-    Eigen::PlainObjectBase<DerivedK> & K);
-}
+    Eigen::PlainObjectBase<DerivedK> & K);}
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "internal_angles.cpp"

include/igl/is_vertex_manifold.cpp

@@ -98,4 +98,5 @@ IGL_INLINE bool igl::is_vertex_manifold(
 #ifdef IGL_STATIC_LIBRARY
 template bool igl::is_vertex_manifold<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 bool igl::is_vertex_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> >&);
+template bool igl::is_vertex_manifold<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> >&);
 #endif

include/igl/line_field_missmatch.cpp

@@ -67,7 +67,7 @@ private:
         double angle_diff = atan2(dir1Rot.dot(PD2.row(f0)),dir1Rot.dot(PD1.row(f0)));
 
         double step=M_PI;
-        int i=(int)floor((angle_diff/step)+0.5);
+        int i=(int)std::floor((angle_diff/step)+0.5);
         assert((i>=-2)&&(i<=2));
         int k=0;
         if (i>=0)

include/igl/line_search.cpp

@@ -15,21 +15,28 @@ IGL_INLINE double igl::line_search(
   double cur_energy)
 {
   double old_energy;
-  if (cur_energy > 0) {
+  if (cur_energy > 0)
+  {
     old_energy = cur_energy;
-  } else {
+  }
+  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) {
+  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) {
+    if ( cur_e >= old_energy)
+    {
       step_size /= 2;
-    } else {
+    }
+    else
+    {
       x = new_x;
       new_energy = cur_e;
     }

include/igl/line_search.h

@@ -7,7 +7,7 @@
 // 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 "igl_inline.h"
 
 #include <Eigen/Dense>
 
@@ -21,7 +21,7 @@ namespace igl
   //   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. 
+  //   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

include/igl/list_to_matrix.cpp

@@ -119,6 +119,8 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Eigen::PlainObject
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 1, -1, -1> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&);
+// generated by autoexplicit.sh
 template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 // generated by autoexplicit.sh
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);

include/igl/loop.cpp

@@ -149,4 +149,4 @@ namespace igl
         }
     }
     
-}
+}

include/igl/mosek/bbw.cpp

@@ -0,0 +1,88 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "bbw.h"
+#include "mosek_quadprog.h"
+#include "../harmonic.h"
+#include "../slice_into.h"
+#include <Eigen/Sparse>
+#include <iostream>
+#include <cstdio>
+
+
+template <
+  typename DerivedV,
+  typename DerivedEle,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedW>
+IGL_INLINE bool igl::mosek::bbw(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedEle> & Ele,
+  const Eigen::PlainObjectBase<Derivedb> & b,
+  const Eigen::PlainObjectBase<Derivedbc> & bc,
+  igl::BBWData & data,
+  igl::mosek::MosekData & mosek_data,
+  Eigen::PlainObjectBase<DerivedW> & W
+  )
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(!data.partition_unity && "partition_unity not implemented yet");
+  // number of domain vertices
+  int n = V.rows();
+  // number of handles
+  int m = bc.cols();
+  // Build biharmonic operator
+  Eigen::SparseMatrix<typename DerivedV::Scalar> Q;
+  harmonic(V,Ele,2,Q);
+  W.derived().resize(n,m);
+  // No linear terms
+  VectorXd c = VectorXd::Zero(n);
+  // No linear constraints
+  SparseMatrix<typename DerivedW::Scalar> A(0,n);
+  VectorXd uc(0,1),lc(0,1);
+  // Upper and lower box constraints (Constant bounds)
+  VectorXd ux = VectorXd::Ones(n);
+  VectorXd lx = VectorXd::Zero(n);
+  // Loop over handles
+  for(int i = 0;i<m;i++)
+  {
+    if(data.verbosity >= 1)
+    {
+      cout<<"BBW: Computing weight for handle "<<i+1<<" out of "<<m<<
+        "."<<endl;
+    }
+    VectorXd bci = bc.col(i);
+    VectorXd Wi;
+    // impose boundary conditions via bounds
+    slice_into(bci,b,ux);
+    slice_into(bci,b,lx);
+    bool r = mosek_quadprog(Q,c,0,A,lc,uc,lx,ux,mosek_data,Wi);
+    if(!r)
+    {
+      return false;
+    }
+    W.col(i) = Wi;
+  }
+#ifndef NDEBUG
+    const double min_rowsum = W.rowwise().sum().array().abs().minCoeff();
+    if(min_rowsum < 0.1)
+    {
+      cerr<<"bbw.cpp: Warning, minimum row sum is very low. Consider more "
+        "active set iterations or enforcing partition of unity."<<endl;
+    }
+#endif
+
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::mosek::bbw<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::BBWData&, igl::mosek::MosekData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif
+

include/igl/mosek/bbw.h

@@ -0,0 +1,60 @@
+// 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_MOSEK_BBW_H
+#define IGL_MOSEK_BBW_H
+#include "../igl_inline.h"
+#include "mosek_quadprog.h"
+#include "../bbw.h"
+#include <Eigen/Dense>
+
+namespace igl
+{
+  namespace mosek
+  {
+    // 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
+    //   mosek_data  object containing mosek options
+    // 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
+    template <
+      typename DerivedV,
+      typename DerivedEle,
+      typename Derivedb,
+      typename Derivedbc,
+      typename DerivedW>
+    IGL_INLINE bool bbw(
+      const Eigen::PlainObjectBase<DerivedV> & V,
+      const Eigen::PlainObjectBase<DerivedEle> & Ele,
+      const Eigen::PlainObjectBase<Derivedb> & b,
+      const Eigen::PlainObjectBase<Derivedbc> & bc,
+      igl::BBWData & data,
+      igl::mosek::MosekData & mosek_data,
+      Eigen::PlainObjectBase<DerivedW> & W);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bbw.cpp"
+#endif
+
+#endif

include/igl/procrustes.cpp

@@ -42,7 +42,7 @@ IGL_INLINE void igl::procrustes(
      double scaleY = YC.norm() / YC.rows();
      scale = scaleY/scaleX;
      XC *= scale;
-     assert (abs(XC.norm() / XC.rows() - scaleY) < 1e-8);
+     assert (std::abs(XC.norm() / XC.rows() - scaleY) < 1e-8);
   }
 
   // Rotation

include/igl/readMESH.cpp

@@ -465,6 +465,8 @@ IGL_INLINE bool igl::readMESH(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template bool igl::readMESH<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template bool igl::readMESH<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh
 template bool igl::readMESH<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);

include/igl/readOFF.cpp

@@ -14,7 +14,8 @@ IGL_INLINE bool igl::readOFF(
   const std::string off_file_name,
   std::vector<std::vector<Scalar > > & V,
   std::vector<std::vector<Index > > & F,
-  std::vector<std::vector<Scalar > > & N)
+  std::vector<std::vector<Scalar > > & N,
+  std::vector<std::vector<Scalar > > & C)
 {
   using namespace std;
   FILE * off_file = fopen(off_file_name.c_str(),"r");
@@ -26,20 +27,25 @@ IGL_INLINE bool igl::readOFF(
   V.clear();
   F.clear();
   N.clear();
+  C.clear();
+
   // First line is always OFF
   char header[1000];
   const std::string OFF("OFF");
   const std::string NOFF("NOFF");
+  const std::string COFF("COFF");
   if(fscanf(off_file,"%s\n",header)!=1
      || !(
        string(header).compare(0, OFF.length(), OFF)==0 ||
+       string(header).compare(0, COFF.length(), COFF)==0 ||
        string(header).compare(0,NOFF.length(),NOFF)==0))
   {
-    printf("Error: %s's first line should be OFF or NOFF not %s...",off_file_name.c_str(),header);
+    printf("Error: readOFF() %s's first line should be OFF or NOFF or COFF, not %s...",off_file_name.c_str(),header);
     fclose(off_file);
     return false;
   }
   bool has_normals = string(header).compare(0,NOFF.length(),NOFF)==0;
+  bool has_vertexColors = string(header).compare(0,COFF.length(),COFF)==0;
   // Second line is #vertices #faces #edges
   int number_of_vertices;
   int number_of_faces;
@@ -56,6 +62,8 @@ IGL_INLINE bool igl::readOFF(
   V.resize(number_of_vertices);
   if (has_normals)
     N.resize(number_of_vertices);
+  if (has_vertexColors)
+    C.resize(number_of_vertices);
   F.resize(number_of_faces);
   //printf("%s %d %d %d\n",(has_normals ? "NOFF" : "OFF"),number_of_vertices,number_of_faces,number_of_edges);
   // Read vertices
@@ -81,6 +89,14 @@ IGL_INLINE bool igl::readOFF(
         normal[2] = nz;
         N[i] = normal;
       }
+
+      if (has_vertexColors)
+      {
+        C[i].resize(3);
+        C[i][0] = nx / 255.0;
+        C[i][1] = ny / 255.0;
+        C[i][2] = nz / 255.0;
+      }
       i++;
     }else if(
         fscanf(off_file,"%[#]",&tic_tac_toe)==1)
@@ -146,10 +162,11 @@ IGL_INLINE bool igl::readOFF(
   std::vector<std::vector<double> > vV;
   std::vector<std::vector<double> > vN;
   std::vector<std::vector<int> > vF;
-  bool success = igl::readOFF(str,vV,vF,vN);
+  std::vector<std::vector<double> > vC;
+  bool success = igl::readOFF(str,vV,vF,vN,vC);
   if(!success)
   {
-    // readOFF(str,vV,vF) should have already printed an error
+    // readOFF(str,vV,vF,vN,vC) should have already printed an error
     // message to stderr
     return false;
   }
@@ -179,10 +196,11 @@ IGL_INLINE bool igl::readOFF(
   std::vector<std::vector<double> > vV;
   std::vector<std::vector<double> > vN;
   std::vector<std::vector<int> > vF;
-  bool success = igl::readOFF(str,vV,vF,vN);
+  std::vector<std::vector<double> > vC;
+  bool success = igl::readOFF(str,vV,vF,vN,vC);
   if(!success)
   {
-    // readOFF(str,vV,vF) should have already printed an error
+    // readOFF(str,vV,vF,vC) should have already printed an error
     // message to stderr
     return false;
   }
@@ -208,6 +226,18 @@ IGL_INLINE bool igl::readOFF(
       return false;
     }
   }
+
+  //Warning: RGB colors will be returned in the N matrix
+  if (vC.size())
+  {
+    bool C_rect = igl::list_to_matrix(vC,N);
+    if(!C_rect)
+    {
+      // igl::list_to_matrix(vC,N) already printed error message to std err
+      return false;
+    }
+  }
+
   return true;
 }
 #endif

include/igl/readOFF.h

@@ -20,7 +20,7 @@
 namespace igl 
 {
   
-  // Read a mesh from an ascii obj file, filling in vertex positions, normals
+  // Read a mesh from an ascii OFF file, filling in vertex positions, normals
   // and texture coordinates. Mesh may have faces of any number of degree
   //
   // Templates:
@@ -32,17 +32,16 @@ namespace igl
   // Outputs:
   //   V  double matrix of vertex positions  #V by 3
   //   F  #F list of face indices into vertex positions
-  //   TC  double matrix of texture coordinats #TC by 2
-  //   FTC  #F list of face indices into vertex texture coordinates
-  //   N  double matrix of corner normals #N by 3
-  //   FN  #F list of face indices into vertex normals
+  //   N  list of vertex normals #V by 3
+  //   C  list of rgb color values per vertex #V by 3
   // Returns true on success, false on errors
   template <typename Scalar, typename Index>
   IGL_INLINE bool readOFF(
     const std::string off_file_name, 
     std::vector<std::vector<Scalar > > & V,
     std::vector<std::vector<Index > > & F,
-    std::vector<std::vector<Scalar > > & N);
+    std::vector<std::vector<Scalar > > & N,
+    std::vector<std::vector<Scalar > > & C);
   
   
 #ifndef IGL_NO_EIGEN

include/igl/readSTL.cpp

@@ -260,6 +260,8 @@ close_true:
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template bool igl::readSTL<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/read_triangle_mesh.cpp

@@ -35,7 +35,7 @@ IGL_INLINE bool igl::read_triangle_mesh(
   pathinfo(str,d,b,e,f);
   // Convert extension to lower case
   std::transform(e.begin(), e.end(), e.begin(), ::tolower);
-  vector<vector<Scalar> > TC, N;
+  vector<vector<Scalar> > TC, N, C;
   vector<vector<Index> > FTC, FN;
   if(e == "obj")
   {
@@ -49,7 +49,7 @@ IGL_INLINE bool igl::read_triangle_mesh(
     return success;
   }else if(e == "off")
   {
-    return readOFF(str,V,F,N);
+    return readOFF(str,V,F,N,C);
   }
   cerr<<"Error: "<<__FUNCTION__<<": "<<
     str<<" is not a recognized mesh file format."<<endl;
@@ -85,7 +85,7 @@ IGL_INLINE bool igl::read_triangle_mesh(
   pathinfo(filename,dir,base,ext,name);
   // Convert extension to lower case
   transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
-  vector<vector<double > > vV,vN,vTC;
+  vector<vector<double > > vV,vN,vTC,vC;
   vector<vector<int > > vF,vFTC,vFN;
   if(ext == "mesh")
   {
@@ -113,7 +113,7 @@ IGL_INLINE bool igl::read_triangle_mesh(
     }
   }else if(ext == "off")
   {
-    if(!readOFF(filename,vV,vF,vN))
+    if(!readOFF(filename,vV,vF,vN,vC))
     {
       return false;
     }
@@ -163,4 +163,5 @@ template bool igl::read_triangle_mesh<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Ei
 template bool igl::read_triangle_mesh<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
 template bool igl::read_triangle_mesh<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&);
 template bool igl::read_triangle_mesh<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >&, std::basic_string<char, std::char_traits<char>, std::allocator<char> >&);
+template bool igl::read_triangle_mesh<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/remove_unreferenced.cpp

@@ -99,8 +99,10 @@ IGL_INLINE void igl::remove_unreferenced(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, 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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > 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::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::remove_unreferenced<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/simplify_polyhedron.cpp

@@ -97,7 +97,8 @@ IGL_INLINE void igl::simplify_polyhedron(
     }
   };
   igl::per_face_normals(OV,OF,N);
+  Eigen::VectorXi I;
   igl::decimate(
-    OV,OF,perfect,igl::infinite_cost_stopping_condition(perfect),V,F,J);
+    OV,OF,perfect,igl::infinite_cost_stopping_condition(perfect),V,F,J,I);
 }
 

include/igl/slim.h

@@ -8,28 +8,21 @@
 #ifndef SLIM_H
 #define SLIM_H
 
+#include "igl_inline.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 {
-
+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 {
+  enum SLIM_ENERGY
+  {
     ARAP,
     LOG_ARAP,
     SYMMETRIC_DIRICHLET,
@@ -40,10 +33,10 @@ struct SLIMData {
   SLIM_ENERGY slim_energy;
 
   // Optional Input
-    // soft constraints
-    Eigen::VectorXi b;
-    Eigen::MatrixXd bc;
-    double soft_const_p;
+  // 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
@@ -60,30 +53,45 @@ struct SLIMData {
   int f_num;
   double proximal_p;
 
-  WeightedGlobalLocal* wGlobalLocal;
+  Eigen::VectorXd WGL_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;
 };
 
-  // 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);
+// 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);
+IGL_INLINE Eigen::MatrixXd slim_solve(SLIMData& data, int iter_num);
 
-}
+} // END NAMESPACE
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "slim.cpp"
 #endif
 
-
 #endif // SLIM_H

include/igl/squared_edge_lengths.cpp

@@ -0,0 +1,105 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "squared_edge_lengths.h"
+#include "parallel_for.h"
+#include <iostream>
+  
+template <typename DerivedV, typename DerivedF, typename DerivedL>
+IGL_INLINE void igl::squared_edge_lengths(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedL>& L)
+{
+  using namespace std;
+  const int m = F.rows();
+  switch(F.cols())
+  {
+    case 2:
+    {
+      L.resize(F.rows(),1);
+      for(int i = 0;i<F.rows();i++)
+      {
+        L(i,0) = (V.row(F(i,1))-V.row(F(i,0))).squaredNorm();
+      }
+      break;
+    }
+    case 3:
+    {
+      L.resize(m,3);
+      // loop over faces
+      parallel_for(
+        m,
+        [&V,&F,&L](const int i)
+        {
+          L(i,0) = (V.row(F(i,1))-V.row(F(i,2))).squaredNorm();
+          L(i,1) = (V.row(F(i,2))-V.row(F(i,0))).squaredNorm();
+          L(i,2) = (V.row(F(i,0))-V.row(F(i,1))).squaredNorm();
+        },
+        1000);
+      break;
+    }
+    case 4:
+    {
+      L.resize(m,6);
+      // loop over faces
+      parallel_for(
+        m,
+        [&V,&F,&L](const int i)
+        {
+          L(i,0) = (V.row(F(i,3))-V.row(F(i,0))).squaredNorm();
+          L(i,1) = (V.row(F(i,3))-V.row(F(i,1))).squaredNorm();
+          L(i,2) = (V.row(F(i,3))-V.row(F(i,2))).squaredNorm();
+          L(i,3) = (V.row(F(i,1))-V.row(F(i,2))).squaredNorm();
+          L(i,4) = (V.row(F(i,2))-V.row(F(i,0))).squaredNorm();
+          L(i,5) = (V.row(F(i,0))-V.row(F(i,1))).squaredNorm();
+        },
+        1000);
+      break;
+    }
+    default:
+    {
+      cerr<< "squared_edge_lengths.h: Error: Simplex size ("<<F.cols()<<
+        ") not supported"<<endl;
+      assert(false);
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(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, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -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::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(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, 6, 0, -1, 6> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<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::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<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> >&);
+#endif

include/igl/squared_edge_lengths.h

@@ -0,0 +1,49 @@
+// 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_SQUARED_EDGE_LENGTHS_H
+#define IGL_SQUARED_EDGE_LENGTHS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Constructs a list of squared lengths of edges opposite each index in a face
+  // (triangle/tet) list
+  //
+  // Templates:
+  //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
+  //   DerivedF derived from face indices matrix type: i.e. MatrixXi
+  //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
+  // Inputs:
+  //   V  eigen matrix #V by 3
+  //   F  #F by 2 list of mesh edges
+  //    or
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  //    or
+  //   T  #T by 4 list of mesh elements (must be tets)
+  // Outputs:
+  //   L  #F by {1|3|6} list of edge lengths squared
+  //     for edges, column of lengths
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //     for tets, columns correspond to edges
+  //     [3 0],[3 1],[3 2],[1 2],[2 0],[0 1]
+  //
+  template <typename DerivedV, typename DerivedF, typename DerivedL>
+  IGL_INLINE void squared_edge_lengths(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedL>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "squared_edge_lengths.cpp"
+#endif
+
+#endif
+

include/igl/unique_edge_map.cpp

@@ -45,6 +45,10 @@ IGL_INLINE void igl::unique_edge_map(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
+// generated by autoexplicit.sh
+template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::unique_edge_map<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);

include/igl/viewer/Viewer.cpp

@@ -80,7 +80,7 @@ static igl::viewer::Viewer * __viewer;
 static double highdpi = 1;
 static double scroll_x = 0;
 static double scroll_y = 0;
-static int global_KMod = 0;
+//static int global_KMod = 0;
 
 static void glfw_mouse_press(GLFWwindow* window, int button, int action, int modifier)
 {

include/igl/viewer/ViewerCore.cpp

@@ -383,8 +383,8 @@ IGL_INLINE void igl::viewer::ViewerCore::draw_buffer(ViewerData& data,
   assert(R.rows() == G.rows() && G.rows() == B.rows() && B.rows() == A.rows());
   assert(R.cols() == G.cols() && G.cols() == B.cols() && B.cols() == A.cols());
 
-  int x = R.rows();
-  int y = R.cols();
+  unsigned x = R.rows();
+  unsigned y = R.cols();
 
   // Create frame buffer
   GLuint frameBuffer;

include/igl/voxel_grid.cpp

@@ -22,7 +22,7 @@ IGL_INLINE void igl::voxel_grid(
   {
     if(i!=si)
     {
-      side(i) = ceil(s * (box.max()(i)-box.min()(i))/s_len);
+      side(i) = std::ceil(s * (box.max()(i)-box.min()(i))/s_len);
     }
   }
   side.array() += 2*pad_count;

include/igl/writeOFF.cpp

@@ -33,6 +33,47 @@ IGL_INLINE bool igl::writeOFF(
   return true;
 }
 
+// write mesh and colors-by-vertex to an ascii off file
+template <typename DerivedV, typename DerivedF, typename DerivedC>
+IGL_INLINE bool igl::writeOFF(
+  const std::string fname,
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedC>& C)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(V.cols() == 3 && "V should have 3 columns");
+  assert(C.cols() == 3 && "C should have 3 columns");
+
+  if(V.rows() != C.rows())
+  {
+    fprintf(stderr,"IOError: writeOFF() Only color per vertex supported. V and C should have same size.\n");
+    return false;
+  }
+
+  ofstream s(fname);
+  if(!s.is_open())
+  {
+    fprintf(stderr,"IOError: writeOFF() could not open %s\n",fname.c_str());
+    return false;
+  }
+
+  //Check if RGB values are in the range [0..1] or [0..255]
+  int rgbScale = (C.maxCoeff() <= 1.0)?255:1;
+  Eigen::MatrixXd RGB = rgbScale * C;
+
+  s<< "COFF\n"<<V.rows()<<" "<<F.rows()<<" 0\n";
+  for (unsigned i=0; i< V.rows(); i++)
+  {
+    s <<V.row(i).format(IOFormat(FullPrecision,DontAlignCols," "," ","","",""," "));
+    s << unsigned(RGB(i,0)) << " " << unsigned(RGB(i,1)) << " " << unsigned(RGB(i,2)) << " 255\n";
+  }
+
+  s<<(F.array()).format(IOFormat(FullPrecision,DontAlignCols," ","\n","3 ","","","\n"));
+  return true;
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
@@ -44,4 +85,6 @@ template bool igl::writeOFF<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix
 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&);
+template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -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&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&);
+template bool igl::writeOFF<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -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<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&);
 #endif

include/igl/writeOFF.h

@@ -14,6 +14,28 @@
 
 namespace igl 
 {
+  //Export geometry and colors-by-vertex
+  // Export a mesh from an ascii OFF file, filling in vertex positions.
+  // Only triangle meshes are supported
+  //
+  // Templates:
+  //   Scalar  type for positions and vectors (will be read as double and cast
+  //     to Scalar)
+  //   Index  type for indices (will be read as int and cast to Index)
+  // Inputs:
+  //  str  path to .off output file
+  //   V  #V by 3 mesh vertex positions
+  //   F  #F by 3 mesh indices into V
+  //   C  double matrix of rgb values per vertex #V by 3
+  // Outputs:
+  // Returns true on success, false on errors
+  template <typename DerivedV, typename DerivedF, typename DerivedC>
+  IGL_INLINE bool writeOFF(
+    const std::string str,
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedC>& C);
+
   template <typename DerivedV, typename DerivedF>
   IGL_INLINE bool writeOFF(
     const std::string str,

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

optional/CMakeLists.txt

@@ -22,7 +22,6 @@ option(LIBIGL_USE_STATIC_LIBRARY "Use LibIGL as static library" ON)
 
 
 option(LIBIGL_WITH_ANTTWEAKBAR      "Use AntTweakBar"    ON)
-option(LIBIGL_WITH_BBW              "Use BBW"            ON)
 find_package(CGAL QUIET)
 option(LIBIGL_WITH_CGAL             "Use CGAL"           "${CGAL_FOUND}")
 option(LIBIGL_WITH_COMISO           "Use CoMiso"         ON)

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

@@ -7,6 +7,7 @@
 
 #include <igl/embree/ambient_occlusion.h>
 #include <igl/embree/reorient_facets_raycast.h>
+#include <igl/embree/line_mesh_intersection.h>
 
 
 void python_export_igl_embree(py::module &me) {
@@ -16,5 +17,6 @@ void python_export_igl_embree(py::module &me) {
 
   #include "../py_igl/embree/py_ambient_occlusion.cpp"
   #include "../py_igl/embree/py_reorient_facets_raycast.cpp"
+  #include "../py_igl/embree/py_line_mesh_intersection.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

@@ -76,6 +76,59 @@ const char *__doc_igl_barycentric_coordinates = R"igl_Qu8mg5v7(// Compute baryce
   // Outputs:
   //   L  #P by 4 list of barycentric coordinates
   //   )igl_Qu8mg5v7";
+const char *__doc_igl_barycentric_to_global = R"igl_Qu8mg5v7(// Converts barycentric coordinates in the embree form to 3D coordinates
+  // Embree stores barycentric coordinates as triples: fid, bc1, bc2
+  // fid is the id of a face, bc1 is the displacement of the point wrt the 
+  // first vertex v0 and the edge v1-v0. Similarly, bc2 is the displacement
+  // wrt v2-v0.
+  // 
+  // Input:
+  // V:  #Vx3 Vertices of the mesh
+  // F:  #Fxe Faces of the mesh
+  // bc: #Xx3 Barycentric coordinates, one row per point
+  //
+  // 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`).
   //
@@ -132,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
@@ -370,6 +430,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:
@@ -398,6 +484,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
   //
@@ -455,6 +550,23 @@ const char *__doc_igl_embree_reorient_facets_raycast = R"igl_Qu8mg5v7(// Orient
     // Outputs:
     //   I  #F list of whether face has been flipped
     //   C  #F list of patch ID (output of bfs_orient > manifold patches))igl_Qu8mg5v7";
+const char *__doc_igl_embree_line_mesh_intersection = R"igl_Qu8mg5v7(// Project the point cloud V_source onto the triangle mesh
+    // V_target,F_target. 
+    // A ray is casted for every vertex in the direction specified by 
+    // N_source and its opposite.
+    //
+    // Input:
+    // V_source: #Vx3 Vertices of the source mesh
+    // N_source: #Vx3 Normals of the point cloud
+    // V_target: #V2x3 Vertices of the target mesh
+    // F_target: #F2x3 Faces of the target mesh
+    //
+    // Output:
+    // #Vx3 matrix of baricentric coordinate. Each row corresponds to 
+    // a vertex of the projected mesh and it has the following format:
+    // id b1 b2. id is the id of a face of the source mesh. b1 and b2 are 
+    // the barycentric coordinates wrt the first two edges of the triangle
+    // To convert to standard global coordinates, see barycentric_to_global.h)igl_Qu8mg5v7";
 const char *__doc_igl_find_cross_field_singularities = R"igl_Qu8mg5v7(// Inputs:
   //   V                #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F                #F by 3 eigen Matrix of face (quad) indices
@@ -493,6 +605,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).
   //
@@ -591,6 +715,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
@@ -702,6 +859,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:
@@ -894,6 +1067,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)
   // 

python/py_doc.h

@@ -4,11 +4,15 @@ extern const char *__doc_igl_arap_solve;
 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;
@@ -23,19 +27,25 @@ 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;
 extern const char *__doc_igl_embree_ambient_occlusion;
 extern const char *__doc_igl_embree_reorient_facets_raycast;
+extern const char *__doc_igl_embree_line_mesh_intersection;
 extern const char *__doc_igl_find_cross_field_singularities;
 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;
@@ -45,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;
@@ -53,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;
@@ -71,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;

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,11 +12,14 @@
 #include <igl/avg_edge_length.h>
 #include <igl/barycenter.h>
 #include <igl/barycentric_coordinates.h>
+#include <igl/barycentric_to_global.h>
+#include <igl/boundary_conditions.h>
 #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>
@@ -23,13 +27,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>
@@ -39,12 +48,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>
@@ -60,6 +72,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>
@@ -82,6 +95,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"
@@ -91,11 +106,14 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_avg_edge_length.cpp"
 #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"
@@ -103,13 +121,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"
@@ -119,12 +142,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"
@@ -140,6 +166,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/embree/py_line_mesh_intersection.cpp

@@ -0,0 +1,13 @@
+
+
+m.def("line_mesh_intersection", []
+(
+  const Eigen::MatrixXd& V_source,
+  const Eigen::MatrixXd& N_source,
+  const Eigen::MatrixXd& V_target,
+  const Eigen::MatrixXi& F_target
+)
+{
+  return igl::embree::line_mesh_intersection(V_source, N_source, V_target, F_target);
+}, __doc_igl_embree_line_mesh_intersection,
+py::arg("V_source"), py::arg("N_source"), py::arg("V_target"), py::arg("F_target"));

python/py_igl/py_barycentric_to_global.cpp

@@ -0,0 +1,12 @@
+
+
+m.def("barycentric_to_global", []
+(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::MatrixXd& bc
+)
+{
+  return igl::barycentric_to_global(V, F, bc);
+}, __doc_igl_barycentric_to_global,
+py::arg("V"), py::arg("F"), py::arg("bc"));

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_internal_angles.cpp

@@ -11,13 +11,23 @@ m.def("internal_angles", []
 }, __doc_igl_internal_angles,
 py::arg("V"), py::arg("F"), py::arg("K"));
 
-m.def("internal_angles", []
+m.def("internal_angles_using_squared_edge_lengths", []
+(
+  const Eigen::MatrixXd& L_sq,
+  Eigen::MatrixXd& K
+)
+{
+  return igl::internal_angles_using_squared_edge_lengths(L_sq, K);
+}, __doc_igl_internal_angles,
+py::arg("L_sq"), py::arg("K"));
+
+m.def("internal_angles_using_edge_lengths", []
 (
   const Eigen::MatrixXd& L,
   Eigen::MatrixXd& K
 )
 {
-  return igl::internal_angles(L, K);
+  return igl::internal_angles_using_edge_lengths(L, K);
 }, __doc_igl_internal_angles,
 py::arg("L"), py::arg("K"));
 

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) {
@@ -57,11 +61,15 @@ PYBIND11_PLUGIN(pyigl) {
            avg_edge_length
            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
@@ -77,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
@@ -86,6 +98,7 @@ PYBIND11_PLUGIN(pyigl) {
            find_cross_field_singularities
            fit_rotations
            floor
+           forward_kinematics
            gaussian_curvature
            get_seconds
            grad
@@ -95,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
@@ -118,6 +134,7 @@ PYBIND11_PLUGIN(pyigl) {
            readMESH
            readOBJ
            readOFF
+           readTGF
            read_triangle_mesh
            remove_duplicate_vertices
            rotate_vectors
@@ -168,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

@@ -4,7 +4,6 @@ project(libigl)
 ### Available options ###
 option(LIBIGL_USE_STATIC_LIBRARY    "Use libigl as static library" OFF)
 option(LIBIGL_WITH_ANTTWEAKBAR      "Use AntTweakBar"    OFF)
-option(LIBIGL_WITH_BBW              "Use BBW"            OFF)
 option(LIBIGL_WITH_CGAL             "Use CGAL"           OFF)
 option(LIBIGL_WITH_COMISO           "Use CoMiso"         OFF)
 option(LIBIGL_WITH_CORK             "Use Cork"           OFF)
@@ -154,20 +153,6 @@ if(LIBIGL_WITH_ANTTWEAKBAR)
   endif()
 endif()
 
-### Compile the BBW part ###
-if(LIBIGL_WITH_BBW)
-  if(LIBIGL_USE_STATIC_LIBRARY)
-    CompileIGL_Module("bbw")
-    if(LIBIGL_WITH_MOSEK)
-      find_package(MOSEK REQUIRED)
-      list(APPEND LIBIGL_INCLUDE_DIRS ${MOSEK_INCLUDE_DIR})
-      target_include_directories(iglbbw PRIVATE ${MOSEK_INCLUDE_DIR})
-    else()
-      target_compile_definitions(iglbbw PRIVATE -DIGL_NO_MOSEK)
-    endif()
-  endif()
-endif()
-
 ### Compile the cgal parts ###
 if(LIBIGL_WITH_CGAL) # to be cleaned
   find_package(CGAL REQUIRED)
@@ -387,6 +372,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 +401,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
@@ -472,7 +458,6 @@ if(NOT ${CMAKE_PROJECT_NAME} STREQUAL ${PROJECT_NAME})
   set(LIBIGL_INCLUDE_DIRS ${LIBIGL_INCLUDE_DIRS} PARENT_SCOPE)
   set(LIBIGL_LIBRARIES ${LIBIGL_LIBRARIES} PARENT_SCOPE)
   set(LIBIGL_ANTTWEAKBAR_EXTRA_LIBRARIES ${LIBIGL_ANTTWEAKBAR_EXTRA_LIBRARIES} PARENT_SCOPE)
-  #set(LIBIGL_BBW_EXTRA_LIBRARIES         ${LIBIGL_BBW_EXTRA_LIBRARIES}         PARENT_SCOPE)
   set(LIBIGL_CGAL_EXTRA_LIBRARIES        ${LIBIGL_CGAL_EXTRA_LIBRARIES}        PARENT_SCOPE)
   set(LIBIGL_COMISO_EXTRA_LIBRARIES      ${LIBIGL_COMISO_EXTRA_LIBRARIES}      PARENT_SCOPE)
   set(LIBIGL_CORK_EXTRA_LIBRARIES        ${LIBIGL_CORK_EXTRA_LIBRARIES}        PARENT_SCOPE)

tutorial/403_BoundedBiharmonicWeights/main.cpp

@@ -21,7 +21,7 @@
 #include <igl/readMESH.h>
 #include <igl/readTGF.h>
 #include <igl/viewer/Viewer.h>
-#include <igl/bbw/bbw.h>
+#include <igl/bbw.h>
 //#include <igl/embree/bone_heat.h>
 
 #include <Eigen/Geometry>
@@ -139,11 +139,11 @@ int main(int argc, char *argv[])
   igl::boundary_conditions(V,T,C,VectorXi(),BE,MatrixXi(),b,bc);
 
   // compute BBW weights matrix
-  igl::bbw::BBWData bbw_data;
+  igl::BBWData bbw_data;
   // only a few iterations for sake of demo
   bbw_data.active_set_params.max_iter = 8;
   bbw_data.verbosity = 2;
-  if(!igl::bbw::bbw(V,T,b,bc,bbw_data,W))
+  if(!igl::bbw(V,T,b,bc,bbw_data,W))
   {
     return false;
   }

tutorial/703_Decimation/main.cpp

@@ -67,6 +67,7 @@ int main(int argc, char * argv[])
 
     C.resize(E.rows(),V.cols());
     VectorXd costs(E.rows());
+    Q.clear();
     for(int e = 0;e<E.rows();e++)
     {
       double cost = e;

tutorial/709_VectorFieldVisualizer/main.cpp

@@ -10,7 +10,6 @@
 #include <igl/slice.h>
 #include <igl/sort_vectors_ccw.h>
 #include <igl/streamlines.h>
-#include <igl/triangle_triangle_adjacency.h>
 #include <igl/copyleft/comiso/nrosy.h>
 #include <igl/viewer/Viewer.h>
 

tutorial/CMakeLists.txt

@@ -7,7 +7,6 @@ option(LIBIGL_WITH_VIEWER      "Use OpenGL viewer"  ON)
 option(LIBIGL_WITH_NANOGUI     "Use Nanogui menu"   OFF)
 
 ### libIGL options: choose your dependencies (by default everything is OFF, in this example we need the viewer) ###
-option(LIBIGL_WITH_BBW              "Use BBW"            ON)
 find_package(CGAL QUIET)
 option(LIBIGL_WITH_CGAL             "Use CGAL"           "${CGAL_FOUND}")
 option(LIBIGL_WITH_COMISO           "Use CoMiso"         ON)
@@ -114,9 +113,7 @@ endif()
 if(TUTORIALS_CHAPTER4)
   add_subdirectory("401_BiharmonicDeformation")
   add_subdirectory("402_PolyharmonicDeformation")
-  if(LIBIGL_WITH_BBW)
-    add_subdirectory("403_BoundedBiharmonicWeights")
-  endif()
+  add_subdirectory("403_BoundedBiharmonicWeights")
   add_subdirectory("404_DualQuaternionSkinning")
   add_subdirectory("405_AsRigidAsPossible")
   add_subdirectory("406_FastAutomaticSkinningTransformations")

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-02cf5efb916ce8e885dcaeee5a89f91fa4e1f440
+cd11b42bda4357cfaae8f62d5036f840df5b627d

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-521b2da65bd08996d25fb1fab30ff66f704da312
+d17b6043e03c640488d1b324beebc066cdc3c830