merge

Former-commit-id: 38a012d8177c490cf7d787d0ac345105941cb9f8

.appveyor.yml

@@ -5,6 +5,7 @@ clone_folder: C:\projects\libigl
 branches:
   only:
     - master
+    - alecjacobson
 install:
   - git submodule update --init --recursive
   - cinstall: python

README.md

@@ -213,7 +213,8 @@ few labs/companies/institutions using libigl:
  - ETH Zurich, [Interactive Geometry Lab](http://igl.ethz.ch/) and [Advanced Technologies Lab](http://ait.inf.ethz.ch/), Swizterland
  - George Mason University, [CraGL](http://cs.gmu.edu/~ygingold/), USA
  - [Hong Kong University of Science and Technology](http://www.ust.hk/), Hong Kong
- - [Inria](Université Grenoble Alpes), France
+ - [Inria, Université Grenoble Alpes](https://www.inria.fr/centre/grenoble/), France
+ - [Jiangnan university](http://english.jiangnan.edu.cn), China
  - [National Institute of Informatics](http://www.nii.ac.jp/en/), Japan
  - New York University, [Media Research Lab](http://mrl.nyu.edu/), USA
  - NYUPoly, [Game Innovation Lab](http://game.engineering.nyu.edu/), USA
@@ -239,7 +240,7 @@ few labs/companies/institutions using libigl:
 ## Contact
 
 Libigl is a group endeavor led by [Alec
-Jacobson](http://www.cs.columbia.edu/~jacobson/) and [Daniele
+Jacobson](http://www.cs.toronto.edu/~jacobson/) and [Daniele
 Panozzo](http://www.inf.ethz.ch/personal/dpanozzo/). Please [contact
 us](mailto:alecjacobson@gmail.com,daniele.panozzo@gmail.com) if you have
 questions or comments. For troubleshooting, please post an

before-submitting-pull-request.md

@@ -54,7 +54,7 @@ and paste it at the bottom of the relevant .cpp file with the word template in f
 
 ```cpp
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::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> >&);
 #endif
 ```

include/igl/AABB.cpp

@@ -11,7 +11,6 @@
 #include "barycentric_coordinates.h"
 #include "colon.h"
 #include "doublearea.h"
-#include "matlab_format.h"
 #include "point_simplex_squared_distance.h"
 #include "project_to_line_segment.h"
 #include "sort.h"
@@ -391,15 +390,16 @@ 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 = [&]()
     {
       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);
+      Scalar sqr_d_right = 
+        m_right->squared_distance(V,Ele,p,sqr_d,i_right,c_right);
+      this->set_min(p,sqr_d_right,i_right,c_right,sqr_d,i,c);
       looked_right = true;
     };
 
@@ -413,8 +413,10 @@ igl::AABB<DerivedV,DIM>::squared_distance(
       look_right();
     }
     // if haven't looked left and could be less than current min, then look
-    Scalar  left_min_sqr_d = m_left->m_box.squaredExteriorDistance(p.transpose());
-    Scalar right_min_sqr_d = m_right->m_box.squaredExteriorDistance(p.transpose());
+    Scalar left_min_sqr_d = 
+      m_left->m_box.squaredExteriorDistance(p.transpose());
+    Scalar right_min_sqr_d = 
+      m_right->m_box.squaredExteriorDistance(p.transpose());
     if(left_min_sqr_d < right_min_sqr_d)
     {
       if(!looked_left && left_min_sqr_d<sqr_d)
@@ -457,7 +459,9 @@ IGL_INLINE void igl::AABB<DerivedV,DIM>::squared_distance(
   assert(P.cols() == V.cols() && "cols in P should match dim of cols in V");
   sqrD.resize(P.rows(),1);
   I.resize(P.rows(),1);
-  C.resize(P.rows(),P.cols());
+  C.resizeLike(P);
+  // O( #P * log #Ele ), where log #Ele is really the depth of this AABB
+  // hierarchy
   for(int p = 0;p<P.rows();p++)
   {
     RowVectorDIMS Pp = P.row(p), c;
@@ -507,7 +511,8 @@ template <
   typename DerivedsqrD, 
   typename DerivedI, 
   typename DerivedC>
-IGL_INLINE typename igl::AABB<DerivedV,DIM>::Scalar igl::AABB<DerivedV,DIM>::squared_distance_helper(
+IGL_INLINE typename igl::AABB<DerivedV,DIM>::Scalar 
+  igl::AABB<DerivedV,DIM>::squared_distance_helper(
   const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::MatrixXi & Ele, 
   const AABB<Derivedother_V,DIM> * other,
@@ -947,7 +952,13 @@ igl::AABB<DerivedV,DIM>::intersect_ray(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&) const;
+// generated by autoexplicit.sh
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -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::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&) const;
+// generated by autoexplicit.sh
+template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::squared_distance<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -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::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&) const;
 // generated by autoexplicit.sh
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::init(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 2>::init(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&);

include/igl/AABB.h

@@ -177,12 +177,10 @@ public:
       // Inputs:
       //   V  #V by dim list of vertex positions
       //   Ele  #Ele by dim list of simplex indices
-      //   P  3 list of query point coordinates
-      //   min_sqr_d  current minimum squared distance (only find distances
-      //   less than this)
+      //   p  dim-long query point 
       // Outputs:
-      //   I  #P list of facet indices corresponding to smallest distances
-      //   C  #P by 3 list of closest points
+      //   i  facet index corresponding to smallest distances
+      //   c  closest point
       // Returns squared distance
       //
       // Known bugs: currently assumes Elements are triangles regardless of
@@ -194,6 +192,22 @@ public:
         int & i,
         RowVectorDIMS & c) const;
 //private:
+      // Compute squared distance to a query point
+      //
+      // Inputs:
+      //   V  #V by dim list of vertex positions
+      //   Ele  #Ele by dim list of simplex indices
+      //   p  dim-long query point 
+      //   min_sqr_d  current minimum squared distance (only consider distances
+      //     less than this), see output.
+      // Outputs:
+      //   min_sqr_d  updated current minimum squared distance
+      //   i  facet index corresponding to smallest distances
+      //   c  closest point
+      // Returns squared distance
+      //
+      // Known bugs: currently assumes Elements are triangles regardless of
+      // dimension.
       IGL_INLINE Scalar squared_distance(
         const Eigen::PlainObjectBase<DerivedV> & V,
         const Eigen::MatrixXi & Ele, 
@@ -226,6 +240,18 @@ public:
 
 
 public:
+      // Compute the squared distance from all query points in P to the
+      // _closest_ points on the primitives stored in the AABB hierarchy for
+      // the mesh (V,Ele).
+      //
+      // Inputs:
+      //   V  #V by dim list of vertex positions
+      //   Ele  #Ele by dim list of simplex indices
+      //   P  #P by dim list of query points
+      // Outputs:
+      //   sqrD  #P list of squared distances
+      //   I  #P list of indices into Ele of closest primitives
+      //   C  #P by dim list of closest points
       template <
         typename DerivedP, 
         typename DerivedsqrD, 
@@ -239,6 +265,21 @@ public:
         Eigen::PlainObjectBase<DerivedI> & I,
         Eigen::PlainObjectBase<DerivedC> & C) const;
 
+      // Compute the squared distance from all query points in P already stored
+      // in its own AABB hierarchy to the _closest_ points on the primitives
+      // stored in the AABB hierarchy for the mesh (V,Ele).
+      //
+      // Inputs:
+      //   V  #V by dim list of vertex positions
+      //   Ele  #Ele by dim list of simplex indices
+      //   other  AABB hierarchy of another set of primitives (must be points)
+      //   other_V  #other_V by dim list of query points
+      //   other_Ele  #other_Ele by ss list of simplex indices into other_V
+      //     (must be simple list of points: ss == 1)
+      // Outputs:
+      //   sqrD  #P list of squared distances
+      //   I  #P list of indices into Ele of closest primitives
+      //   C  #P by dim list of closest points
       template < 
         typename Derivedother_V,
         typename DerivedsqrD, 
@@ -269,7 +310,21 @@ private:
         Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
         Eigen::PlainObjectBase<DerivedI> & I,
         Eigen::PlainObjectBase<DerivedC> & C) const;
-      // Helper function for leaves: works in-place on sqr_d
+      // Compute the squared distance to the primitive in this node: assumes
+      // that this is indeed a leaf node.
+      //
+      // Inputs:
+      //   V  #V by dim list of vertex positions
+      //   Ele  #Ele by dim list of simplex indices
+      //   p  dim-long query point
+      //   sqr_d  current minimum distance for this query, see output
+      //   i  current index into Ele of closest point, see output
+      //   c  dim-long current closest point, see output
+      // Outputs:
+      //   sqr_d   minimum of initial value and squared distance to this
+      //     primitive
+      //   i  possibly updated index into Ele of closest point
+      //   c  dim-long possibly updated closest point
       IGL_INLINE void leaf_squared_distance(
         const Eigen::PlainObjectBase<DerivedV> & V,
         const Eigen::MatrixXi & Ele, 
@@ -277,6 +332,23 @@ private:
         Scalar & sqr_d,
         int & i,
         RowVectorDIMS & c) const;
+      // If new distance (sqr_d_candidate) is less than current distance
+      // (sqr_d), then update this distance and its associated values
+      // _in-place_:
+      //
+      // Inputs:
+      //   p  dim-long query point (only used in DEBUG mode)
+      //   sqr_d  candidate minimum distance for this query, see output
+      //   i  candidate index into Ele of closest point, see output
+      //   c  dim-long candidate closest point, see output
+      //   sqr_d  current minimum distance for this query, see output
+      //   i  current index into Ele of closest point, see output
+      //   c  dim-long current closest point, see output
+      // Outputs:
+      //   sqr_d   minimum of initial value and squared distance to this
+      //     primitive
+      //   i  possibly updated index into Ele of closest point
+      //   c  dim-long possibly updated closest point
       IGL_INLINE void set_min(
         const RowVectorDIMS & p,
         const Scalar sqr_d_candidate,

include/igl/EPS.cpp

@@ -26,5 +26,5 @@ template <> IGL_INLINE double igl::EPS_SQ()
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/Timer.h

@@ -21,6 +21,7 @@
 #else
 #include <sys/time.h>
 #endif
+#include <cstddef>
 
 namespace igl
 {

include/igl/active_set.cpp

@@ -10,7 +10,7 @@
 #include "slice.h"
 #include "slice_into.h"
 #include "cat.h"
-#include "matlab_format.h"
+//#include "matlab_format.h"
 
 #include <iostream>
 #include <limits>
@@ -56,7 +56,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   // Discard const qualifiers
   //if(B.size() == 0)
   //{
-  //  B = Eigen::PlainObjectBase<DerivedB>::Zero(n,1);
+  //  B = DerivedB::Zero(n,1);
   //}
   assert(n == B.rows() && "B.rows() must match A.rows()");
   assert(B.cols() == 1 && "B must be a column vector");
@@ -71,7 +71,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   Eigen::Matrix<typename Derivedux::Scalar,Eigen::Dynamic,1> ux;
   if(p_lx.size() == 0)
   {
-    lx = Eigen::PlainObjectBase<Derivedlx>::Constant(
+    lx = Derivedlx::Constant(
       n,1,-numeric_limits<typename Derivedlx::Scalar>::max());
   }else
   {
@@ -79,7 +79,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   }
   if(p_ux.size() == 0)
   {
-    ux = Eigen::PlainObjectBase<Derivedux>::Constant(
+    ux = Derivedux::Constant(
       n,1,numeric_limits<typename Derivedux::Scalar>::max());
   }else
   {
@@ -109,8 +109,8 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   Matrix<BOOL,Dynamic,1> as_ieq = Matrix<BOOL,Dynamic,1>::Constant(Aieq.rows(),1,FALSE);
 
   // Keep track of previous Z for comparison
-  PlainObjectBase<DerivedZ> old_Z;
-  old_Z = PlainObjectBase<DerivedZ>::Constant(
+  DerivedZ old_Z;
+  old_Z = DerivedZ::Constant(
       n,1,numeric_limits<typename DerivedZ::Scalar>::max());
 
   int iter = 0;
@@ -143,7 +143,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       }
       if(Aieq.rows() > 0)
       {
-        PlainObjectBase<DerivedZ> AieqZ;
+        DerivedZ AieqZ;
         AieqZ = Aieq*Z;
         for(int a = 0;a<Aieq.rows();a++)
         {
@@ -190,9 +190,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 #endif
 
     // PREPARE FIXED VALUES
-    PlainObjectBase<Derivedknown> known_i;
+    Derivedknown known_i;
     known_i.resize(nk + as_lx_count + as_ux_count,1);
-    PlainObjectBase<DerivedY> Y_i;
+    DerivedY Y_i;
     Y_i.resize(nk + as_lx_count + as_ux_count,1);
     {
       known_i.block(0,0,known.rows(),known.cols()) = known;
@@ -225,7 +225,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     // PREPARE EQUALITY CONSTRAINTS
     VectorXi as_ieq_list(as_ieq_count,1);
     // Gather active constraints and resp. rhss
-    PlainObjectBase<DerivedBeq> Beq_i;
+    DerivedBeq Beq_i;
     Beq_i.resize(Beq.rows()+as_ieq_count,1);
     Beq_i.head(Beq.rows()) = Beq;
     {
@@ -262,7 +262,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     }
 #endif
 
-    Eigen::PlainObjectBase<DerivedZ> sol;
+    DerivedZ sol;
     if(known_i.size() == A.rows())
     {
       // Everything's fixed?
@@ -311,7 +311,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     SparseMatrix<AT> Ak;
     // Slow
     slice(A,known_i,1,Ak);
-    Eigen::PlainObjectBase<DerivedB> Bk;
+    DerivedB Bk;
     slice(B,known_i,Bk);
     MatrixXd Lambda_known_i = -(0.5*Ak*Z + 0.5*Bk);
     // reverse the lambda values for lx
@@ -364,7 +364,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template igl::SolverStatus igl::active_set<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, 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::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, igl::active_set_params const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template igl::SolverStatus igl::active_set<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, 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::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::active_set_params const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/adjacency_list.cpp

@@ -159,7 +159,9 @@ IGL_INLINE void igl::adjacency_list(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::adjacency_list<Eigen::Matrix<int, -1, 2, 0, -1, 2>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, bool);
 // generated by autoexplicit.sh
 template void igl::adjacency_list<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(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> > > >&, bool);
 template void igl::adjacency_list<Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, bool);

include/igl/adjacency_matrix.cpp

@@ -13,7 +13,7 @@
 
 template <typename DerivedF, typename T>
 IGL_INLINE void igl::adjacency_matrix(
-  const Eigen::PlainObjectBase<DerivedF> & F, 
+  const Eigen::MatrixBase<DerivedF> & F, 
   Eigen::SparseMatrix<T>& A)
 {
   using namespace std;
@@ -65,7 +65,9 @@ IGL_INLINE void igl::adjacency_matrix(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
-template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<int, 0, int>&);
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, bool>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<bool, 0, int>&);
+template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::adjacency_matrix<Eigen::Matrix<int, -1, -1, 0, -1, -1>, int>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<int, 0, int>&);
 #endif

include/igl/adjacency_matrix.h

@@ -40,7 +40,7 @@ namespace igl
   // See also: edges, cotmatrix, diag
   template <typename DerivedF, typename T>
   IGL_INLINE void adjacency_matrix(
-    const Eigen::PlainObjectBase<DerivedF> & F, 
+    const Eigen::MatrixBase<DerivedF> & F, 
     Eigen::SparseMatrix<T>& A);
 }
 

include/igl/all.cpp

@@ -0,0 +1,26 @@
+// 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 "all.h"
+#include "redux.h"
+
+
+template <typename AType, typename DerivedB>
+IGL_INLINE void igl::all(
+  const Eigen::SparseMatrix<AType> & A, 
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  typedef typename DerivedB::Scalar Scalar;
+  igl::redux(A,dim,[](Scalar a, Scalar b){ return a && b!=0;},B);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+#endif
+
+

include/igl/all.h

@@ -0,0 +1,36 @@
+// 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/.
+#ifndef IGL_ALL_H
+#define IGL_ALL_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // For Dense matrices use: A.rowwise().all() or A.colwise().all()
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   dim  dimension along which to check for all (1 or 2)
+  // Output:
+  //   B  n-long vector (if dim == 1) 
+  //   or
+  //   B  m-long vector (if dim == 2)
+  //
+  template <typename AType, typename DerivedB>
+  IGL_INLINE void all(
+    const Eigen::SparseMatrix<AType> & A, 
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "all.cpp"
+#endif
+#endif
+
+

include/igl/all_edges.cpp

@@ -13,38 +13,39 @@ IGL_INLINE void igl::all_edges(
   Eigen::PlainObjectBase<DerivedE> & E)
 {
   E.resize(F.rows()*F.cols(),F.cols()-1);
+  typedef typename DerivedE::Scalar EScalar;
   switch(F.cols())
   {
     case 4:
-      E.block(0*F.rows(),0,F.rows(),1) = F.col(1);
-      E.block(0*F.rows(),1,F.rows(),1) = F.col(3);
-      E.block(0*F.rows(),2,F.rows(),1) = F.col(2);
+      E.block(0*F.rows(),0,F.rows(),1) = F.col(1).template cast<EScalar>();
+      E.block(0*F.rows(),1,F.rows(),1) = F.col(3).template cast<EScalar>();
+      E.block(0*F.rows(),2,F.rows(),1) = F.col(2).template cast<EScalar>();
 
-      E.block(1*F.rows(),0,F.rows(),1) = F.col(0);
-      E.block(1*F.rows(),1,F.rows(),1) = F.col(2);
-      E.block(1*F.rows(),2,F.rows(),1) = F.col(3);
+      E.block(1*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(1*F.rows(),1,F.rows(),1) = F.col(2).template cast<EScalar>();
+      E.block(1*F.rows(),2,F.rows(),1) = F.col(3).template cast<EScalar>();
 
-      E.block(2*F.rows(),0,F.rows(),1) = F.col(0);
-      E.block(2*F.rows(),1,F.rows(),1) = F.col(3);
-      E.block(2*F.rows(),2,F.rows(),1) = F.col(1);
+      E.block(2*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(2*F.rows(),1,F.rows(),1) = F.col(3).template cast<EScalar>();
+      E.block(2*F.rows(),2,F.rows(),1) = F.col(1).template cast<EScalar>();
 
-      E.block(3*F.rows(),0,F.rows(),1) = F.col(0);
-      E.block(3*F.rows(),1,F.rows(),1) = F.col(1);
-      E.block(3*F.rows(),2,F.rows(),1) = F.col(2);
+      E.block(3*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(3*F.rows(),1,F.rows(),1) = F.col(1).template cast<EScalar>();
+      E.block(3*F.rows(),2,F.rows(),1) = F.col(2).template cast<EScalar>();
       return;
     case 3:
-      E.block(0*F.rows(),0,F.rows(),1) = F.col(1);
-      E.block(0*F.rows(),1,F.rows(),1) = F.col(2);
-      E.block(1*F.rows(),0,F.rows(),1) = F.col(2);
-      E.block(1*F.rows(),1,F.rows(),1) = F.col(0);
-      E.block(2*F.rows(),0,F.rows(),1) = F.col(0);
-      E.block(2*F.rows(),1,F.rows(),1) = F.col(1);
+      E.block(0*F.rows(),0,F.rows(),1) = F.col(1).template cast<EScalar>();
+      E.block(0*F.rows(),1,F.rows(),1) = F.col(2).template cast<EScalar>();
+      E.block(1*F.rows(),0,F.rows(),1) = F.col(2).template cast<EScalar>();
+      E.block(1*F.rows(),1,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(2*F.rows(),0,F.rows(),1) = F.col(0).template cast<EScalar>();
+      E.block(2*F.rows(),1,F.rows(),1) = F.col(1).template cast<EScalar>();
       return;
   }
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::all_edges<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::all_edges<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
 template void igl::all_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);

include/igl/all_pairs_distances.cpp

@@ -33,7 +33,7 @@ IGL_INLINE void igl::all_pairs_distances(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template void igl::all_pairs_distances<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, bool, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/ambient_occlusion.cpp

@@ -127,7 +127,7 @@ IGL_INLINE void igl::ambient_occlusion(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template void igl::ambient_occlusion<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -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> >&);
 // generated by autoexplicit.sh

include/igl/angle_bound_frame_fields.cpp

@@ -38,8 +38,7 @@ namespace igl {
       Eigen::VectorXi indInteriorToFull;
       Eigen::VectorXi indFullToInterior;
 
-//#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
-      Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
+      DerivedV B1, B2, FN;
 
       //laplacians
       Eigen::SparseMatrix<std::complex<typename DerivedV::Scalar>> DDA, DDB;
@@ -756,5 +755,5 @@ IGL_INLINE bool igl::angle_bound_frame_fields(const igl::AngleBoundFFSolverData<
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/any.cpp

@@ -0,0 +1,26 @@
+// 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 "any.h"
+#include "redux.h"
+
+
+template <typename AType, typename DerivedB>
+IGL_INLINE void igl::any(
+  const Eigen::SparseMatrix<AType> & A, 
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  typedef typename DerivedB::Scalar Scalar;
+  igl::redux(A,dim,[](Scalar a, Scalar b){ return a || b!=0;},B);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::any<bool, Eigen::Array<bool, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<bool, 0, int> const&, int, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/any.h

@@ -0,0 +1,35 @@
+// 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/.
+#ifndef IGL_ANY_H
+#define IGL_ANY_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // For Dense matrices use: A.rowwise().any() or A.colwise().any()
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   dim  dimension along which to check for any (1 or 2)
+  // Output:
+  //   B  n-long vector (if dim == 1) 
+  //   or
+  //   B  m-long vector (if dim == 2)
+  //
+  template <typename AType, typename DerivedB>
+  IGL_INLINE void any(
+    const Eigen::SparseMatrix<AType> & A, 
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "any.cpp"
+#endif
+#endif
+

include/igl/any_of.cpp

@@ -14,7 +14,7 @@ IGL_INLINE bool igl::any_of(const Mat & S)
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template bool igl::any_of<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 #endif
 

include/igl/arap.cpp

@@ -55,8 +55,8 @@ IGL_INLINE bool igl::arap_precomputation(
   assert(data.dim <= V.cols() && "solve dim should be <= embedding");
   bool flat = (V.cols() - data.dim)==1;
 
-  PlainObjectBase<DerivedV> plane_V;
-  PlainObjectBase<DerivedF> plane_F;
+  DerivedV plane_V;
+  DerivedF plane_F;
   typedef SparseMatrix<Scalar> SparseMatrixS;
   SparseMatrixS ref_map,ref_map_dim;
   if(flat)

include/igl/arap_dof.cpp

@@ -149,7 +149,7 @@ IGL_INLINE bool igl::arap_dof_precomputation(
   // Apply group sum to each dimension's block of covariance scatter matrix
   SparseMatrix<double> G_sum_dim;
   repdiag(G_sum,data.dim,G_sum_dim);
-  CSM = G_sum_dim * CSM;
+  CSM = (G_sum_dim * CSM).eval();
 #ifdef EXTREME_VERBOSE
   cout<<"CSMIJV=["<<endl;print_ijv(CSM,1);cout<<endl<<"];"<<
     endl<<"CSM=sparse(CSMIJV(:,1),CSMIJV(:,2),CSMIJV(:,3),"<<
@@ -874,7 +874,7 @@ IGL_INLINE bool igl::arap_dof_update(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit instanciation
+// Explicit template instantiation
 template bool igl::arap_dof_update<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>(ArapDOFData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, int, double, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 template bool igl::arap_dof_recomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>(Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int> const&, ArapDOFData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>&);
 template bool igl::arap_dof_precomputation<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, ArapDOFData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, double>&);

include/igl/arap_linear_block.cpp

@@ -248,6 +248,6 @@ IGL_INLINE void igl::arap_linear_block_elements(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template IGL_INLINE void igl::arap_linear_block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, int, igl::ARAPEnergyType, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/average_onto_faces.cpp

@@ -25,5 +25,5 @@ IGL_INLINE void igl::average_onto_faces(const Eigen::Matrix<T, Eigen::Dynamic, E
 };
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/average_onto_vertices.cpp

@@ -29,5 +29,5 @@ IGL_INLINE void igl::average_onto_vertices(const Eigen::MatrixBase<DerivedV> &V,
 };
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/avg_edge_length.cpp

@@ -30,7 +30,7 @@ IGL_INLINE double igl::avg_edge_length(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template double igl::avg_edge_length<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template double igl::avg_edge_length<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);

include/igl/axis_angle_to_quat.cpp

@@ -34,7 +34,7 @@ IGL_INLINE void igl::axis_angle_to_quat(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template void igl::axis_angle_to_quat<double>(double const*, double, double*);
 // generated by autoexplicit.sh

include/igl/barycenter.cpp

@@ -32,7 +32,7 @@ IGL_INLINE void igl::barycenter(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit instanciation
+// Explicit template instantiation
 template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(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, 4, 0, -1, 4> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
 template void igl::barycenter<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> >&);

include/igl/barycentric_coordinates.cpp

@@ -7,7 +7,6 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "barycentric_coordinates.h"
 #include "volume.h"
-#include "doublearea.h"
 
 template <
   typename DerivedP,

include/igl/bbw.cpp

@@ -137,7 +137,7 @@ IGL_INLINE bool igl::bbw(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 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/bbw.cpp

@@ -1,201 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
-// obtain one at http://mozilla.org/MPL/2.0/.
-#include "bbw.h"
-
-#include <igl/cotmatrix.h>
-#include <igl/massmatrix.h>
-#include <igl/invert_diag.h>
-#include <igl/speye.h>
-#include <igl/slice_into.h>
-#include <igl/min_quad_with_fixed.h>
-
-#include <Eigen/Sparse>
-
-#include <iostream>
-#include <cstdio>
-
-igl::bbw::BBWData::BBWData():
-  partition_unity(false),
-  W0(),
-#ifndef IGL_NO_MOSEK
-  mosek_data(),
-#endif
-  active_set_params(),
-  qp_solver(QP_SOLVER_IGL_ACTIVE_SET),
-  verbosity(0)
-{
-  // We know that the Bilaplacian is positive semi-definite
-  active_set_params.Auu_pd = true;
-}
-
-void igl::bbw::BBWData::print()
-{
-  using namespace std;
-  cout<<"partition_unity: "<<partition_unity<<endl;
-  cout<<"W0=["<<endl<<W0<<endl<<"];"<<endl;
-  cout<<"qp_solver: "<<QPSolverNames[qp_solver]<<endl;
-}
-
-
-template <
-  typename DerivedV,
-  typename DerivedEle,
-  typename Derivedb,
-  typename Derivedbc,
-  typename DerivedW>
-IGL_INLINE bool igl::bbw::bbw(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedEle> & Ele,
-  const Eigen::PlainObjectBase<Derivedb> & b,
-  const Eigen::PlainObjectBase<Derivedbc> & bc,
-  igl::bbw::BBWData & data,
-  Eigen::PlainObjectBase<DerivedW> & W
-  )
-{
-  using namespace std;
-  using namespace Eigen;
-
-  // number of domain vertices
-  int n = V.rows();
-  // number of handles
-  int m = bc.cols();
-  // Build biharmonic operator
-  SparseMatrix<typename DerivedW::Scalar> L;
-  cotmatrix(V,Ele,L);
-  SparseMatrix<typename DerivedW::Scalar> M;
-  SparseMatrix<typename DerivedW::Scalar> Mi;
-  massmatrix(V,Ele,MASSMATRIX_TYPE_DEFAULT,M);
-  invert_diag(M,Mi);
-  SparseMatrix<typename DerivedW::Scalar> Q = L.transpose() * Mi * L;
-  assert(!data.partition_unity && "partition_unity not implemented yet");
-
-  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 uc(0,1),Beq(0,1),Bieq(0,1),lc(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;
-    switch(data.qp_solver)
-    {
-      case QP_SOLVER_IGL_ACTIVE_SET:
-      {
-        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
-#pragma omp parallel for
-        for(int i = 0;i<m;i++)
-        {
-          // Quicker exit for openmp
-          if(error)
-          {
-            continue;
-          }
-          if(data.verbosity >= 1)
-          {
-#pragma omp 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;
-        }
-        if(error)
-        {
-          return false;
-        }
-        break;
-      }
-      case QP_SOLVER_MOSEK:
-      {
-#ifdef IGL_NO_MOSEK
-        assert(false && "Use another QPSolver. Recompile without IGL_NO_MOSEK defined.");
-        cerr<<"Use another QPSolver. Recompile without IGL_NO_MOSEK defined."<<endl;
-        return false;
-#else
-        // 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,data.mosek_data,Wi);
-          if(!r)
-          {
-            return false;
-          }
-          W.col(i) = Wi;
-        }
-#endif
-        break;
-      }
-      default:
-      {
-        assert(false && "Unknown qp_solver");
-        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::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::bbw::BBWData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-#endif

include/igl/bbw/bbw.h

@@ -1,99 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-//
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-//
-// This Source Code Form is subject to the terms of the Mozilla Public License
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_BBW_BBW_H
-#define IGL_BBW_BBW_H
-#include "../igl_inline.h"
-
-#include <Eigen/Dense>
-#ifndef IGL_NO_MOSEK
-#  include <igl/mosek/mosek_quadprog.h>
-#endif
-#include <igl/active_set.h>
-
-namespace igl
-{
-  namespace bbw
-  {
-    enum QPSolver
-    {
-      QP_SOLVER_IGL_ACTIVE_SET = 0,
-      QP_SOLVER_MOSEK = 1,
-      NUM_QP_SOLVERS = 2
-    };
-    const char * const QPSolverNames[NUM_QP_SOLVERS] =
-    {
-      "QP_SOLVER_IGL_ACTIVE_SET",
-      "QP_SOLVER_MOSEK"
-    };
-    // 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;
-#ifndef IGL_NO_MOSEK
-        igl::mosek::MosekData mosek_data;
-#endif
-        igl::active_set_params active_set_params;
-        // Which solver
-        QPSolver qp_solver;
-        // 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/bfs.cpp

@@ -0,0 +1,93 @@
+#include "bfs.h"
+#include "list_to_matrix.h"
+#include <vector>
+#include <queue>
+
+template <
+  typename AType,
+  typename DerivedD,
+  typename DerivedP>
+IGL_INLINE void igl::bfs(
+  const AType & A,
+  const size_t s,
+  Eigen::PlainObjectBase<DerivedD> & D,
+  Eigen::PlainObjectBase<DerivedP> & P)
+{
+  std::vector<typename DerivedD::Scalar> vD;
+  std::vector<typename DerivedP::Scalar> vP;
+  bfs(A,s,vD,vP);
+  list_to_matrix(vD,D);
+  list_to_matrix(vP,P);
+}
+
+template <
+  typename AType,
+  typename DType,
+  typename PType>
+IGL_INLINE void igl::bfs(
+  const std::vector<std::vector<AType> > & A,
+  const size_t s,
+  std::vector<DType> & D,
+  std::vector<PType> & P)
+{
+  // number of nodes
+  int N = s+1;
+  for(const auto & Ai : A) for(const auto & a : Ai) N = std::max(N,a+1);
+  std::vector<bool> seen(N,false);
+  P.resize(N,-1);
+  std::queue<std::pair<int,int> > Q;
+  Q.push({s,-1});
+  while(!Q.empty())
+  {
+    const int f = Q.front().first;
+    const int p = Q.front().second;
+    Q.pop();
+    if(seen[f])
+    {
+      continue;
+    }
+    D.push_back(f);
+    P[f] = p;
+    seen[f] = true;
+    for(const auto & n : A[f]) Q.push({n,f});
+  }
+}
+
+
+template <
+  typename AType,
+  typename DType,
+  typename PType>
+IGL_INLINE void igl::bfs(
+  const Eigen::SparseMatrix<AType> & A,
+  const size_t s,
+  std::vector<DType> & D,
+  std::vector<PType> & P)
+{
+  // number of nodes
+  int N = A.rows();
+  assert(A.rows() == A.cols());
+  std::vector<bool> seen(N,false);
+  P.resize(N,-1);
+  std::queue<std::pair<int,int> > Q;
+  Q.push({s,-1});
+  while(!Q.empty())
+  {
+    const int f = Q.front().first;
+    const int p = Q.front().second;
+    Q.pop();
+    if(seen[f])
+    {
+      continue;
+    }
+    D.push_back(f);
+    P[f] = p;
+    seen[f] = true;
+    for(typename Eigen::SparseMatrix<AType>::InnerIterator it (A,f); it; ++it)
+    {
+      if(it.value()) Q.push({it.index(),f});
+    }
+  }
+
+}
+

include/igl/bfs.h

@@ -0,0 +1,54 @@
+#ifndef IGL_BFS_H
+#define IGL_BFS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Traverse a **directed** graph represented by an adjacency list using
+  // breadth first search
+  //
+  // Inputs:
+  //   A  #V list of adjacency lists  or #V by #V adjacency matrix
+  //   s  starting node (index into A)
+  // Outputs:
+  //   D  #V list of indices into rows of A in the order in which graph nodes
+  //     are discovered.
+  //   P  #V list of indices into rows of A of predecessor in resulting
+  //     spanning tree {-1 indicates root/not discovered), order corresponds to
+  //     V **not** D.
+  template <
+    typename AType,
+    typename DerivedD,
+    typename DerivedP>
+  IGL_INLINE void bfs(
+    const AType & A,
+    const size_t s,
+    Eigen::PlainObjectBase<DerivedD> & D,
+    Eigen::PlainObjectBase<DerivedP> & P);
+
+  template <
+    typename AType,
+    typename DType,
+    typename PType>
+  IGL_INLINE void bfs(
+    const std::vector<std::vector<AType> > & A,
+    const size_t s,
+    std::vector<DType> & D,
+    std::vector<PType> & P);
+  template <
+    typename AType,
+    typename DType,
+    typename PType>
+  IGL_INLINE void bfs(
+    const Eigen::SparseMatrix<AType> & A,
+    const size_t s,
+    std::vector<DType> & D,
+    std::vector<PType> & P);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "bfs.cpp"
+#endif
+#endif
+

include/igl/bfs_orient.cpp

@@ -95,6 +95,6 @@ IGL_INLINE void igl::bfs_orient(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::bfs_orient<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<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> >&);
 #endif

include/igl/biharmonic_coordinates.cpp

@@ -146,6 +146,6 @@ IGL_INLINE bool igl::biharmonic_coordinates(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template bool igl::biharmonic_coordinates<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, int, 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&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/bijective_composite_harmonic_mapping.cpp

@@ -0,0 +1,115 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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 "bijective_composite_harmonic_mapping.h"
+
+#include "slice.h"
+#include "doublearea.h"
+#include "harmonic.h"
+//#include "matlab/MatlabWorkspace.h"
+#include <iostream>
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedU>
+IGL_INLINE bool igl::bijective_composite_harmonic_mapping(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedb> & b,
+  const Eigen::MatrixBase<Derivedbc> & bc,
+  Eigen::PlainObjectBase<DerivedU> & U)
+{
+  return bijective_composite_harmonic_mapping(V,F,b,bc,1,200,20,true,U);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedb,
+  typename Derivedbc,
+  typename DerivedU>
+IGL_INLINE bool igl::bijective_composite_harmonic_mapping(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedb> & b,
+  const Eigen::MatrixBase<Derivedbc> & bc,
+  const int min_steps,
+  const int max_steps,
+  const int num_inner_iters,
+  const bool test_for_flips,
+  Eigen::PlainObjectBase<DerivedU> & U)
+{
+  typedef typename Derivedbc::Scalar Scalar;
+  assert(V.cols() == 2 && bc.cols() == 2 && "Input should be 2D");
+  assert(F.cols() == 3 && "F should contain triangles");
+  int tries = 0;
+  int nsteps = min_steps;
+  Derivedbc bc0;
+  slice(V,b,1,bc0);
+
+  // It's difficult to check for flips "robustly" in the sense that the input
+  // mesh might not have positive/consistent sign to begin with. 
+
+  while(nsteps<=max_steps)
+  {
+    U = V;
+    int flipped = 0;
+    int nans = 0;
+    int step = 0;
+    for(;step<=nsteps;step++)
+    {
+      const Scalar t = ((Scalar)step)/((Scalar)nsteps);
+      // linearly interpolate boundary conditions
+      // TODO: replace this with something that guarantees a homotopic "morph"
+      // of the boundary conditions. Something like "Homotopic Morphing of
+      // Planar Curves" [Dym et al. 2015] but also handling multiple connected
+      // components.
+      Derivedbc bct = bc0 + t*(bc - bc0);
+      // Compute dsicrete harmonic map using metric of previous step
+      for(int iter = 0;iter<num_inner_iters;iter++)
+      {
+        //std::cout<<nsteps<<" t: "<<t<<" iter: "<<iter;
+        //igl::matlab::MatlabWorkspace mw;
+        //mw.save(U,"U");
+        //mw.save_index(F,"F");
+        //mw.save_index(b,"b");
+        //mw.save(bct,"bct");
+        //mw.write("numerical.mat");
+        harmonic(DerivedU(U),F,b,bct,1,U);
+        igl::slice(U,b,1,bct);
+        nans = (U.array() != U.array()).count();
+        if(test_for_flips)
+        {
+          Eigen::Matrix<Scalar,Eigen::Dynamic,1> A;
+          doublearea(U,F,A);
+          flipped = (A.array() < 0 ).count();
+          //std::cout<<"  "<<flipped<<"  nan? "<<(U.array() != U.array()).any()<<std::endl;
+          if(flipped == 0 && nans == 0) break;
+        }
+      }
+      if(flipped > 0 || nans>0) break;
+    }
+    if(flipped == 0 && nans == 0)
+    {
+      return step == nsteps+1;
+    }
+    nsteps *= 2;
+  }
+  //std::cout<<"failed to finish in "<<nsteps<<"..."<<std::endl;
+  return false;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::bijective_composite_harmonic_mapping<Eigen::Matrix<double, -1, -1, 1, -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, 1, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&);
+// generated by autoexplicit.sh
+template bool igl::bijective_composite_harmonic_mapping<Eigen::Matrix<double, -1, -1, 1, -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, 1, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, int, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> >&);
+#endif

include/igl/bijective_composite_harmonic_mapping.h

@@ -0,0 +1,79 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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_BIJECTIVE_COMPOSITE_HARMONIC_MAPPING_H
+#define IGL_BIJECTIVE_COMPOSITE_HARMONIC_MAPPING_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl 
+{
+  // Compute a planar mapping of a triangulated polygon (V,F) subjected to
+  // boundary conditions (b,bc). The mapping should be bijective in the sense
+  // that no triangles' areas become negative (this assumes they started
+  // positive). This mapping is computed by "composing" harmonic mappings
+  // between incremental morphs of the boundary conditions. This is a bit like
+  // a discrete version of "Bijective Composite Mean Value Mappings" [Schneider
+  // et al. 2013] but with a discrete harmonic map (cf. harmonic coordinates)
+  // instead of mean value coordinates. This is inspired by "Embedding a
+  // triangular graph within a given boundary" [Xu et al. 2011].
+  //
+  // Inputs:
+  //   V  #V by 2 list of triangle mesh vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  //   b  #b list of boundary indices into V
+  //   bc  #b by 2 list of boundary conditions corresponding to b
+  // Outputs:
+  //   U  #V by 2 list of output mesh vertex locations
+  // Returns true if and only if U contains a successfull bijectie mapping
+  //
+  // 
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Derivedb,
+    typename Derivedbc,
+    typename DerivedU>
+  IGL_INLINE bool bijective_composite_harmonic_mapping(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<Derivedb> & b,
+    const Eigen::MatrixBase<Derivedbc> & bc,
+    Eigen::PlainObjectBase<DerivedU> & U);
+  //
+  // Inputs:
+  //   min_steps  mininum number of steps to take from V(b,:) to bc
+  //   max_steps  mininum number of steps to take from V(b,:) to bc (if
+  //     max_steps == min_steps then no further number of steps will be tried)
+  //   num_inner_iters  number of iterations of harmonic solves to run after
+  //     for each morph step (to try to push flips back in)
+  //   test_for_flips  wether to check if flips occured (and trigger more
+  //     steps). if test_for_flips = false then this function always returns
+  //     true
+  // 
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Derivedb,
+    typename Derivedbc,
+    typename DerivedU>
+  IGL_INLINE bool bijective_composite_harmonic_mapping(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<Derivedb> & b,
+    const Eigen::MatrixBase<Derivedbc> & bc,
+    const int min_steps,
+    const int max_steps,
+    const int num_inner_iters,
+    const bool test_for_flips,
+    Eigen::PlainObjectBase<DerivedU> & U);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bijective_composite_harmonic_mapping.cpp"
+#endif
+#endif

include/igl/boundary_conditions.h

@@ -32,7 +32,11 @@ namespace igl
   //     (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
+  // Returns false if boundary conditions are suspicious:
+  //   P and BE are empty
+  //   bc is empty
+  //   some column of bc doesn't have a 0 (assuming bc has >1 columns)
+  //   some column of bc doesn't have a 1 (assuming bc has >1 columns)
   IGL_INLINE bool boundary_conditions(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & Ele,

include/igl/boundary_facets.cpp

@@ -110,9 +110,9 @@ IGL_INLINE void igl::boundary_facets(
   using namespace std;
   using namespace Eigen;
   // Cop out: use vector of vectors version
-  vector<vector<typename Eigen::PlainObjectBase<DerivedT>::Scalar> > vT;
+  vector<vector<typename DerivedT::Scalar> > vT;
   matrix_to_list(T,vT);
-  vector<vector<typename Eigen::PlainObjectBase<DerivedF>::Scalar> > vF;
+  vector<vector<typename DerivedF::Scalar> > vF;
   boundary_facets(vT,vF);
   list_to_matrix(vF,F);
 }
@@ -128,7 +128,7 @@ Ret igl::boundary_facets(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
 // generated by autoexplicit.sh

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)
   {
@@ -141,6 +148,6 @@ IGL_INLINE void igl::boundary_loop(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::boundary_loop<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/bounding_box.cpp

@@ -10,15 +10,25 @@
 
 template <typename DerivedV, typename DerivedBV, typename DerivedBF>
 IGL_INLINE void igl::bounding_box(
-  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedBV>& BV,
+  Eigen::PlainObjectBase<DerivedBF>& BF)
+{
+  return bounding_box(V,0.,BV,BF);
+}
+
+template <typename DerivedV, typename DerivedBV, typename DerivedBF>
+IGL_INLINE void igl::bounding_box(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const typename DerivedV::Scalar pad,
   Eigen::PlainObjectBase<DerivedBV>& BV,
   Eigen::PlainObjectBase<DerivedBF>& BF)
 {
   using namespace std;
 
   const int dim = V.cols();
-  const auto & minV = V.colwise().minCoeff();
-  const auto & maxV = V.colwise().maxCoeff();
+  const auto & minV = V.colwise().minCoeff().array()-pad;
+  const auto & maxV = V.colwise().maxCoeff().array()+pad;
   // 2^n vertices
   BV.resize((1<<dim),dim);
 
@@ -81,7 +91,13 @@ IGL_INLINE void igl::bounding_box(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-template void igl::bounding_box<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::bounding_box<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::bounding_box<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::bounding_box<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 1, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::bounding_box<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::bounding_box<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, -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> >&);
+template void igl::bounding_box<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/bounding_box.h

@@ -20,7 +20,13 @@ namespace igl
   //   BF  #BF by dim list of simplex facets 
   template <typename DerivedV, typename DerivedBV, typename DerivedBF>
   IGL_INLINE void bounding_box(
-    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::MatrixBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedBV>& BV,
+    Eigen::PlainObjectBase<DerivedBF>& BF);
+  template <typename DerivedV, typename DerivedBV, typename DerivedBF>
+  IGL_INLINE void bounding_box(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const typename DerivedV::Scalar pad,
     Eigen::PlainObjectBase<DerivedBV>& BV,
     Eigen::PlainObjectBase<DerivedBF>& BF);
 }

include/igl/bounding_box_diagonal.cpp

@@ -22,5 +22,5 @@ IGL_INLINE double igl::bounding_box_diagonal(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/canonical_quaternions.cpp

@@ -17,5 +17,5 @@ template <> IGL_INLINE double igl::CANONICAL_VIEW_QUAT<double>(int i, int j)
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/cat.cpp

@@ -252,7 +252,7 @@ IGL_INLINE void igl::cat(const std::vector<std::vector< Mat > > & A, Mat & C)
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::cat<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(int, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&);
 // generated by autoexplicit.sh

include/igl/cat.h

@@ -30,7 +30,7 @@ namespace igl
   // Inputs:
   //   A  first input matrix
   //   B  second input matrix
-  //   dim  dimension along which to concatenate, 0 or 1
+  //   dim  dimension along which to concatenate, 1 or 2
   // Outputs:
   //   C  output matrix
   //   

include/igl/ceil.cpp

@@ -14,7 +14,7 @@ IGL_INLINE void igl::ceil(
   Eigen::PlainObjectBase<DerivedY>& Y)
 {
   using namespace std;
-  //Y = Eigen::PlainObjectBase<DerivedY>::Zero(m,n);
+  //Y = DerivedY::Zero(m,n);
 //#pragma omp parallel for
   //for(int i = 0;i<m;i++)
   //{
@@ -28,6 +28,6 @@ IGL_INLINE void igl::ceil(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit instanciation
+// Explicit template instantiation
 template void igl::ceil<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/circulation.h

@@ -28,7 +28,7 @@ namespace igl
   //     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
   //     e=(j->i)
   //   EI  #E by 2 list of edge flap corners (see above).
-  // Returns list of faces touched by circulation.
+  // Returns list of faces touched by circulation (in cyclically order).
   //   
   IGL_INLINE std::vector<int> circulation(
     const int e,

include/igl/circumradius.cpp

@@ -0,0 +1,26 @@
+// 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 "circumradius.h"
+#include "edge_lengths.h"
+#include "doublearea.h"
+template <
+  typename DerivedV, 
+  typename DerivedF,
+  typename DerivedR>
+IGL_INLINE void igl::circumradius(
+  const Eigen::PlainObjectBase<DerivedV> & V, 
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedR> & R)
+{
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> l;
+  igl::edge_lengths(V,F,l);
+  DerivedR A;
+  igl::doublearea(l,0.,A);
+  // use formula: R=abc/(4*area) to compute the circum radius
+  R = l.col(0).array() * l.col(1).array() * l.col(2).array() / (2.0*A.array());
+}

include/igl/circumradius.h

@@ -0,0 +1,34 @@
+// 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/.
+#ifndef IGL_CIRCUMRADIUS_H
+#define IGL_CIRCUMRADIUS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Compute the circumradius of each triangle in a mesh (V,F)
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  // Outputs:
+  //   R  #F list of circumradii
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedF,
+    typename DerivedR>
+  IGL_INLINE void circumradius(
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedR> & R);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "circumradius.cpp"
+#endif
+#endif

include/igl/collapse_edge.cpp

@@ -176,8 +176,96 @@ IGL_INLINE bool igl::collapse_edge(
   Eigen::MatrixXd & C)
 {
   int e,e1,e2,f1,f2;
+  const auto always_try = [](
+    const Eigen::MatrixXd &                                         ,/*V*/
+    const Eigen::MatrixXi &                                         ,/*F*/
+    const Eigen::MatrixXi &                                         ,/*E*/
+    const Eigen::VectorXi &                                         ,/*EMAP*/
+    const Eigen::MatrixXi &                                         ,/*EF*/
+    const Eigen::MatrixXi &                                         ,/*EI*/
+    const std::set<std::pair<double,int> > &                        ,/*Q*/
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
+    const Eigen::MatrixXd &                                         ,/*C*/
+    const int                                                        /*e*/
+    ) -> bool { return true;};
+  const auto never_care = [](
+    const Eigen::MatrixXd &                                         ,   /*V*/
+    const Eigen::MatrixXi &                                         ,   /*F*/
+    const Eigen::MatrixXi &                                         ,   /*E*/
+    const Eigen::VectorXi &                                         ,/*EMAP*/
+    const Eigen::MatrixXi &                                         ,  /*EF*/
+    const Eigen::MatrixXi &                                         ,  /*EI*/
+    const std::set<std::pair<double,int> > &                        ,   /*Q*/
+    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
+    const Eigen::MatrixXd &                                         ,   /*C*/
+    const int                                                       ,   /*e*/
+    const int                                                       ,  /*e1*/
+    const int                                                       ,  /*e2*/
+    const int                                                       ,  /*f1*/
+    const int                                                       ,  /*f2*/
+    const bool                                                  /*collapsed*/
+    )-> void { };
   return 
-    collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+    collapse_edge(
+      cost_and_placement,always_try,never_care,
+      V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
+}
+
+IGL_INLINE bool igl::collapse_edge(
+  const std::function<void(
+    const int,
+    const Eigen::MatrixXd &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::VectorXi &,
+    const Eigen::MatrixXi &,
+    const Eigen::MatrixXi &,
+    double &,
+    Eigen::RowVectorXd &)> & cost_and_placement,
+  const std::function<bool(
+    const Eigen::MatrixXd &                                         ,/*V*/
+    const Eigen::MatrixXi &                                         ,/*F*/
+    const Eigen::MatrixXi &                                         ,/*E*/
+    const Eigen::VectorXi &                                         ,/*EMAP*/
+    const Eigen::MatrixXi &                                         ,/*EF*/
+    const Eigen::MatrixXi &                                         ,/*EI*/
+    const std::set<std::pair<double,int> > &                        ,/*Q*/
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
+    const Eigen::MatrixXd &                                         ,/*C*/
+    const int                                                        /*e*/
+    )> & pre_collapse,
+  const std::function<void(
+    const Eigen::MatrixXd &                                         ,   /*V*/
+    const Eigen::MatrixXi &                                         ,   /*F*/
+    const Eigen::MatrixXi &                                         ,   /*E*/
+    const Eigen::VectorXi &                                         ,/*EMAP*/
+    const Eigen::MatrixXi &                                         ,  /*EF*/
+    const Eigen::MatrixXi &                                         ,  /*EI*/
+    const std::set<std::pair<double,int> > &                        ,   /*Q*/
+    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
+    const Eigen::MatrixXd &                                         ,   /*C*/
+    const int                                                       ,   /*e*/
+    const int                                                       ,  /*e1*/
+    const int                                                       ,  /*e2*/
+    const int                                                       ,  /*f1*/
+    const int                                                       ,  /*f2*/
+    const bool                                                  /*collapsed*/
+    )> & post_collapse,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E,
+  Eigen::VectorXi & EMAP,
+  Eigen::MatrixXi & EF,
+  Eigen::MatrixXi & EI,
+  std::set<std::pair<double,int> > & Q,
+  std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+  Eigen::MatrixXd & C)
+{
+  int e,e1,e2,f1,f2;
+  return 
+    collapse_edge(
+      cost_and_placement,pre_collapse,post_collapse,
+      V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2);
 }
 
 
@@ -192,6 +280,35 @@ IGL_INLINE bool igl::collapse_edge(
     const Eigen::MatrixXi &,
     double &,
     Eigen::RowVectorXd &)> & cost_and_placement,
+  const std::function<bool(
+    const Eigen::MatrixXd &                                         ,/*V*/
+    const Eigen::MatrixXi &                                         ,/*F*/
+    const Eigen::MatrixXi &                                         ,/*E*/
+    const Eigen::VectorXi &                                         ,/*EMAP*/
+    const Eigen::MatrixXi &                                         ,/*EF*/
+    const Eigen::MatrixXi &                                         ,/*EI*/
+    const std::set<std::pair<double,int> > &                        ,/*Q*/
+    const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
+    const Eigen::MatrixXd &                                         ,/*C*/
+    const int                                                        /*e*/
+    )> & pre_collapse,
+  const std::function<void(
+    const Eigen::MatrixXd &                                         ,   /*V*/
+    const Eigen::MatrixXi &                                         ,   /*F*/
+    const Eigen::MatrixXi &                                         ,   /*E*/
+    const Eigen::VectorXi &                                         ,/*EMAP*/
+    const Eigen::MatrixXi &                                         ,  /*EF*/
+    const Eigen::MatrixXi &                                         ,  /*EI*/
+    const std::set<std::pair<double,int> > &                        ,   /*Q*/
+    const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
+    const Eigen::MatrixXd &                                         ,   /*C*/
+    const int                                                       ,   /*e*/
+    const int                                                       ,  /*e1*/
+    const int                                                       ,  /*e2*/
+    const int                                                       ,  /*f1*/
+    const int                                                       ,  /*f2*/
+    const bool                                                  /*collapsed*/
+    )> & post_collapse,
   Eigen::MatrixXd & V,
   Eigen::MatrixXi & F,
   Eigen::MatrixXi & E,
@@ -225,8 +342,16 @@ IGL_INLINE bool igl::collapse_edge(
   std::vector<int> N  = circulation(e, true,F,E,EMAP,EF,EI);
   std::vector<int> Nd = circulation(e,false,F,E,EMAP,EF,EI);
   N.insert(N.begin(),Nd.begin(),Nd.end());
-  const bool collapsed =
-    collapse_edge(e,C.row(e),V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+  bool collapsed = true;
+  if(pre_collapse(V,F,E,EMAP,EF,EI,Q,Qit,C,e))
+  {
+    collapsed = collapse_edge(e,C.row(e),V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+  }else
+  {
+    // Aborted by pre collapse callback
+    collapsed = false;
+  }
+  post_collapse(V,F,E,EMAP,EF,EI,Q,Qit,C,e,e1,e2,f1,f2,collapsed);
   if(collapsed)
   {
     // Erase the two, other collapsed edges

include/igl/collapse_edge.h

@@ -70,6 +70,9 @@ namespace igl
   //   cost_and_placement  function computing cost of collapsing an edge and 3d
   //     position where it should be placed:
   //     cost_and_placement(V,F,E,EMAP,EF,EI,cost,placement);
+  //     **If the edges is collapsed** then this function will be called on all
+  //     edges of all faces previously incident on the endpoints of the
+  //     collapsed edge.
   //   Q  queue containing pairs of costs and edge indices
   //   Qit  list of iterators so that Qit[e] --> iterator of edge e in Q
   //   C  #E by dim list of stored placements
@@ -93,6 +96,64 @@ namespace igl
     std::set<std::pair<double,int> > & Q,
     std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
     Eigen::MatrixXd & C);
+  // Inputs:
+  //   pre_collapse  callback called with index of edge whose collapse is about
+  //     to be attempted. This function should return whether to **proceed**
+  //     with the collapse: returning true means "yes, try to collapse",
+  //     returning false means "No, consider this edge 'uncollapsable', behave
+  //     as if collapse_edge(e) returned false.
+  //   post_collapse  callback called with index of edge whose collapse was
+  //     just attempted and a flag revealing whether this was successful.
+  IGL_INLINE bool collapse_edge(
+    const std::function<void(
+      const int,
+      const Eigen::MatrixXd &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::VectorXi &,
+      const Eigen::MatrixXi &,
+      const Eigen::MatrixXi &,
+      double &,
+      Eigen::RowVectorXd &)> & cost_and_placement,
+    const std::function<bool(
+      const Eigen::MatrixXd &                                         ,/*V*/
+      const Eigen::MatrixXi &                                         ,/*F*/
+      const Eigen::MatrixXi &                                         ,/*E*/
+      const Eigen::VectorXi &                                         ,/*EMAP*/
+      const Eigen::MatrixXi &                                         ,/*EF*/
+      const Eigen::MatrixXi &                                         ,/*EI*/
+      const std::set<std::pair<double,int> > &                        ,/*Q*/
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
+      const Eigen::MatrixXd &                                         ,/*C*/
+      const int                                                        /*e*/
+      )> & pre_collapse,
+    const std::function<void(
+      const Eigen::MatrixXd &                                         ,   /*V*/
+      const Eigen::MatrixXi &                                         ,   /*F*/
+      const Eigen::MatrixXi &                                         ,   /*E*/
+      const Eigen::VectorXi &                                         ,/*EMAP*/
+      const Eigen::MatrixXi &                                         ,  /*EF*/
+      const Eigen::MatrixXi &                                         ,  /*EI*/
+      const std::set<std::pair<double,int> > &                        ,   /*Q*/
+      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
+      const Eigen::MatrixXd &                                         ,   /*C*/
+      const int                                                       ,   /*e*/
+      const int                                                       ,  /*e1*/
+      const int                                                       ,  /*e2*/
+      const int                                                       ,  /*f1*/
+      const int                                                       ,  /*f2*/
+      const bool                                                  /*collapsed*/
+      )> & post_collapse,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E,
+    Eigen::VectorXi & EMAP,
+    Eigen::MatrixXi & EF,
+    Eigen::MatrixXi & EI,
+    std::set<std::pair<double,int> > & Q,
+    std::vector<std::set<std::pair<double,int> >::iterator > & Qit,
+    Eigen::MatrixXd & C);
+
   IGL_INLINE bool collapse_edge(
     const std::function<void(
       const int,
@@ -104,6 +165,35 @@ namespace igl
       const Eigen::MatrixXi &,
       double &,
       Eigen::RowVectorXd &)> & cost_and_placement,
+    const std::function<bool(
+      const Eigen::MatrixXd &                                         ,/*V*/
+      const Eigen::MatrixXi &                                         ,/*F*/
+      const Eigen::MatrixXi &                                         ,/*E*/
+      const Eigen::VectorXi &                                         ,/*EMAP*/
+      const Eigen::MatrixXi &                                         ,/*EF*/
+      const Eigen::MatrixXi &                                         ,/*EI*/
+      const std::set<std::pair<double,int> > &                        ,/*Q*/
+      const std::vector<std::set<std::pair<double,int> >::iterator > &,/*Qit*/
+      const Eigen::MatrixXd &                                         ,/*C*/
+      const int                                                        /*e*/
+      )> & pre_collapse,
+    const std::function<void(
+      const Eigen::MatrixXd &                                         ,   /*V*/
+      const Eigen::MatrixXi &                                         ,   /*F*/
+      const Eigen::MatrixXi &                                         ,   /*E*/
+      const Eigen::VectorXi &                                         ,/*EMAP*/
+      const Eigen::MatrixXi &                                         ,  /*EF*/
+      const Eigen::MatrixXi &                                         ,  /*EI*/
+      const std::set<std::pair<double,int> > &                        ,   /*Q*/
+      const std::vector<std::set<std::pair<double,int> >::iterator > &, /*Qit*/
+      const Eigen::MatrixXd &                                         ,   /*C*/
+      const int                                                       ,   /*e*/
+      const int                                                       ,  /*e1*/
+      const int                                                       ,  /*e2*/
+      const int                                                       ,  /*f1*/
+      const int                                                       ,  /*f2*/
+      const bool                                                  /*collapsed*/
+      )> & post_collapse,
     Eigen::MatrixXd & V,
     Eigen::MatrixXi & F,
     Eigen::MatrixXi & E,

include/igl/collapse_small_triangles.cpp

@@ -31,7 +31,7 @@ void igl::collapse_small_triangles(
   MatrixXd l;
   edge_lengths(V,F,l);
   VectorXd dblA;
-  doublearea(l,dblA);
+  doublearea(l,0.,dblA);
 
   // Minimum area tolerance
   const double min_dblarea = 2.0*eps*bbd*bbd;
@@ -91,7 +91,7 @@ void igl::collapse_small_triangles(
     }
   }
 
-  FF.resize(rF.rows(),rF.cols());
+  FF.resizeLike(rF);
   int num_face_collapses=0;
   // Only keep uncollapsed faces
   {

include/igl/colon.cpp

@@ -40,7 +40,7 @@ IGL_INLINE Eigen::Matrix<T,Eigen::Dynamic,1> igl::colon(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template Eigen::Matrix<int,-1,1,0,-1,1> igl::colon<int,int,int>(int, int);
 template Eigen::Matrix<int,-1,1,0,-1,1> igl::colon<int,int,long>(int,long);

include/igl/columnize.cpp

@@ -55,7 +55,7 @@ IGL_INLINE void igl::columnize(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::columnize<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&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::columnize<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&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::columnize<Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, int, int, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);

include/igl/comb_cross_field.cpp

@@ -142,7 +142,7 @@ IGL_INLINE void igl::comb_cross_field(const Eigen::PlainObjectBase<DerivedV> &V,
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::comb_cross_field<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::comb_cross_field<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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

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.));
     }
@@ -75,7 +75,7 @@ IGL_INLINE void igl::comb_frame_field(const Eigen::PlainObjectBase<DerivedV> &V,
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::comb_frame_field<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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::comb_frame_field<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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/comb_line_field.cpp

@@ -128,5 +128,5 @@ IGL_INLINE void igl::comb_line_field(const Eigen::PlainObjectBase<DerivedV> &V,
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/combine.cpp

@@ -0,0 +1,57 @@
+// 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 "combine.h"
+#include <cassert>
+template <
+  typename DerivedVV, 
+  typename DerivedFF, 
+  typename DerivedV, 
+  typename DerivedF>
+IGL_INLINE void igl::combine(
+  const std::vector<DerivedVV> & VV,
+  const std::vector<DerivedFF> & FF,
+  Eigen::PlainObjectBase<DerivedV> & V,
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  assert(VV.size() == FF.size() && 
+    "Lists of verex lists and face lists should be same size");
+  // Dimension of vertex positions
+  const int dim = VV.size() > 0 ? VV[0].cols() : 0;
+  // Simplex/element size
+  const int ss = FF.size() > 0 ? FF[0].cols() : 0;
+  int n = 0;
+  int m = 0;
+  for(int i = 0;i<VV.size();i++)
+  {
+    const auto & Vi = VV[i];
+    const auto & Fi = FF[i];
+    n+=Vi.rows();
+    assert(dim == Vi.cols() && "All vertex lists should have same #columns");
+    m+=Fi.rows();
+    assert(ss == Fi.cols() && "All face lists should have same #columns");
+  }
+  V.resize(n,dim);
+  F.resize(m,ss);
+  {
+    int kv = 0;
+    int kf = 0;
+    for(int i = 0;i<VV.size();i++)
+    {
+      const auto & Vi = VV[i];
+      const int ni = Vi.rows();
+      const auto & Fi = FF[i];
+      const int mi = Fi.rows();
+      F.block(kf,0,mi,ss) = Fi.array()+kv;
+      kf+=mi;
+      V.block(kv,0,ni,dim) = Vi;
+      kv+=ni;
+    }
+    assert(kv == V.rows());
+    assert(kf == F.rows());
+  }
+}

include/igl/combine.h

@@ -0,0 +1,48 @@
+// 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/.
+#ifndef IGL_COMBINE_H
+#define IGL_COMBINE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // Combine k meshes into a single >=k connected component mesh with a single
+  // vertex list and face list. Similar to Maya's Combine operation.
+  //
+  // Inputs:
+  //   VV  k-long list of lists of mesh vertex positions
+  //   FF  k-long list of lists of mesh face indices so that FF[i] indexes
+  //     VV[i]
+  // Outputs:
+  //   V   VV[0].rows()+...+VV[k-1].rows() by VV[0].cols() list of mesh
+  //     vertex positions
+  //   F   FF[0].rows()+...+FF[k-1].rows() by FF[0].cols() list of mesh faces
+  //     indices into V
+  // Example:
+  //   // Suppose you have mesh A (VA,FA) and mesh B (VB,FB)
+  //   igl::combine<Eigen::MatrixXd,Eigen::MatrixXi>({VA,VB},{FA,FB},V,F);
+  //
+  template <
+    typename DerivedVV, 
+    typename DerivedFF, 
+    typename DerivedV, 
+    typename DerivedF>
+  IGL_INLINE void combine(
+    const std::vector<DerivedVV> & VV,
+    const std::vector<DerivedFF> & FF,
+    Eigen::PlainObjectBase<DerivedV> & V,
+    Eigen::PlainObjectBase<DerivedF> & F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "combine.cpp"
+#endif
+#endif

include/igl/components.cpp

@@ -24,6 +24,7 @@ IGL_INLINE void igl::components(
   C.resize(n,1);
   typename DerivedC::Scalar id = 0;
   vector<typename Derivedcounts::Scalar> vcounts;
+  // breadth first search
   for(int k=0; k<A.outerSize(); ++k)
   {
     if(seen(k))
@@ -48,7 +49,7 @@ IGL_INLINE void igl::components(
       for(typename SparseMatrix<AScalar>::InnerIterator it (A,f); it; ++it)
       {
         const int g = it.index();
-        if(!seen(g))
+        if(!seen(g) && it.value())
         {
           Q.push(g);
         }
@@ -77,7 +78,7 @@ IGL_INLINE void igl::components(
 
 template <typename DerivedF, typename DerivedC>
 IGL_INLINE void igl::components(
-  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedF> & F,
   Eigen::PlainObjectBase<DerivedC> & C)
 {
   Eigen::SparseMatrix<typename DerivedC::Scalar> A;
@@ -86,9 +87,12 @@ IGL_INLINE void igl::components(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-template void igl::components<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::components<bool, Eigen::Array<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<bool, 0, int> const&, Eigen::PlainObjectBase<Eigen::Array<int, -1, 1, 0, -1, 1> >&);
+template void igl::components<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::components<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::components<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<int, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::components<double, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::components<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/components.h

@@ -39,7 +39,7 @@ namespace igl
   //   C  max(F) list of component ids
   template <typename DerivedF, typename DerivedC>
   IGL_INLINE void components(
-    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedF> & F,
     Eigen::PlainObjectBase<DerivedC> & C);
 
 }

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.));
 
@@ -78,7 +78,7 @@ IGL_INLINE void igl::compute_frame_field_bisectors(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::compute_frame_field_bisectors<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::compute_frame_field_bisectors<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<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/conjugate_frame_fields.cpp

@@ -428,6 +428,6 @@ IGL_INLINE typename DerivedV::Scalar igl::conjugate_frame_fields(const igl::Conj
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar igl::conjugate_frame_fields<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(igl::ConjugateFFSolverData<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, bool);
 #endif

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -72,13 +72,10 @@ namespace igl
           typedef CGAL::Segment_2<Kernel>  Segment_2; 
           typedef CGAL::Triangle_2<Kernel> Triangle_2; 
           // 2D Constrained Delaunay Triangulation types
-          typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
-          typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
-          typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
+          //typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
+          //typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
+          //typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
           typedef CGAL::Exact_intersections_tag Itag;
-          typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
-            CDT_2;
-          typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
           // Axis-align boxes for all-pairs self-intersection detection
           typedef std::vector<Triangle_3> Triangles;
           typedef typename Triangles::iterator TrianglesIterator;
@@ -206,18 +203,6 @@ namespace igl
           //   b  box containing a triangle
           inline void box_intersect(const Box& a, const Box& b);
           inline void process_intersecting_boxes();
-        private:
-          // Compute 2D delaunay triangulation of a given 3d triangle and a list of
-          // intersection objects (points,segments,triangles). CGAL uses an affine
-          // projection rather than an isometric projection, so we're not
-          // guaranteed that the 2D delaunay triangulation here will be a delaunay
-          // triangulation in 3D.
-          //
-          // Inputs:
-          //   A  triangle in 3D
-          //   A_objects_3  updated list of intersection objects for A
-          // Outputs:
-          //   cdt  Contrained delaunay triangulation in projected 2D plane
         public:
           // Getters:
           //const IndexList& get_lIF() const{ return lIF;}

include/igl/copyleft/cgal/assign.cpp

@@ -0,0 +1,54 @@
+// 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 "assign.h"
+#include "assign_scalar.h"
+
+template <typename DerivedC, typename DerivedD>
+IGL_INLINE void igl::copyleft::cgal::assign(
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  D.resizeLike(C);
+  for(int i = 0;i<C.rows();i++)
+  {
+    for(int j = 0;j<C.cols();j++)
+    {
+      assign_scalar(C(i,j),D(i,j));
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -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<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -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<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::assign<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<double, 8, 3, 0, 8, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 8, 3, 0, 8, 3> >&);
+#endif

include/igl/copyleft/cgal/assign.h

@@ -0,0 +1,31 @@
+// 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/.
+#ifndef IGL_COPYLEFT_CGAL_ASSIGN_H
+#define IGL_COPYLEFT_CGAL_ASSIGN_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      template <typename DerivedC, typename DerivedD>
+      IGL_INLINE void assign(
+        const Eigen::PlainObjectBase<DerivedC> & C,
+        Eigen::PlainObjectBase<DerivedD> & D);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "assign.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/assign_scalar.cpp

@@ -8,18 +8,18 @@
 #include "assign_scalar.h"
 
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(
-  const typename CGAL::Epeck::FT & cgal,
+  const CGAL::Epeck::FT & cgal,
   CGAL::Epeck::FT & d)
 {
   d = cgal;
 }
 
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(
-  const typename CGAL::Epeck::FT & _cgal,
+  const CGAL::Epeck::FT & _cgal,
   double & d)
 {
   // FORCE evaluation of the exact type otherwise interval might be huge.
-  const typename CGAL::Epeck::FT cgal = _cgal.exact();
+  const CGAL::Epeck::FT cgal = _cgal.exact();
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {
@@ -30,11 +30,11 @@ IGL_INLINE void igl::copyleft::cgal::assign_scalar(
 }
 
 IGL_INLINE void igl::copyleft::cgal::assign_scalar(
-  const typename CGAL::Epeck::FT & _cgal,
+  const CGAL::Epeck::FT & _cgal,
   float& d)
 {
   // FORCE evaluation of the exact type otherwise interval might be huge.
-  const typename CGAL::Epeck::FT cgal = _cgal.exact();
+  const CGAL::Epeck::FT cgal = _cgal.exact();
   const auto interval = CGAL::to_interval(cgal);
   d = interval.first;
   do {

include/igl/copyleft/cgal/assign_scalar.h

@@ -20,13 +20,13 @@ namespace igl
       // Outputs:
       //   d  output scalar
       IGL_INLINE void assign_scalar(
-        const typename CGAL::Epeck::FT & cgal,
+        const CGAL::Epeck::FT & cgal,
         CGAL::Epeck::FT & d);
       IGL_INLINE void assign_scalar(
-        const typename CGAL::Epeck::FT & cgal,
+        const CGAL::Epeck::FT & cgal,
         double & d);
       IGL_INLINE void assign_scalar(
-        const typename CGAL::Epeck::FT & cgal,
+        const CGAL::Epeck::FT & cgal,
         float& d);
       IGL_INLINE void assign_scalar(
         const double & c,

include/igl/copyleft/cgal/barycenter.cpp

@@ -0,0 +1,15 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 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 "../../barycenter.h"
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#ifdef IGL_STATIC_LIBRARY
+#undef IGL_STATIC_LIBRARY
+#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> >&);
+#endif

include/igl/copyleft/cgal/closest_facet.cpp

@@ -15,7 +15,7 @@
 #include "order_facets_around_edge.h"
 #include "submesh_aabb_tree.h"
 #include "../../vertex_triangle_adjacency.h"
-#include "../../writePLY.h"
+//#include "../../writePLY.h"
 
 template<
   typename DerivedV,
@@ -402,7 +402,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
       for (size_t i=0; i<adj_faces.size(); i++) {
         tmp_faces.row(i) = F.row(adj_faces[i]);
       }
-      igl::writePLY("debug.ply", tmp_vertices, tmp_faces, false);
+      //igl::writePLY("debug.ply", tmp_vertices, tmp_faces, false);
       throw std::runtime_error("Invalid vertex neighborhood");
     }
     const auto itr = v2f.equal_range(d);

include/igl/copyleft/cgal/component_inside_component.cpp

@@ -7,19 +7,18 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "component_inside_component.h"
 
-#include <cassert>
-#include <list>
-#include <limits>
-#include <vector>
+#include "order_facets_around_edge.h"
+#include "points_inside_component.h"
 
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
 #include <CGAL/AABB_triangle_primitive.h>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
-#include "order_facets_around_edge.h"
-#include "assign_scalar.h"
-#include "points_inside_component.h"
+#include <cassert>
+#include <list>
+#include <limits>
+#include <vector>
 
 template <typename DerivedV, typename DerivedF, typename DerivedI>
 IGL_INLINE bool igl::copyleft::cgal::component_inside_component(
@@ -61,7 +60,7 @@ IGL_INLINE bool igl::copyleft::cgal::component_inside_component(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template bool igl::copyleft::cgal::component_inside_component<
 Eigen::Matrix<double, -1, -1, 0, -1, -1>,
 Eigen::Matrix<   int, -1, -1, 0, -1, -1>,

include/igl/copyleft/cgal/extract_feature.cpp

@@ -0,0 +1,123 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "extract_feature.h"
+#include <igl/unique_edge_map.h>
+#include <CGAL/Kernel/global_functions.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedE >
+IGL_INLINE void igl::copyleft::cgal::extract_feature(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const double tol,
+    Eigen::PlainObjectBase<DerivedE>& feature_edges) {
+
+  using IndexType = typename DerivedE::Scalar;
+  DerivedE E, uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<IndexType> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+  igl::copyleft::cgal::extract_feature(V, F, tol, E, uE, uE2E, feature_edges);
+}
+
+template<
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedE >
+IGL_INLINE void igl::copyleft::cgal::extract_feature(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const double tol,
+    const Eigen::PlainObjectBase<DerivedE>& E,
+    const Eigen::PlainObjectBase<DerivedE>& uE,
+    const std::vector<std::vector<typename DerivedE::Scalar> >& uE2E,
+    Eigen::PlainObjectBase<DerivedE>& feature_edges) {
+
+  assert(V.cols() == 3);
+  assert(F.cols() == 3);
+  using Scalar = typename DerivedV::Scalar;
+  using IndexType = typename DerivedE::Scalar;
+  using Vertex = Eigen::Matrix<Scalar, 3, 1>;
+  using Kernel = typename CGAL::Exact_predicates_exact_constructions_kernel;
+  using Point = typename Kernel::Point_3;
+
+  const size_t num_unique_edges = uE.rows();
+  const size_t num_faces = F.rows();
+  // NOTE: CGAL's definition of dihedral angle measures the angle between two
+  // facets instead of facet normals.
+  const double cos_tol = cos(M_PI - tol);
+  std::vector<size_t> result; // Indices into uE
+
+  auto is_non_manifold = [&uE2E](size_t ei) -> bool {
+    return uE2E[ei].size() > 2;
+  };
+
+  auto is_boundary = [&uE2E](size_t ei) -> bool {
+    return uE2E[ei].size() == 1;
+  };
+
+  auto opposite_vertex = [&uE, &F](size_t ei, size_t fi) -> IndexType {
+    const size_t v0 = uE(ei, 0);
+    const size_t v1 = uE(ei, 1);
+    for (size_t i=0; i<3; i++) {
+      const size_t v = F(fi, i);
+      if (v != v0 && v != v1) { return v; }
+    }
+    throw "Input face must be topologically degenerate!";
+  };
+
+  auto is_feature = [&V, &F, &uE, &uE2E, &opposite_vertex, num_faces](
+      size_t ei, double cos_tol) -> bool {
+    auto adj_faces = uE2E[ei];
+    assert(adj_faces.size() == 2);
+    const Vertex v0 = V.row(uE(ei, 0));
+    const Vertex v1 = V.row(uE(ei, 1));
+    const Vertex v2 = V.row(opposite_vertex(ei, adj_faces[0] % num_faces));
+    const Vertex v3 = V.row(opposite_vertex(ei, adj_faces[1] % num_faces));
+    const Point p0(v0[0], v0[1], v0[2]);
+    const Point p1(v1[0], v1[1], v1[2]);
+    const Point p2(v2[0], v2[1], v2[2]);
+    const Point p3(v3[0], v3[1], v3[2]);
+    const auto ori = CGAL::orientation(p0, p1, p2, p3);
+    switch (ori) {
+      case CGAL::POSITIVE:
+        return CGAL::compare_dihedral_angle(p0, p1, p2, p3, cos_tol) ==
+          CGAL::SMALLER;
+      case CGAL::NEGATIVE:
+        return CGAL::compare_dihedral_angle(p0, p1, p3, p2, cos_tol) ==
+          CGAL::SMALLER;
+      case CGAL::COPLANAR:
+        if (!CGAL::collinear(p0, p1, p2) && !CGAL::collinear(p0, p1, p3)) {
+          return CGAL::compare_dihedral_angle(p0, p1, p2, p3, cos_tol) ==
+            CGAL::SMALLER;
+        } else {
+          throw "Dihedral angle (and feature edge) is not well defined for"
+              " degenerate triangles!";
+        }
+      default:
+        throw "Unknown CGAL orientation";
+    }
+  };
+
+  for (size_t i=0; i<num_unique_edges; i++) {
+    if (is_boundary(i) || is_non_manifold(i) || is_feature(i, cos_tol)) {
+      result.push_back(i);
+    }
+  }
+
+  const size_t num_feature_edges = result.size();
+  feature_edges.resize(num_feature_edges, 2);
+  for (size_t i=0; i<num_feature_edges; i++) {
+    feature_edges.row(i) = uE.row(result[i]);
+  }
+}

include/igl/copyleft/cgal/extract_feature.h

@@ -0,0 +1,90 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Qingnan Zhou <qnzhou@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_COPYLEFT_CGAL_EXTRACT_FEATURE_H
+#define IGL_COPYLEFT_CGAL_EXTRACT_FEATURE_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Extract feature edges based on dihedral angle.
+      // Here, dihedral angle is defined as the angle between surface
+      // __normals__ as described in
+      // http://mathworld.wolfram.com/DihedralAngle.html
+      //
+      // Non-manifold and boundary edges are automatically considered as
+      // features.
+      //
+      // Inputs:
+      //   V   #V by 3 array of vertices.
+      //   F   #F by 3 array of faces.
+      //   tol Edges with dihedral angle larger than this are considered
+      //       as features.  Angle is measured in radian.
+      //
+      // Output:
+      //   feature_edges: #E by 2 array of edges.  Each edge satisfies at
+      //      least one of the following criteria:
+      //
+      //      * Edge has dihedral angle larger than tol.
+      //      * Edge is boundary.
+      //      * Edge is non-manifold (i.e. it has more than 2 adjacent
+      //        faces).
+      template <
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedE>
+      IGL_INLINE void extract_feature(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const double tol,
+            Eigen::PlainObjectBase<DerivedE>& feature_edges);
+
+      // Inputs:
+      //   V    #V by 3 array of vertices.
+      //   F    #F by 3 array of faces.
+      //   tol  Edges with dihedral angle larger than this are considered
+      //        as features.  Angle is measured in radian.
+      //   E    #E by 2 array of directed edges.
+      //   uE   #uE by 2 array of undirected edges.
+      //   uE2E #uE list of lists mapping undirected edges to all corresponding
+      //        directed edges.
+      //
+      // Output:
+      //   feature_edges: #E by 2 array of edges.  Each edge satisfies at
+      //      least one of the following criteria:
+      //
+      //      * Edge has dihedral angle larger than tol.
+      //      * Edge is boundary.
+      //      * Edge is non-manifold (i.e. it has more than 2 adjacent
+      //        faces).
+      template <
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedE>
+      IGL_INLINE void extract_feature(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const double tol,
+            const Eigen::PlainObjectBase<DerivedE>& E,
+            const Eigen::PlainObjectBase<DerivedE>& uE,
+            const std::vector<std::vector<typename DerivedE::Scalar> >& uE2E,
+            Eigen::PlainObjectBase<DerivedE>& feature_edges);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "extract_feature.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/half_space_box.cpp

@@ -111,7 +111,7 @@ IGL_INLINE void igl::copyleft::cgal::half_space_box(
   return half_space_box(P,V,BV,BF);
 }
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template void igl::copyleft::cgal::half_space_box<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 1, 4, 1, 1, 4>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 1, 4, 1, 1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>&, Eigen::Matrix<int, 12, 3, 0, 12, 3>&);
 // generated by autoexplicit.sh

include/igl/copyleft/cgal/insert_into_cdt.cpp

@@ -0,0 +1,69 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// 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 "insert_into_cdt.h"
+#include <CGAL/Point_3.h>
+#include <CGAL/Segment_3.h>
+#include <CGAL/Triangle_3.h>
+
+template <typename Kernel>
+IGL_INLINE void igl::copyleft::cgal::insert_into_cdt(
+  const CGAL::Object & obj,
+  const CGAL::Plane_3<Kernel> & P,
+  CGAL::Constrained_triangulation_plus_2<
+    CGAL::Constrained_Delaunay_triangulation_2<
+      Kernel,
+      CGAL::Triangulation_data_structure_2<
+        CGAL::Triangulation_vertex_base_2<Kernel>,
+        CGAL::Constrained_triangulation_face_base_2< Kernel>
+      >,
+      CGAL::Exact_intersections_tag
+    > 
+  > 
+  & cdt)
+{
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Segment_3<Kernel>  Segment_3; 
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef CGAL::Plane_3<Kernel>    Plane_3;
+
+  if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj)) 
+  {
+    // Add segment constraint
+    cdt.insert_constraint( P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
+  }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj)) 
+  {
+    // Add point
+    cdt.insert(P.to_2d(*ipoint));
+  } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj)) 
+  {
+    // Add 3 segment constraints
+    cdt.insert_constraint( P.to_2d(itri->vertex(0)),P.to_2d(itri->vertex(1)));
+    cdt.insert_constraint( P.to_2d(itri->vertex(1)),P.to_2d(itri->vertex(2)));
+    cdt.insert_constraint( P.to_2d(itri->vertex(2)),P.to_2d(itri->vertex(0)));
+  } else if(const std::vector<Point_3 > *polyp = 
+      CGAL::object_cast< std::vector<Point_3 > >(&obj)) 
+  {
+    const std::vector<Point_3 > & poly = *polyp;
+    const size_t m = poly.size();
+    assert(m>=2);
+    for(size_t p = 0;p<m;p++)
+    {
+      const size_t np = (p+1)%m;
+      cdt.insert_constraint(P.to_2d(poly[p]),P.to_2d(poly[np]));
+    }
+  }else {
+    throw std::runtime_error("Unknown intersection object!");
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+template void igl::copyleft::cgal::insert_into_cdt<CGAL::Epick>(CGAL::Object const&, CGAL::Plane_3<CGAL::Epick> const&, CGAL::Constrained_triangulation_plus_2<CGAL::Constrained_Delaunay_triangulation_2<CGAL::Epick, CGAL::Triangulation_data_structure_2<CGAL::Triangulation_vertex_base_2<CGAL::Epick, CGAL::Triangulation_ds_vertex_base_2<void> >, CGAL::Constrained_triangulation_face_base_2<CGAL::Epick, CGAL::Triangulation_face_base_2<CGAL::Epick, CGAL::Triangulation_ds_face_base_2<void> > > >, CGAL::Exact_intersections_tag> >&);
+template void igl::copyleft::cgal::insert_into_cdt<CGAL::Epeck>(CGAL::Object const&, CGAL::Plane_3<CGAL::Epeck> const&, CGAL::Constrained_triangulation_plus_2<CGAL::Constrained_Delaunay_triangulation_2<CGAL::Epeck, CGAL::Triangulation_data_structure_2<CGAL::Triangulation_vertex_base_2<CGAL::Epeck, CGAL::Triangulation_ds_vertex_base_2<void> >, CGAL::Constrained_triangulation_face_base_2<CGAL::Epeck, CGAL::Triangulation_face_base_2<CGAL::Epeck, CGAL::Triangulation_ds_face_base_2<void> > > >, CGAL::Exact_intersections_tag> >&);
+#endif

include/igl/copyleft/cgal/insert_into_cdt.h

@@ -0,0 +1,59 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_COPYLEFT_CGAL_INSERT_INTO_CDT_H
+#define IGL_COPYLEFT_CGAL_INSERT_INTO_CDT_H
+#include "../../igl_inline.h"
+
+#include <CGAL/Plane_3.h>
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Constrained_triangulation_plus_2.h>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given a current 2D constrained Delaunay triangulation (cdt), insert a
+      // 3D "object" (e.g., resulting from intersecting two triangles) into the
+      // cdt, by projecting it via the given plane (P) and adding appropriate
+      // constraints.
+      //
+      // Inputs:
+      //   obj  CGAL::Object representing a vertex, segment, or (convex)
+      //     polygon. All vertices should lie on the plane P. If not, then this
+      //     adds the _projection_ of this object to the cdt and that might not
+      //     be what you wanted to do.
+      //   P  plane obj lies on and upon which the cdt is being performed
+      //   cdt  current CDT, see output
+      // Outputs:
+      //   cdt  CDT updated to contain constraints for the given object
+      // 
+      template <typename Kernel>
+      IGL_INLINE void insert_into_cdt(
+        const CGAL::Object & obj,
+        const CGAL::Plane_3<Kernel> & P,
+        CGAL::Constrained_triangulation_plus_2<
+          CGAL::Constrained_Delaunay_triangulation_2<
+            Kernel,
+            CGAL::Triangulation_data_structure_2<
+              CGAL::Triangulation_vertex_base_2<Kernel>,
+              CGAL::Constrained_triangulation_face_base_2< Kernel>
+            >,
+            CGAL::Exact_intersections_tag
+          > 
+        > 
+        & cdt);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "insert_into_cdt.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/intersect_other.cpp

@@ -96,9 +96,6 @@ namespace igl
         typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTAB_2;
         typedef CGAL::Triangulation_data_structure_2<TVB_2,CTAB_2> TDS_2;
         typedef CGAL::Exact_intersections_tag Itag;
-        typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
-          CDT_2;
-        typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
         // Axis-align boxes for all-pairs self-intersection detection
         typedef std::vector<Triangle_3> Triangles;
         typedef typename Triangles::iterator TrianglesIterator;
@@ -285,5 +282,5 @@ IGL_INLINE bool igl::copyleft::cgal::intersect_other(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/copyleft/cgal/intersect_with_half_space.cpp

@@ -76,7 +76,7 @@ IGL_INLINE bool igl::copyleft::cgal::intersect_with_half_space(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 template bool igl::copyleft::cgal::intersect_with_half_space<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 1, 4, 1, 1, 4>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 1, 4, 1, 1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh

include/igl/copyleft/cgal/mesh_boolean.cpp

@@ -8,7 +8,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 //
 #include "mesh_boolean.h"
-#include "assign_scalar.h"
+#include "assign.h"
 #include "extract_cells.h"
 #include "mesh_boolean_type_to_funcs.h"
 #include "propagate_winding_numbers.h"
@@ -113,7 +113,8 @@ IGL_INLINE bool igl::copyleft::cgal::mesh_boolean(
   {
     for(int d = 0;d<3;d++) VV(VA.rows()+b,d) = VB(b,d);
   }
-  FF << FA, FB.array() + VA.rows();
+  FF.block(0, 0, FA.rows(), 3) = FA;
+  FF.block(FA.rows(), 0, FB.rows(), 3) = FB.array() + VA.rows();
   return mesh_boolean(VV,FF,sizes,wind_num_op,keep,VC,FC,J);
 }
 
@@ -393,14 +394,8 @@ IGL_INLINE bool igl::copyleft::cgal::mesh_boolean(
     }
 #endif
 
-    MatrixX3S Vs(V.rows(), V.cols());
-    for (size_t i=0; i<(size_t)V.rows(); i++)
-    {
-      for (size_t j=0; j<(size_t)V.cols(); j++)
-      {
-        igl::copyleft::cgal::assign_scalar(V(i,j), Vs(i,j));
-      }
-    }
+    MatrixX3S Vs;
+    assign(V,Vs);
     Eigen::VectorXi newIM;
     igl::remove_unreferenced(Vs,G,VC,FC,newIM);
   }
@@ -431,10 +426,8 @@ IGL_INLINE bool igl::copyleft::cgal::mesh_boolean(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-// generated by autoexplicit.sh
+// Explicit template instantiation
 template bool igl::copyleft::cgal::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, 12, 3, 0, 12, 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> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, 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::copyleft::cgal::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 4, 0, -1, 4>, 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, 12, 3, 0, 12, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 4, 0, -1, 4> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template bool igl::copyleft::cgal::mesh_boolean<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<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> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template bool igl::copyleft::cgal::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<long, -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, -1, 0, -1, -1> > 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> > const&, igl::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
@@ -445,9 +438,4 @@ template bool igl::copyleft::cgal::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_n
 template bool igl::copyleft::cgal::mesh_boolean<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3>, Eigen::Matrix<int, 12, 3, 0, 12, 3>, 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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, 8, 3, 0, 8, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, 12, 3, 0, 12, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType 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 bool igl::copyleft::cgal::mesh_boolean<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::basic_string<char, std::char_traits<char>, std::allocator<char> > 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 bool igl::copyleft::cgal::mesh_boolean<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::MeshBooleanType 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> >&);
-#undef IGL_STATIC_LIBRARY
-#include "../../remove_unreferenced.cpp"
-template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-#include "../../slice.cpp"
-template void igl::slice<Eigen::PlainObjectBase<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> > >(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> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/copyleft/cgal/mesh_boolean_type_to_funcs.cpp

@@ -35,5 +35,5 @@ IGL_INLINE void igl::copyleft::cgal::mesh_boolean_type_to_funcs(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #endif

include/igl/copyleft/cgal/mesh_to_cgal_triangle_list.cpp

@@ -6,7 +6,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 "mesh_to_cgal_triangle_list.h"
-#include "assign_scalar.h"
+#include "assign.h"
 
 #include <cassert>
 
@@ -30,14 +30,7 @@ IGL_INLINE void igl::copyleft::cgal::mesh_to_cgal_triangle_list(
     DerivedV::RowsAtCompileTime,
     DerivedV::ColsAtCompileTime> 
     KV(V.rows(),V.cols());
-  // Just use f'ing for loops. What if V and KV don't use same ordering?
-  for(int i = 0;i<V.rows();i++)
-  {
-    for(int j = 0;j<V.cols();j++)
-    {
-      assign_scalar(V(i,j),KV(i,j));
-    }
-  }
+  assign(V,KV);
   // Must be triangles
   assert(F.cols() == 3);
   T.reserve(F.rows());
@@ -53,7 +46,7 @@ IGL_INLINE void igl::copyleft::cgal::mesh_to_cgal_triangle_list(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // 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

include/igl/copyleft/cgal/mesh_to_polyhedron.cpp

@@ -47,7 +47,7 @@ IGL_INLINE bool igl::copyleft::cgal::mesh_to_polyhedron(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Polyhedron_items_with_id_3.h>
 template bool igl::copyleft::cgal::mesh_to_polyhedron<CGAL::Polyhedron_3<CGAL::Simple_cartesian<double>, CGAL::Polyhedron_items_with_id_3, CGAL::HalfedgeDS_default, std::allocator<int> > >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::Polyhedron_3<CGAL::Simple_cartesian<double>, CGAL::Polyhedron_items_with_id_3, CGAL::HalfedgeDS_default, std::allocator<int> >&);

include/igl/copyleft/cgal/minkowski_sum.cpp

@@ -12,7 +12,6 @@
 #include "../../unique.h"
 #include "../../get_seconds.h"
 #include "../../edges.h"
-#include "assign_scalar.h"
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <cassert>
 #include <vector>

include/igl/copyleft/cgal/order_facets_around_edge.cpp

@@ -201,7 +201,7 @@ void igl::copyleft::cgal::order_facets_around_edge(
     return order;
   };
 
-  Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
+  DerivedI positive_order, negative_order;
   order_facets_around_edge(V, F, s, d, positive_side, positive_order, debug);
   order_facets_around_edge(V, F, s, d, negative_side, negative_order, debug);
   std::vector<size_t> tie_positive_order = index_sort(tie_positive_oriented);
@@ -399,7 +399,7 @@ void igl::copyleft::cgal::order_facets_around_edge(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 template void igl::copyleft::cgal::order_facets_around_edge<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, 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&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool);
 template void igl::copyleft::cgal::order_facets_around_edge<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&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, bool);
 template void igl::copyleft::cgal::order_facets_around_edge<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, 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&, unsigned long, unsigned long, std::vector<int, std::allocator<int> > 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> >&);

include/igl/copyleft/cgal/order_facets_around_edges.cpp

@@ -320,7 +320,7 @@ IGL_INLINE void igl::copyleft::cgal::order_facets_around_edges(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // 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

include/igl/copyleft/cgal/outer_element.cpp

@@ -207,7 +207,7 @@ IGL_INLINE void igl::copyleft::cgal::outer_edge(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 template void igl::copyleft::cgal::outer_edge<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -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<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<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<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>, 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<int, -1, 1, 0, -1, 1> > const&, long&, long&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);

include/igl/copyleft/cgal/outer_facet.cpp

@@ -148,7 +148,7 @@ IGL_INLINE void igl::copyleft::cgal::outer_facet(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // 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/igl/copyleft/cgal/outer_hull.cpp

@@ -1,14 +1,14 @@
 // 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"
 #include "remesh_self_intersections.h"
-#include "assign_scalar.h"
+#include "assign.h"
 #include "../../remove_unreferenced.h"
 
 #include <CGAL/AABB_tree.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;
@@ -94,14 +94,8 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull(
   // Remove unreferenced vertices and re-index faces
   {
     // Cast to output type
-    DerivedHV Vr_cast(Vr.rows(),Vr.cols());
-    for(int i = 0;i<Vr.rows();i++)
-    {
-      for(int j = 0;j<Vr.cols();j++)
-      {
-        assign_scalar(Vr(i,j), Vr_cast(i,j));
-      }
-    }
+    DerivedHV Vr_cast;
+    assign(Vr,Vr_cast);
     Eigen::VectorXi I;
     remove_unreferenced(Vr_cast,DerivedHF(HF),HV,HF,I);
   }
@@ -532,12 +526,10 @@ IGL_INLINE void igl::copyleft::cgal::outer_hull_legacy(
 
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 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> >&);
 #endif

include/igl/copyleft/cgal/peel_outer_hull_layers.cpp

@@ -114,7 +114,7 @@ IGL_INLINE size_t igl::copyleft::cgal::peel_outer_hull_layers(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 // generated by autoexplicit.sh
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

include/igl/copyleft/cgal/point_mesh_squared_distance.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "point_mesh_squared_distance.h"
 #include "mesh_to_cgal_triangle_list.h"
+#include "assign_scalar.h"
 
 template <
   typename Kernel,
@@ -131,5 +132,6 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
+template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Epeck, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<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<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<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/copyleft/cgal/point_segment_squared_distance.cpp

@@ -0,0 +1,41 @@
+// 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 "point_segment_squared_distance.h"
+
+template < typename Kernel>
+IGL_INLINE void igl::copyleft::cgal::point_segment_squared_distance(
+  const CGAL::Point_3<Kernel> & P1,
+  const CGAL::Segment_3<Kernel> & S2,
+  CGAL::Point_3<Kernel> & P2,
+  typename Kernel::FT & d)
+{
+  if(S2.is_degenerate())
+  {
+    P2 = S2.source();
+    d = (P1-P2).squared_length();
+    return;
+  }
+  // http://stackoverflow.com/a/1501725/148668
+  const auto sqr_len = S2.squared_length();
+  assert(sqr_len != 0);
+  const auto & V = S2.source();
+  const auto & W = S2.target();
+  const auto t = (P1-V).dot(W-V)/sqr_len;
+  if(t<0)
+  {
+    P2 = V;
+  }else if(t>1)
+  {
+    P2 = W;
+  }else
+  {
+    P2 = V  + t*(W-V);
+  }
+  d = (P1-P2).squared_length();
+}
+

include/igl/copyleft/cgal/point_segment_squared_distance.h

@@ -0,0 +1,44 @@
+// 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/.
+#ifndef IGL_COPYLEFT_CGAL_POINT_SEGMENT_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_POINT_SEGMENT_SQUARED_DISTANCE_H
+#include "../../igl_inline.h"
+#include <CGAL/Segment_3.h>
+#include <CGAL/Point_3.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given a point P1 and segment S2 find the points on each of closest
+      // approach and the squared distance thereof.
+      // 
+      // Inputs:
+      //   P1  point
+      //   S2  segment
+      // Outputs:
+      //   P2  point on S2 closest to P1
+      //   d  distance betwee P1 and S2
+      template < typename Kernel>
+      IGL_INLINE void point_segment_squared_distance(
+          const CGAL::Point_3<Kernel> & P1,
+          const CGAL::Segment_3<Kernel> & S2,
+          CGAL::Point_3<Kernel> & P2,
+          typename Kernel::FT & d
+          );
+
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_segment_squared_distance.cpp"
+#endif
+
+#endif
+

include/igl/copyleft/cgal/point_triangle_squared_distance.cpp

@@ -0,0 +1,47 @@
+// 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 "point_triangle_squared_distance.h"
+#include <CGAL/Segment_3.h>
+template < typename Kernel>
+IGL_INLINE void point_triangle_squared_distance(
+  const CGAL::Point_3<Kernel> & P1,
+  const CGAL::Triangle_3<Kernel> & T2,
+  CGAL::Point_3<Kernel> & P2,
+  typename Kernel::FT & d)
+{
+  assert(!T2.is_degenerate());
+  if(T2.has_on(P1))
+  {
+    P2 = P1;
+    d = 0;
+    return;
+  }
+  const auto proj_1 = T2.supporting_plane().projection(P2);
+  if(T2.has_on(proj_1))
+  {
+    P2 = proj_1;
+    d = (proj_1-P1).squared_length();
+    return;
+  }
+  // closest point must be on the boundary
+  bool first = true;
+  // loop over edges
+  for(int i=0;i<3;i++)
+  {
+    CGAL::Point_3<Kernel> P2i;
+    typename Kernel::FT di;
+    const CGAL::Segment_3<Kernel> si( T2.vertex(i+1), T2.vertex(i+2));
+    point_segment_squared_distance(P1,si,P2i,di);
+    if(first || di < d)
+    {
+      first = false;
+      d = di;
+      P2 = P2i;
+    }
+  }
+}

include/igl/copyleft/cgal/point_triangle_squared_distance.h

@@ -0,0 +1,45 @@
+// 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/.
+#ifndef IGL_COPYLEFT_CGAL_POINT_TRIANGLE_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_POINT_TRIANGLE_SQUARED_DISTANCE_H
+#include "../../igl_inline.h"
+#include <CGAL/Triangle_3.h>
+#include <CGAL/Point_3.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given a point P1 and triangle T2 find the points on each of closest
+      // approach and the squared distance thereof.
+      // 
+      // Inputs:
+      //   P1  point
+      //   T2  triangle
+      // Outputs:
+      //   P2  point on T2 closest to P1
+      //   d  distance betwee P1 and T2
+      template < typename Kernel>
+      IGL_INLINE void point_triangle_squared_distance(
+        const CGAL::Point_3<Kernel> & P1,
+        const CGAL::Triangle_3<Kernel> & T2,
+        CGAL::Point_3<Kernel> & P2,
+        typename Kernel::FT & d
+        );
+
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_triangle_squared_distance.cpp"
+#endif
+
+#endif
+
+

include/igl/copyleft/cgal/points_inside_component.cpp

@@ -338,7 +338,7 @@ IGL_INLINE void igl::copyleft::cgal::points_inside_component(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 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<   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<   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<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, 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> >&);

include/igl/copyleft/cgal/polyhedron_to_mesh.cpp

@@ -57,7 +57,7 @@ IGL_INLINE void igl::copyleft::cgal::polyhedron_to_mesh(
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
+// Explicit template instantiation
 #include <CGAL/Simple_cartesian.h>
 #include <CGAL/Polyhedron_items_with_id_3.h>
 template void igl::copyleft::cgal::polyhedron_to_mesh<CGAL::Polyhedron_3<CGAL::Simple_cartesian<double>, CGAL::Polyhedron_items_with_id_3, CGAL::HalfedgeDS_default, std::allocator<int> > >(CGAL::Polyhedron_3<CGAL::Simple_cartesian<double>, CGAL::Polyhedron_items_with_id_3, CGAL::HalfedgeDS_default, std::allocator<int> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);

include/igl/copyleft/cgal/projected_cdt.cpp

@@ -0,0 +1,82 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// 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 "projected_cdt.h"
+#include "insert_into_cdt.h"
+#include "assign_scalar.h"
+#include "../../list_to_matrix.h"
+template <typename Kernel, typename Index>
+IGL_INLINE void igl::copyleft::cgal::projected_cdt(
+  const std::vector<CGAL::Object> & objects,
+  const CGAL::Plane_3<Kernel> & P,
+  std::vector<CGAL::Point_3<Kernel> >& vertices,
+  std::vector<std::vector<Index> >& faces)
+{
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Segment_3<Kernel>  Segment_3; 
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef CGAL::Plane_3<Kernel>    Plane_3;
+  typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
+  typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
+  typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
+  typedef CGAL::Exact_intersections_tag Itag;
+  typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> CDT_2;
+  typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
+  CDT_plus_2 cdt;
+  for(const auto & obj : objects) insert_into_cdt(obj,P,cdt);
+  // Read off vertices of the cdt, remembering index
+  std::map<typename CDT_plus_2::Vertex_handle,Index> v2i;
+  size_t count=0;
+  for (
+    auto itr = cdt.finite_vertices_begin();
+    itr != cdt.finite_vertices_end(); 
+    itr++) 
+  {
+    vertices.push_back(P.to_3d(itr->point()));
+    v2i[itr] = count;
+    count++;
+  }
+  // Read off faces and store index triples
+  for (
+    auto itr = cdt.finite_faces_begin();
+    itr != cdt.finite_faces_end(); 
+    itr++)
+  {
+    faces.push_back( 
+      { v2i[itr->vertex(0)], v2i[itr->vertex(1)], v2i[itr->vertex(2)] });
+  }
+}
+
+template < typename Kernel, typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::copyleft::cgal::projected_cdt(
+  const std::vector<CGAL::Object> & objects,
+  const CGAL::Plane_3<Kernel> & P,
+  Eigen::PlainObjectBase<DerivedV> & V,
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  std::vector<CGAL::Point_3<Kernel> > vertices;
+  std::vector<std::vector<typename DerivedF::Scalar> > faces;
+  projected_cdt(objects,P,vertices,faces);
+  V.resize(vertices.size(),3);
+  for(int v = 0;v<vertices.size();v++)
+  {
+    for(int d = 0;d<3;d++)
+    {
+      assign_scalar(vertices[v][d], V(v,d));
+    }
+  }
+  list_to_matrix(faces,F);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::projected_cdt<CGAL::Epick, long>(std::vector<CGAL::Object, std::allocator<CGAL::Object> > const&, CGAL::Plane_3<CGAL::Epick> const&, std::vector<CGAL::Point_3<CGAL::Epick>, std::allocator<CGAL::Point_3<CGAL::Epick> > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::projected_cdt<CGAL::Epeck, long>(std::vector<CGAL::Object, std::allocator<CGAL::Object> > const&, CGAL::Plane_3<CGAL::Epeck> const&, std::vector<CGAL::Point_3<CGAL::Epeck>, std::allocator<CGAL::Point_3<CGAL::Epeck> > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#endif

include/igl/copyleft/cgal/projected_cdt.h

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_COPYLEFT_CGAL_PROJECTED_CDT_H
+#define IGL_COPYLEFT_CGAL_PROJECTED_CDT_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <CGAL/Plane_3.h>
+#include <CGAL/Point_3.h>
+#include <vector>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given a list of objects (e.g., resulting from intersecting a triangle
+      // with many other triangles), construct a constrained Delaunay
+      // triangulation on a given plane (P), by inersting constraints for each
+      // object projected onto that plane.
+      //
+      // Inputs:
+      //   objects  list of objects. This should lie on the given plane (P),
+      //     otherwise they are added to the cdt _after_ their non-trivial
+      //     projection
+      //   P  plane upon which all objects lie and upon which the CDT is
+      //     conducted
+      // Outputs:
+      //   vertices  list of vertices of the CDT mesh _back on the 3D plane_
+      //   faces  list of list of triangle indices into vertices
+      //   
+      template <typename Kernel, typename Index>
+      IGL_INLINE void projected_cdt(
+        const std::vector<CGAL::Object> & objects,
+        const CGAL::Plane_3<Kernel> & P,
+        std::vector<CGAL::Point_3<Kernel> >& vertices,
+        std::vector<std::vector<Index> >& faces);
+      // Outputs:
+      //   V  #V by 3 list of vertices of the CDT mesh _back on the 3D plane_,
+      //     **cast** from the number type of Kernel to the number type of
+      //     DerivedV
+      //   F  #F by 3 list of triangle indices into V
+      template < typename Kernel, typename DerivedV, typename DerivedF>
+      IGL_INLINE void projected_cdt(
+        const std::vector<CGAL::Object> & objects,
+        const CGAL::Plane_3<Kernel> & P,
+        Eigen::PlainObjectBase<DerivedV> & V,
+        Eigen::PlainObjectBase<DerivedF> & F);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "projected_cdt.cpp"
+#endif
+#endif