merge

Former-commit-id: 1c72eafdd18eb45f4aeb03a28534c6b44892673f

README.md

@@ -190,6 +190,8 @@ Libigl is used by many research groups around the world. In 2015, it won the
 Eurographics/ACM Symposium on Geometry Processing software award. Here are a
 few labs/companies/institutions using libigl:
 
+ - [Adobe Research](http://www.adobe.com/technology/)  
+ - [Pixar Research](http://graphics.pixar.com/research/)
  - [Spine by Esoteric Software](http://esotericsoftware.com/) is an animation tool dedicated to 2D characters.
  - Columbia University, [Columbia Computer Graphics Group](http://www.cs.columbia.edu/cg/), USA
  - [Cornell University](http://www.graphics.cornell.edu/), USA
@@ -200,13 +202,18 @@ few labs/companies/institutions using libigl:
  - [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
- - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
+ - [TU Berlin](https://www.cg.tu-berlin.de), Germany
  - [TU Delft](http://www.tudelft.nl/en/), Netherlands
+ - [Telecom ParisTech](http://www.telecom-paristech.fr/en/formation-et-innovation-dans-le-numerique.html), Paris, France
  - [Universidade Federal de Santa Catarina](http://mtm.ufsc.br/~leo/), Brazil
- - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
  - [University College London](http://vecg.cs.ucl.ac.uk/), England
+ - [University of California Berkeley](http://vis.berkeley.edu/), USA
  - [University of Cambridge](http://www.cam.ac.uk/), England
  - [University of Pennsylvania](http://cg.cis.upenn.edu/), USA
+ - [University of Texas at Austin](http://www.cs.utexas.edu/users/evouga/), USA
+ - [University of Victoria](https://www.csc.uvic.ca/Research/graphics/), Canada
+ - [Università della Svizzera Italiana](http://www.usi.ch/en), Switzerland
+ - [Université Toulouse III Paul Sabatier](http://www.univ-tlse3.fr/), France
 
 
 ## Contact

include/igl/AABB.h

@@ -8,6 +8,7 @@
 #ifndef IGL_AABB_H
 #define IGL_AABB_H
 
+#include "Hit.h"
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 #include <vector>
@@ -134,7 +135,8 @@ public:
           const Eigen::VectorXi & I);
       // Return whether at leaf node
       inline bool is_leaf() const;
-      // Find the indices of elements containing given point.
+      // Find the indices of elements containing given point: this makes sense
+      // when Ele is a co-dimension 0 simplex (tets in 3D, triangles in 2D).
       //
       // Inputs:
       //   V  #V by dim list of mesh vertex positions. **Should be same as used to
@@ -185,19 +187,43 @@ public:
       // Known bugs: currently assumes Elements are triangles regardless of
       // dimension.
       inline Scalar squared_distance(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::MatrixXi & Ele, 
-          const RowVectorDIMS & p,
-          int & i,
-          RowVectorDIMS & c) const;
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & p,
+        int & i,
+        RowVectorDIMS & c) const;
 //private:
       inline Scalar squared_distance(
-          const Eigen::PlainObjectBase<DerivedV> & V,
-          const Eigen::MatrixXi & Ele, 
-          const RowVectorDIMS & p,
-          const Scalar min_sqr_d,
-          int & i,
-          RowVectorDIMS & c) const;
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & p,
+        const Scalar min_sqr_d,
+        int & i,
+        RowVectorDIMS & c) const;
+      // All hits
+      inline bool intersect_ray(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        std::vector<igl::Hit> & hits) const;
+      // First hit
+      inline bool intersect_ray(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        igl::Hit & hit) const;
+//private:
+      inline bool intersect_ray(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::MatrixXi & Ele, 
+        const RowVectorDIMS & origin,
+        const RowVectorDIMS & dir,
+        const Scalar min_t,
+        igl::Hit & hit) const;
+
+
 public:
       template <
         typename DerivedP, 
@@ -281,10 +307,14 @@ public:
 #include "project_to_line_segment.h"
 #include "sort.h"
 #include "volume.h"
+#include "ray_box_intersect.h"
+#include "ray_mesh_intersect.h"
 #include <iostream>
 #include <iomanip>
 #include <limits>
 #include <list>
+#include <queue>
+#include <stack>
 
 template <typename DerivedV, int DIM>
   template <typename Derivedbb_mins, typename Derivedbb_maxs>
@@ -1149,4 +1179,160 @@ igl::AABB<DerivedV,DIM>::barycentric_coordinates(
   bary(0) = 1.0f - bary(1) - bary(2);
 }
 
+template <typename DerivedV, int DIM>
+inline bool 
+igl::AABB<DerivedV,DIM>::intersect_ray(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  std::vector<igl::Hit> & hits) const
+{
+  hits.clear();
+  const Scalar t0 = 0;
+  const Scalar t1 = std::numeric_limits<Scalar>::infinity();
+  {
+    Scalar _1,_2;
+    if(!ray_box_intersect(origin,dir,m_box,t0,t1,_1,_2))
+    {
+      return false;
+    }
+  }
+  if(this->is_leaf())
+  {
+    // Actually process elements
+    assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
+    // Cheesecake way of hitting element
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hits);
+  }
+  std::vector<igl::Hit> left_hits;
+  std::vector<igl::Hit> right_hits;
+  const bool left_ret = m_left->intersect_ray(V,Ele,origin,dir,left_hits);
+  const bool right_ret = m_right->intersect_ray(V,Ele,origin,dir,right_hits);
+  hits.insert(hits.end(),left_hits.begin(),left_hits.end());
+  hits.insert(hits.end(),right_hits.begin(),right_hits.end());
+  return left_ret || right_ret;
+}
+
+template <typename DerivedV, int DIM>
+inline bool 
+igl::AABB<DerivedV,DIM>::intersect_ray(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  igl::Hit & hit) const
+{
+#if false
+  // BFS
+  std::queue<const AABB *> Q;
+  // Or DFS
+  //std::stack<const AABB *> Q;
+  Q.push(this);
+  bool any_hit = false;
+  hit.t = std::numeric_limits<Scalar>::infinity();
+  while(!Q.empty())
+  {
+    const AABB * tree = Q.front();
+    //const AABB * tree = Q.top();
+    Q.pop();
+    {
+      Scalar _1,_2;
+      if(!ray_box_intersect(
+        origin,dir,tree->m_box,Scalar(0),Scalar(hit.t),_1,_2))
+      {
+        continue;
+      }
+    }
+    if(tree->is_leaf())
+    {
+      // Actually process elements
+      assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
+      igl::Hit leaf_hit;
+      if(
+        ray_mesh_intersect(origin,dir,V,Ele.row(tree->m_primitive).eval(),leaf_hit)&&
+        leaf_hit.t < hit.t)
+      {
+        hit = leaf_hit;
+      }
+      continue;
+    }
+    // Add children to queue
+    Q.push(tree->m_left);
+    Q.push(tree->m_right);
+  }
+  return any_hit;
+#else
+  // DFS
+  return intersect_ray(
+    V,Ele,origin,dir,std::numeric_limits<Scalar>::infinity(),hit);
+#endif
+}
+
+template <typename DerivedV, int DIM>
+inline bool 
+igl::AABB<DerivedV,DIM>::intersect_ray(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::MatrixXi & Ele, 
+  const RowVectorDIMS & origin,
+  const RowVectorDIMS & dir,
+  const Scalar _min_t,
+  igl::Hit & hit) const
+{
+  //// Naive, slow
+  //std::vector<igl::Hit> hits;
+  //intersect_ray(V,Ele,origin,dir,hits);
+  //if(hits.size() > 0)
+  //{
+  //  hit = hits.front();
+  //  return true;
+  //}else
+  //{
+  //  return false;
+  //}
+  Scalar min_t = _min_t;
+  const Scalar t0 = 0;
+  {
+    Scalar _1,_2;
+    if(!ray_box_intersect(origin,dir,m_box,t0,min_t,_1,_2))
+    {
+      return false;
+    }
+  }
+  if(this->is_leaf())
+  {
+    // Actually process elements
+    assert((Ele.size() == 0 || Ele.cols() == 3) && "Elements should be triangles");
+    // Cheesecake way of hitting element
+    return ray_mesh_intersect(origin,dir,V,Ele.row(m_primitive).eval(),hit);
+  }
+
+  // Doesn't seem like smartly choosing left before/after right makes a
+  // differnce
+  igl::Hit left_hit;
+  igl::Hit right_hit;
+  bool left_ret = m_left->intersect_ray(V,Ele,origin,dir,min_t,left_hit);
+  if(left_ret && left_hit.t<min_t)
+  {
+    // It's scary that this line doesn't seem to matter....
+    min_t = left_hit.t;
+    hit = left_hit;
+    left_ret = true;
+  }else
+  {
+    left_ret = false;
+  }
+  bool right_ret = m_right->intersect_ray(V,Ele,origin,dir,min_t,right_hit);
+  if(right_ret && right_hit.t<min_t)
+  {
+    min_t = right_hit.t;
+    hit = right_hit;
+    right_ret = true;
+  }else
+  {
+    right_ret = false;
+  }
+  return left_ret || right_ret;
+}
+
 #endif

include/igl/Camera.h

@@ -66,11 +66,13 @@ namespace igl
       //     glMultMatrixd(projection().data());
       //
       inline Eigen::Matrix4d projection() const;
-      // Return an Affine transformation (rigid actually) that takes a world 3d coordinate and
-      // transforms it into the relative camera coordinates.
+      // Return an Affine transformation (rigid actually) that 
+      // takes relative coordinates and tramsforms them into world 3d
+      // coordinates: moves the camera into the scene.
       inline Eigen::Affine3d affine() const;
-      // Return an Affine transformation (rigid actually) that takes relative
-      // coordinates and tramsforms them into world 3d coordinates.
+      // Return an Affine transformation (rigid actually) that puts the takes a
+      // world 3d coordinate and transforms it into the relative camera
+      // coordinates: moves the scene in front of the camera.
       //
       // Note:
       //
@@ -81,7 +83,7 @@ namespace igl
       //
       // Is equivalent to
       //
-      //     glMultMatrixd(camera.affine().matrix().data());
+      //     glMultMatrixd(camera.inverse().matrix().data());
       //
       // See also: affine, eye, at, up
       inline Eigen::Affine3d inverse() const;

include/igl/Hit.h

@@ -0,0 +1,25 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
+//               2014 Christian Schüller <schuellchr@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_HIT_H
+#define IGL_HIT_H
+
+namespace igl
+{
+  // Reimplementation of the embree::Hit struct from embree1.0
+  // 
+  // TODO: template on floating point type
+  struct Hit
+  {
+    int id; // primitive id
+    int gid; // geometry id
+    float u,v; // barycentric coordinates
+    float t; // distance = direction*t to intersection
+  };
+}
+#endif 

include/igl/Singular_Value_Decomposition_Preamble.hpp

@@ -19,6 +19,17 @@
 #include <iostream>
 #endif
 
+// Prevent warnings
+#ifdef ENABLE_SCALAR_IMPLEMENTATION
+#  undef ENABLE_SCALAR_IMPLEMENTATION
+#endif
+#ifdef ENABLE_SSE_IMPLEMENTATION
+#  undef ENABLE_SSE_IMPLEMENTATION
+#endif
+#ifdef ENABLE_AVX_IMPLEMENTATION
+#  undef ENABLE_AVX_IMPLEMENTATION
+#endif
+
 #ifdef USE_SCALAR_IMPLEMENTATION
 #define ENABLE_SCALAR_IMPLEMENTATION(X) X
 #else

include/igl/active_set.cpp

@@ -67,8 +67,8 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   assert((Aieq.size() == 0 && Bieq.size() == 0) || Aieq.cols() == n);
   assert((Aieq.size() == 0 && Bieq.size() == 0) || Aieq.rows() == Bieq.rows());
   assert((Aieq.size() == 0 && Bieq.size() == 0) || Bieq.cols() == 1);
-  Eigen::PlainObjectBase<Derivedlx> lx;
-  Eigen::PlainObjectBase<Derivedux> ux;
+  Eigen::Matrix<typename Derivedlx::Scalar,Eigen::Dynamic,1> lx;
+  Eigen::Matrix<typename Derivedux::Scalar,Eigen::Dynamic,1> ux;
   if(p_lx.size() == 0)
   {
     lx = Eigen::PlainObjectBase<Derivedlx>::Constant(

include/igl/ambient_occlusion.cpp

@@ -0,0 +1,135 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "ambient_occlusion.h"
+#include "random_dir.h"
+#include "ray_mesh_intersect.h"
+#include "EPS.h"
+#include "Hit.h"
+
+template <
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::ambient_occlusion(
+  const std::function<
+    bool(
+      const Eigen::Vector3f&,
+      const Eigen::Vector3f&)
+      > & shoot_ray,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  using namespace Eigen;
+  using namespace igl;
+  const int n = P.rows();
+  // Resize output
+  S.resize(n,1);
+  VectorXi hits = VectorXi::Zero(n,1);
+  // Embree seems to be parallel when constructing but not when tracing rays
+#pragma omp parallel for
+  const MatrixXf D = random_dir_stratified(num_samples).cast<float>();
+  // loop over mesh vertices
+  for(int p = 0;p<n;p++)
+  {
+    const Vector3f origin = P.row(p).template cast<float>();
+    const Vector3f normal = N.row(p).template cast<float>();
+    int num_hits = 0;
+    for(int s = 0;s<num_samples;s++)
+    {
+//      //Vector3d d = random_dir();
+      Vector3f d = D.row(s);
+      if(d.dot(normal) < 0)
+      {
+        // reverse ray
+        d *= -1;
+      }
+      if(shoot_ray(origin,d))
+      {
+        num_hits++;
+      }
+    }
+    S(p) = (double)num_hits/(double)num_samples;
+  }
+}
+
+template <
+  typename DerivedV,
+  int DIM,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::ambient_occlusion(
+  const igl::AABB<DerivedV,DIM> & aabb,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  const auto & shoot_ray = [&aabb,&V,&F](
+    const Eigen::Vector3f& _s,
+    const Eigen::Vector3f& dir)->bool
+  {
+    Eigen::Vector3f s = _s+1e-4*dir;
+    igl::Hit hit;
+    return aabb.intersect_ray(
+      V,
+      F,
+      s  .cast<typename DerivedV::Scalar>().eval(),
+      dir.cast<typename DerivedV::Scalar>().eval(),
+      hit);
+  };
+  return ambient_occlusion(shoot_ray,P,N,num_samples,S);
+
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::ambient_occlusion(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  if(F.rows() < 100)
+  {
+    // Super naive
+    const auto & shoot_ray = [&V,&F](
+      const Eigen::Vector3f& _s,
+      const Eigen::Vector3f& dir)->bool
+    {
+      Eigen::Vector3f s = _s+1e-4*dir;
+      igl::Hit hit;
+      return ray_mesh_intersect(s,dir,V,F,hit);
+    };
+    return ambient_occlusion(shoot_ray,P,N,num_samples,S);
+  }
+  AABB<DerivedV,3> aabb;
+  aabb.init(V,F);
+  return ambient_occlusion(aabb,V,F,P,N,num_samples,S);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// 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
+template void igl::ambient_occlusion<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, 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, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template void igl::ambient_occlusion<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, 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, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/ambient_occlusion.h

@@ -0,0 +1,80 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_AMBIENT_OCCLUSION_H
+#define IGL_AMBIENT_OCCLUSION_H
+#include "igl_inline.h"
+#include "AABB.h"
+#include <Eigen/Core>
+#include <functional>
+namespace igl
+{
+  // Compute ambient occlusion per given point
+  //
+  // Inputs:
+  //    shoot_ray  function handle that outputs hits of a given ray against a
+  //      mesh (embedded in function handles as captured variable/data)
+  //    P  #P by 3 list of origin points
+  //    N  #P by 3 list of origin normals
+  // Outputs:
+  //    S  #P list of ambient occlusion values between 1 (fully occluded) and
+  //      0 (not occluded)
+  //
+  template <
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void ambient_occlusion(
+    const std::function<
+      bool(
+        const Eigen::Vector3f&,
+        const Eigen::Vector3f&)
+        > & shoot_ray,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  // Inputs:
+  //   AABB  axis-aligned bounding box hierarchy around (V,F)
+  template <
+    typename DerivedV,
+    int DIM,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void ambient_occlusion(
+    const igl::AABB<DerivedV,DIM> & aabb,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  // Inputs:
+  //    V  #V by 3 list of mesh vertex positions
+  //    F  #F by 3 list of mesh face indices into V
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void ambient_occlusion(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+
+};
+#ifndef IGL_STATIC_LIBRARY
+#  include "ambient_occlusion.cpp"
+#endif
+
+#endif

include/igl/angle_bound_frame_fields.cpp

@@ -38,6 +38,7 @@ namespace igl {
       Eigen::VectorXi indInteriorToFull;
       Eigen::VectorXi indFullToInterior;
 
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
       Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
 
       //laplacians

include/igl/arap_dof.cpp

@@ -14,7 +14,6 @@
 #include "repmat.h"
 #include "slice.h"
 #include "colon.h"
-#include "full.h"
 #include "is_sparse.h"
 #include "mode.h"
 #include "is_symmetric.h"
@@ -203,9 +202,8 @@ IGL_INLINE bool igl::arap_dof_precomputation(
       if(is_sparse(CSMjM_i))
       {
         // Convert to full
-        MatrixXd CSMjM_ifull;
         //printf("CSM_M(): full\n");
-        full(CSMjM_i,CSMjM_ifull);
+        MatrixXd CSMjM_ifull(CSMjM_i);
 //        printf("CSM_M[%d]: %d %d\n",i,data.CSM_M[i].rows(),data.CSM_M[i].cols());
 //        printf("CSM_M[%d].block(%d*%d=%d,0,%d,%d): %d %d\n",i,j,k,CSMjM_i.rows(),CSMjM_i.cols(),
 //            data.CSM_M[i].block(j*k,0,CSMjM_i.rows(),CSMjM_i.cols()).rows(),
@@ -546,8 +544,8 @@ IGL_INLINE bool igl::arap_dof_recomputation(
 
   // Compute dense solve matrix (alternative of matrix factorization)
   //printf("min_quad_dense_precompute()\n");
-  MatrixXd Qfull; full(*Q, Qfull);  
-  MatrixXd A_eqfull; full(A_eq, A_eqfull);
+  MatrixXd Qfull(*Q);
+  MatrixXd A_eqfull(A_eq);
   MatrixXd M_Solve;
 
   double timer0_start = get_seconds_hires();

include/igl/average_onto_vertices.cpp

@@ -7,15 +7,15 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "average_onto_vertices.h"
 
-template <typename T, typename I>
-IGL_INLINE void igl::average_onto_vertices(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
-            const Eigen::Matrix<I, Eigen::Dynamic, Eigen::Dynamic> &F,
-            const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &S,
-            Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &SV)
+template<typename DerivedV,typename DerivedF,typename DerivedS>
+IGL_INLINE void igl::average_onto_vertices(const Eigen::MatrixBase<DerivedV> &V,
+  const Eigen::MatrixBase<DerivedF> &F,
+  const Eigen::MatrixBase<DerivedS> &S,
+  Eigen::MatrixBase<DerivedS> &SV)
 {
-
-  SV = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>::Zero(V.rows(),S.cols());
-  Eigen::Matrix<T, Eigen::Dynamic, 1> COUNT = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(V.rows());
+  SV = DerivedS::Zero(V.rows(),S.cols());
+  Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> COUNT(V.rows());
+  COUNT.setZero();
   for (int i = 0; i <F.rows(); ++i)
   {
     for (int j = 0; j<F.cols(); ++j)
@@ -26,10 +26,8 @@ IGL_INLINE void igl::average_onto_vertices(const Eigen::Matrix<T, Eigen::Dynamic
   }
   for (int i = 0; i <V.rows(); ++i)
     SV.row(i) /= COUNT[i];
-
 };
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::average_onto_vertices<double, int>(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<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/average_onto_vertices.h

@@ -21,12 +21,11 @@ namespace igl
   // 
   // Output:
   // SV: scalar field defined on vertices
-  template <typename T, typename I>
-  IGL_INLINE void average_onto_vertices(
-    const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
-    const Eigen::Matrix<I, Eigen::Dynamic, Eigen::Dynamic> &F,
-    const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &S,
-    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &SV);
+  template<typename DerivedV,typename DerivedF,typename DerivedS>
+  IGL_INLINE void average_onto_vertices(const Eigen::MatrixBase<DerivedV> &V,
+    const Eigen::MatrixBase<DerivedF> &F,
+    const Eigen::MatrixBase<DerivedS> &S,
+    Eigen::MatrixBase<DerivedS> &SV);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/boundary_loop.cpp

@@ -23,10 +23,10 @@ IGL_INLINE void igl::boundary_loop(
   if(F.rows() == 0)
     return;
 
-  MatrixXd Vdummy(F.maxCoeff()+1,1);
-  MatrixXi TT,TTi;
+  VectorXd Vdummy(F.maxCoeff()+1,1);
+  DerivedF TT,TTi;
   vector<std::vector<int> > VF, VFi;
-  triangle_triangle_adjacency(Vdummy,F,TT,TTi);
+  triangle_triangle_adjacency(F,TT,TTi);
   vertex_triangle_adjacency(Vdummy,F,VF,VFi);
 
   vector<bool> unvisited = is_border_vertex(Vdummy,F);

include/igl/comb_cross_field.cpp

@@ -27,6 +27,7 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
   private:
@@ -72,7 +73,7 @@ namespace igl {
     PD2(_PD2)
     {
       igl::per_face_normals(V,F,N);
-      igl::triangle_triangle_adjacency(V,F,TT,TTi);
+      igl::triangle_triangle_adjacency(F,TT,TTi);
     }
     inline void comb(Eigen::PlainObjectBase<DerivedV> &PD1out,
               Eigen::PlainObjectBase<DerivedV> &PD2out)

include/igl/comb_frame_field.cpp

@@ -24,7 +24,7 @@ IGL_INLINE void igl::comb_frame_field(const Eigen::PlainObjectBase<DerivedV> &V,
                                       Eigen::PlainObjectBase<DerivedP> &PD1_combed,
                                       Eigen::PlainObjectBase<DerivedP> &PD2_combed)
 {
-  Eigen::PlainObjectBase<DerivedV> B1, B2, B3;
+  DerivedV B1, B2, B3;
   igl::local_basis(V,F,B1,B2,B3);
 
   PD1_combed.resize(BIS1_combed.rows(),3);

include/igl/comb_line_field.cpp

@@ -25,10 +25,12 @@ public:
     const Eigen::PlainObjectBase<DerivedV> &V;
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
 private:
     // internal
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 
@@ -65,7 +67,7 @@ public:
         PD1(_PD1)
     {
         igl::per_face_normals(V,F,N);
-        igl::triangle_triangle_adjacency(V,F,TT,TTi);
+        igl::triangle_triangle_adjacency(F,TT,TTi);
     }
 
     inline void comb(Eigen::PlainObjectBase<DerivedV> &PD1out)

include/igl/comiso/frame_field.cpp

@@ -115,7 +115,7 @@ FrameInterpolator::FrameInterpolator(const Eigen::MatrixXd& _V, const Eigen::Mat
 
 
   // Generate topological relations
-  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  igl::triangle_triangle_adjacency(F,TT,TTi);
   igl::edge_topology(V,F, EV, FE, EF);
 
   // Flag border edges

include/igl/comiso/miq.cpp

@@ -1191,7 +1191,7 @@ F(F_)
   igl::cut_mesh(V, F, Handle_Seams, Vcut, Fcut);
 
   igl::local_basis(V,F,B1,B2,B3);
-  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  igl::triangle_triangle_adjacency(F,TT,TTi);
 
   // Prepare indexing for the linear system
   VertexIndexing<DerivedV, DerivedF> VInd(V, F, Vcut, Fcut, TT, TTi, Handle_MMatch, Handle_Singular, Handle_Seams);

include/igl/comiso/nrosy.cpp

@@ -165,7 +165,7 @@ igl::comiso::NRosyField::NRosyField(const Eigen::MatrixXd& _V, const Eigen::Matr
 
 
   // Generate topological relations
-  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  igl::triangle_triangle_adjacency(F,TT,TTi);
   igl::edge_topology(V,F, EV, FE, EF);
 
   // Flag border edges

include/igl/compute_frame_field_bisectors.cpp

@@ -70,7 +70,7 @@ IGL_INLINE void igl::compute_frame_field_bisectors(
                                                    Eigen::PlainObjectBase<DerivedV>& BIS1,
                                                    Eigen::PlainObjectBase<DerivedV>& BIS2)
 {
-  Eigen::PlainObjectBase<DerivedV> B1, B2, B3;
+  DerivedV B1, B2, B3;
   igl::local_basis(V,F,B1,B2,B3);
 
   compute_frame_field_bisectors( V, F, B1, B2, PD1, PD2, BIS1, BIS2);

include/igl/copyleft/boolean/minkowski_sum.cpp

@@ -0,0 +1,229 @@
+#include "minkowski_sum.h"
+#include "mesh_boolean.h"
+
+#include "../../slice_mask.h"
+#include "../../unique.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename Deriveds,
+  typename Derivedd,
+  typename DerivedW,
+  typename DerivedG,
+  typename DerivedJ>
+IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
+  const Eigen::PlainObjectBase<DerivedVA> & VA,
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<Deriveds> & s,
+  const Eigen::PlainObjectBase<Derivedd> & d,
+  const bool resolve_overlaps, 
+  Eigen::PlainObjectBase<DerivedW> & W,
+  Eigen::PlainObjectBase<DerivedG> & G,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  using namespace Eigen;
+  using namespace std;
+  // silly base case
+  if(FA.size() == 0)
+  {
+    W.resize(0,3);
+    G.resize(0,3);
+    return;
+  }
+  const int dim = VA.cols();
+  assert(dim == 3 && "dim must be 3D");
+  assert(s.size() == 3 && "s must be 3D point");
+  assert(d.size() == 3 && "d must be 3D point");
+  // segment vector
+  const CGAL::Vector_3<CGAL::Epeck> v(d(0)-s(0),d(1)-s(1),d(2)-s(2));
+  // number of vertices
+  const int n = VA.rows();
+  // duplicate vertices at s and d, we'll remove unreferernced later
+  DerivedW WW(2*n,dim);
+  for(int i = 0;i<n;i++)
+  {
+    for(int j = 0;j<dim;j++)
+    {
+      WW  (i,j) = VA(i,j) + s(j);
+      WW(i+n,j) = VA(i,j) + d(j);
+    }
+  }
+  // number of faces
+  const int m = FA.rows();
+  // Mask whether positive dot product, or negative: because of exactly zero,
+  // these are not necessarily complementary
+  Matrix<bool,Dynamic,1> P(m,1),N(m,1);
+  // loop over faces
+  int mp = 0,mn = 0;
+  for(int f = 0;f<m;f++)
+  {
+    const CGAL::Plane_3<CGAL::Epeck> plane(
+      CGAL::Point_3<CGAL::Epeck>(VA(FA(f,0),0),VA(FA(f,0),1),VA(FA(f,0),2)),
+      CGAL::Point_3<CGAL::Epeck>(VA(FA(f,1),0),VA(FA(f,1),1),VA(FA(f,1),2)),
+      CGAL::Point_3<CGAL::Epeck>(VA(FA(f,2),0),VA(FA(f,2),1),VA(FA(f,2),2)));
+    const auto normal = plane.orthogonal_vector();
+    const auto dt = normal * v;
+    if(dt > 0)
+    {
+      P(f) = true;
+      N(f) = false;
+      mp++;
+    }else if(dt < 0)
+    {
+      P(f) = false;
+      N(f) = true;
+      mn++;
+    }else
+    {
+      P(f) = false;
+      N(f) = false;
+    }
+  }
+
+  typedef Matrix<typename DerivedG::Scalar,Dynamic,Dynamic> MatrixXI;
+  typedef Matrix<typename DerivedG::Scalar,Dynamic,1> VectorXI;
+  MatrixXI GT(mp+mn,3);
+  GT<< slice_mask(FA,N,1), slice_mask((FA.array()+n).eval(),P,1);
+  // J indexes FA for parts at s and m+FA for parts at d
+  J = DerivedJ::LinSpaced(m,0,m-1);
+  DerivedJ JT(mp+mn);
+  JT << slice_mask(J,P,1), slice_mask(J,N,1);
+  JT.block(mp,0,mn,1).array()+=m;
+
+  // Original non-co-planar faces with positively oriented reversed
+  MatrixXI BA(mp+mn,3);
+  BA << slice_mask(FA,P,1).rowwise().reverse(), slice_mask(FA,N,1);
+  // Quads along **all** sides
+  MatrixXI GQ((mp+mn)*3,4);
+  GQ<< 
+    BA.col(1), BA.col(0), BA.col(0).array()+n, BA.col(1).array()+n,
+    BA.col(2), BA.col(1), BA.col(1).array()+n, BA.col(2).array()+n,
+    BA.col(0), BA.col(2), BA.col(2).array()+n, BA.col(0).array()+n;
+
+  MatrixXI uGQ;
+  VectorXI S,sI,sJ;
+  //const auto & total_signed_distance = 
+  [](
+      const MatrixXI & F,
+      VectorXI & S,
+      MatrixXI & uF,
+      VectorXI & I,
+      VectorXI & J)
+  {
+    const int m = F.rows();
+    const int d = F.cols();
+    MatrixXI sF = F;
+    const auto MN = sF.rowwise().minCoeff().eval();
+    // rotate until smallest index is first
+    for(int p = 0;p<d;p++)
+    {
+      for(int f = 0;f<m;f++)
+      {
+        if(sF(f,0) != MN(f))
+        {
+          for(int r = 0;r<d-1;r++)
+          {
+            std::swap(sF(f,r),sF(f,r+1));
+          }
+        }
+      }
+    }
+    // swap orienation
+    for(int f = 0;f<m;f++)
+    {
+      if(sF(f,d-1) < sF(f,1))
+      {
+        sF.block(f,1,1,d-1) = sF.block(f,1,1,d-1).reverse().eval();
+      }
+    }
+    Matrix<bool,Dynamic,1> M = Matrix<bool,Dynamic,1>::Zero(m,1);
+    {
+      VectorXI P = VectorXI::LinSpaced(d,0,d-1);
+      for(int p = 0;p<d;p++)
+      {
+        for(int f = 0;f<m;f++)
+        {
+          bool all = true;
+          for(int r = 0;r<d;r++)
+          {
+            all = all && (sF(f,P(r)) == F(f,r));
+          }
+          M(f) = M(f) || all;
+        }
+        for(int r = 0;r<d-1;r++)
+        {
+          std::swap(P(r),P(r+1));
+        }
+      }
+    }
+    unique_rows(sF,uF,I,J);
+    S = VectorXI::Zero(uF.rows(),1);
+    assert(m == J.rows());
+    for(int f = 0;f<m;f++)
+    {
+      S(J(f)) += M(f) ? 1 : -1;
+    }
+  }(MatrixXI(GQ),S,uGQ,sI,sJ);
+  assert(S.rows() == uGQ.rows());
+  const int nq = (S.array().abs()==2).count();
+  GQ.resize(nq,4);
+  {
+    int k = 0;
+    for(int q = 0;q<uGQ.rows();q++)
+    {
+      switch(S(q))
+      {
+        case -2:
+          GQ.row(k++) = uGQ.row(q).reverse().eval();
+          break;
+        case 2:
+          GQ.row(k++) = uGQ.row(q);
+          break;
+        default:
+        // do not add
+          break;
+      }
+    }
+    assert(nq == k);
+  }
+
+  MatrixXI GG(GT.rows()+2*GQ.rows(),3);
+  GG<< 
+    GT,
+    GQ.col(0), GQ.col(1), GQ.col(2), 
+    GQ.col(0), GQ.col(2), GQ.col(3);
+  J.resize(JT.rows()+2*GQ.rows(),1);
+  J<<JT,DerivedJ::Constant(2*GQ.rows(),1,2*m+1);
+  if(resolve_overlaps)
+  {
+    DerivedJ SJ;
+    mesh_boolean(
+      WW,GG,
+      Matrix<typename DerivedVA::Scalar,Dynamic,Dynamic>(),MatrixXI(),
+      MESH_BOOLEAN_TYPE_UNION,
+      W,G,SJ);
+    J = slice(DerivedJ(J),SJ,1);
+  }
+}
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename Deriveds,
+  typename Derivedd,
+  typename DerivedW,
+  typename DerivedG,
+  typename DerivedJ>
+IGL_INLINE void igl::copyleft::boolean::minkowski_sum(
+  const Eigen::PlainObjectBase<DerivedVA> & VA,
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<Deriveds> & s,
+  const Eigen::PlainObjectBase<Derivedd> & d,
+  Eigen::PlainObjectBase<DerivedW> & W,
+  Eigen::PlainObjectBase<DerivedG> & G,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  return minkowski_sum(VA,FA,s,d,true,W,G,J);
+}

include/igl/copyleft/boolean/minkowski_sum.h

@@ -0,0 +1,71 @@
+#ifndef IGL_COPYLEFT_CGAL_MINKOWSKI_SUM_H
+#define IGL_COPYLEFT_CGAL_MINKOWSKI_SUM_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace boolean
+    {
+      // Compute the Minkowski sum of a closed triangle mesh (V,F) and a
+      // segment [s,d] in 3D.
+      //
+      // Inputs:
+      //   VA  #VA by 3 list of mesh vertices in 3D
+      //   FA  #FA by 3 list of triangle indices into VA
+      //   s  segment source endpoint in 3D
+      //   d  segment source endpoint in 3D
+      //   resolve_overlaps  whether or not to resolve self-union. If false
+      //     then result may contain self-intersections if input mesh is
+      //     non-convex.
+      // Outputs:
+      //   W  #W by 3 list of mesh vertices in 3D
+      //   G  #G by 3 list of triangle indices into W
+      //   J  #G list of indices into [F;#V+F;[s d]] of birth parents
+      //
+      template <
+        typename DerivedVA,
+        typename DerivedFA,
+        typename Deriveds,
+        typename Derivedd,
+        typename DerivedW,
+        typename DerivedG,
+        typename DerivedJ>
+      IGL_INLINE void minkowski_sum(
+        const Eigen::PlainObjectBase<DerivedVA> & VA,
+        const Eigen::PlainObjectBase<DerivedFA> & FA,
+        const Eigen::PlainObjectBase<Deriveds> & s,
+        const Eigen::PlainObjectBase<Derivedd> & d,
+        const bool resolve_overlaps,
+        Eigen::PlainObjectBase<DerivedW> & W,
+        Eigen::PlainObjectBase<DerivedG> & G,
+        Eigen::PlainObjectBase<DerivedJ> & J);
+      template <
+        typename DerivedVA,
+        typename DerivedFA,
+        typename Deriveds,
+        typename Derivedd,
+        typename DerivedW,
+        typename DerivedG,
+        typename DerivedJ>
+      IGL_INLINE void minkowski_sum(
+        const Eigen::PlainObjectBase<DerivedVA> & VA,
+        const Eigen::PlainObjectBase<DerivedFA> & FA,
+        const Eigen::PlainObjectBase<Deriveds> & s,
+        const Eigen::PlainObjectBase<Derivedd> & d,
+        Eigen::PlainObjectBase<DerivedW> & W,
+        Eigen::PlainObjectBase<DerivedG> & G,
+        Eigen::PlainObjectBase<DerivedJ> & J);
+
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "minkowski_sum.cpp"
+#endif
+
+#endif

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -413,6 +413,10 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
 
   remesh_intersections(V,F,T,offending,edge2faces,VV,FF,J,IM);
 
+#ifdef IGL_SELFINTERSECTMESH_DEBUG
+  cout<<"remesh intersection: "<<tictoc()<<endl;
+#endif
+
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.
 

include/igl/copyleft/cgal/closest_facet.cpp

@@ -226,21 +226,42 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
     return process_edge_case(query_idx, f[0], f[1], I(fid, 0), orientation);
   };
 
-  auto process_vertex_case = [&](const size_t query_idx, size_t s,
-      size_t preferred_facet, bool& orientation) {
+  // Given that the closest point to query point P(query_idx,:) on (V,F(I,:))
+  // is the vertex at V(s,:) which is incident at least on triangle
+  // F(preferred_facet,:), determine a facet incident on V(s,:) that is
+  // _exposed_ to the query point and determine whether that facet is facing
+  // _toward_ or _away_ from the query point.
+  //
+  // Inputs:
+  //   query_idx  index into P of query point
+  //   s  index into V of closest point at vertex
+  //   preferred_facet  facet incident on s
+  // Outputs:
+  //   orientation  whether returned face is facing toward or away from
+  //     query (parity unclear)
+  // Returns face guaranteed to be "exposed" to P(query_idx,:)
+  auto process_vertex_case = [&](
+    const size_t query_idx, 
+    size_t s,
+    size_t preferred_facet, 
+    bool& orientation) 
+  {
     const Point_3 query_point(
         P(query_idx, 0), P(query_idx, 1), P(query_idx, 2));
     const Point_3 closest_point(V(s, 0), V(s, 1), V(s, 2));
     std::vector<size_t> adj_faces;
     std::vector<size_t> adj_face_corners;
     {
-      // Gether adj faces to s within I.
+      // Gather adj faces to s within I.
       const auto& all_adj_faces = VF[s];
       const auto& all_adj_face_corners = VFi[s];
       const size_t num_all_adj_faces = all_adj_faces.size();
-      for (size_t i=0; i<num_all_adj_faces; i++) {
+      for (size_t i=0; i<num_all_adj_faces; i++) 
+      {
         const size_t fid = all_adj_faces[i];
-        if (in_I[fid]) {
+        // Shouldn't this always be true if I is a full connected component?
+        if (in_I[fid]) 
+        {
           adj_faces.push_back(fid);
           adj_face_corners.push_back(all_adj_face_corners[i]);
         }
@@ -251,7 +272,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
 
     std::set<size_t> adj_vertices_set;
     std::unordered_multimap<size_t, size_t> v2f;
-    for (size_t i=0; i<num_adj_faces; i++) {
+    for (size_t i=0; i<num_adj_faces; i++) 
+    {
       const size_t fid = adj_faces[i];
       const size_t cid = adj_face_corners[i];
       const auto& f = F.row(adj_faces[i]);
@@ -267,7 +289,8 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
         adj_vertices.begin());
 
     std::vector<Point_3> adj_points;
-    for (size_t vid : adj_vertices) {
+    for (size_t vid : adj_vertices) 
+    {
       adj_points.emplace_back(V(vid,0), V(vid,1), V(vid,2));
     }
 

include/igl/copyleft/cgal/extract_cells.cpp

@@ -11,6 +11,7 @@
 #include "../../facet_components.h"
 #include "../../triangle_triangle_adjacency.h"
 #include "../../unique_edge_map.h"
+#include "../../get_seconds.h"
 #include "closest_facet.h"
 #include "order_facets_around_edge.h"
 #include "outer_facet.h"
@@ -18,6 +19,8 @@
 #include <vector>
 #include <queue>
 
+//#define EXTRACT_CELLS_DEBUG
+
 template<
 typename DerivedV,
 typename DerivedF,
@@ -173,6 +176,16 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
         const std::vector<std::vector<uE2EType> >& uE2E,
         const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
         Eigen::PlainObjectBase<DerivedC>& cells) {
+#ifdef EXTRACT_CELLS_DEBUG
+  const auto & tictoc = []()
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  tictoc();
+#endif
     const size_t num_faces = F.rows();
     typedef typename DerivedF::Scalar Index;
     const size_t num_patches = P.maxCoeff()+1;
@@ -181,12 +194,21 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
     const size_t num_raw_cells =
         igl::copyleft::cgal::extract_cells_single_component(
                 V, F, P, uE, uE2E, EMAP, raw_cells);
+#ifdef EXTRACT_CELLS_DEBUG
+    std::cout << "Extract single component cells: " << tictoc() << std::endl;
+#endif
 
     std::vector<std::vector<std::vector<Index > > > TT,_1;
     igl::triangle_triangle_adjacency(E, EMAP, uE2E, false, TT, _1);
+#ifdef EXTRACT_CELLS_DEBUG
+    std::cout << "face adj: " << tictoc() << std::endl;
+#endif
 
     Eigen::VectorXi C, counts;
     igl::facet_components(TT, C, counts);
+#ifdef EXTRACT_CELLS_DEBUG
+    std::cout << "face comp: " << tictoc() << std::endl;
+#endif
 
     const size_t num_components = counts.size();
     std::vector<std::vector<size_t> > components(num_components);
@@ -209,6 +231,9 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
         outer_facet_orientation[i] = flipped?1:0;
         outer_cells[i] = raw_cells(P[outer_facets[i]], outer_facet_orientation[i]);
     }
+#ifdef EXTRACT_CELLS_DEBUG
+    std::cout << "Per comp outer facet: " << tictoc() << std::endl;
+#endif
 
     auto get_triangle_center = [&](const size_t fid) {
         return ((V.row(F(fid, 0)) + V.row(F(fid, 1)) + V.row(F(fid, 2)))
@@ -218,12 +243,47 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
     std::vector<std::vector<size_t> > ambient_cells(num_raw_cells);
     std::vector<std::vector<size_t> > ambient_comps(num_components);
     if (num_components > 1) {
+        DerivedV bbox_min(num_components, 3);
+        DerivedV bbox_max(num_components, 3);
+        bbox_min.rowwise() = V.colwise().maxCoeff().eval();
+        bbox_max.rowwise() = V.colwise().minCoeff().eval();
+        for (size_t i=0; i<num_faces; i++) {
+          const auto comp_id = C[i];
+          const auto& f = F.row(i);
+          for (size_t j=0; j<3; j++) {
+            bbox_min(comp_id, 0) = std::min(bbox_min(comp_id, 0), V(f[j], 0));
+            bbox_min(comp_id, 1) = std::min(bbox_min(comp_id, 1), V(f[j], 1));
+            bbox_min(comp_id, 2) = std::min(bbox_min(comp_id, 2), V(f[j], 2));
+            bbox_max(comp_id, 0) = std::max(bbox_max(comp_id, 0), V(f[j], 0));
+            bbox_max(comp_id, 1) = std::max(bbox_max(comp_id, 1), V(f[j], 1));
+            bbox_max(comp_id, 2) = std::max(bbox_max(comp_id, 2), V(f[j], 2));
+          }
+        }
+        auto bbox_intersects = [&](size_t comp_i, size_t comp_j) {
+          return !(
+              bbox_max(comp_i,0) < bbox_min(comp_j,0) ||
+              bbox_max(comp_i,1) < bbox_min(comp_j,1) ||
+              bbox_max(comp_i,2) < bbox_min(comp_j,2) ||
+              bbox_max(comp_j,0) < bbox_min(comp_i,0) ||
+              bbox_max(comp_j,1) < bbox_min(comp_i,1) ||
+              bbox_max(comp_j,2) < bbox_min(comp_i,2));
+        };
+
         for (size_t i=0; i<num_components; i++) {
-            DerivedV queries(num_components-1, 3);
+            std::vector<size_t> candidate_comps;
+            candidate_comps.reserve(num_components);
             for (size_t j=0; j<num_components; j++) {
                 if (i == j) continue;
-                size_t index = j<i ? j:j-1;
-                queries.row(index) = get_triangle_center(outer_facets[j]);
+                if (bbox_intersects(i,j)) candidate_comps.push_back(j);
+            }
+
+            const size_t num_candidate_comps = candidate_comps.size();
+            if (num_candidate_comps == 0) continue;
+
+            DerivedV queries(num_candidate_comps, 3);
+            for (size_t j=0; j<num_candidate_comps; j++) {
+                const size_t index = candidate_comps[j];
+                queries.row(j) = get_triangle_center(outer_facets[index]);
             }
 
             const auto& I = Is[i];
@@ -231,22 +291,24 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
             igl::copyleft::cgal::closest_facet(V, F, I, queries,
                     uE2E, EMAP, closest_facets, closest_facet_orientations);
 
-            for (size_t j=0; j<num_components; j++) {
-                if (i == j) continue;
-                size_t index = j<i ? j:j-1;
-                const size_t closest_patch = P[closest_facets[index]];
-                const size_t closest_patch_side = closest_facet_orientations[index]
+            for (size_t j=0; j<num_candidate_comps; j++) {
+                const size_t index = candidate_comps[j];
+                const size_t closest_patch = P[closest_facets[j]];
+                const size_t closest_patch_side = closest_facet_orientations[j]
                     ? 0:1;
                 const size_t ambient_cell = raw_cells(closest_patch,
                         closest_patch_side);
                 if (ambient_cell != (size_t)outer_cells[i]) {
-                    nested_cells[ambient_cell].push_back(outer_cells[j]);
-                    ambient_cells[outer_cells[j]].push_back(ambient_cell);
-                    ambient_comps[j].push_back(i);
+                    nested_cells[ambient_cell].push_back(outer_cells[index]);
+                    ambient_cells[outer_cells[index]].push_back(ambient_cell);
+                    ambient_comps[index].push_back(i);
                 }
             }
         }
     }
+#ifdef EXTRACT_CELLS_DEBUG
+    std::cout << "Determine nested relaitonship: " << tictoc() << std::endl;
+#endif
 
     const size_t INVALID = std::numeric_limits<size_t>::max();
     const size_t INFINITE_CELL = num_raw_cells;
@@ -308,6 +370,9 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
         raw_cells(i, 1) = negative_cell_id;
     }
     cells = raw_cells;
+#ifdef EXTRACT_CELLS_DEBUG
+    std::cout << "Finalize and output: " << tictoc() << std::endl;
+#endif
     return count;
 }
 

include/igl/copyleft/cgal/piecewise_constant_winding_number.cpp

@@ -0,0 +1,29 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "piecewise_constant_winding_number.h"
+#include "../../piecewise_constant_winding_number.h"
+#include "remesh_self_intersections.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <algorithm>
+
+template < typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::piecewise_constant_winding_number(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F)
+{
+  Eigen::Matrix<CGAL::Epeck::FT,Eigen::Dynamic,3> VV;
+  Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,3> FF;
+  Eigen::Matrix<typename DerivedF::Index,Eigen::Dynamic,2> IF;
+  Eigen::Matrix<typename DerivedF::Index,Eigen::Dynamic,1> J;
+  Eigen::Matrix<typename DerivedV::Index,Eigen::Dynamic,1> UIM,IM;
+  // resolve intersections
+  remesh_self_intersections(V,F,{},VV,FF,IF,J,IM);
+  // _apply_ duplicate vertex mapping IM to FF
+  for_each(FF.data(),FF.data()+FF.size(),[&IM](int & a){a=IM(a);});
+  return piecewise_constant_winding_number(FF);
+}

include/igl/copyleft/cgal/piecewise_constant_winding_number.h

@@ -0,0 +1,41 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_PIECEWISE_CONSTANT_WINDING_NUMBER_H
+#define IGL_COPYLEFT_CGAL_PIECEWISE_CONSTANT_WINDING_NUMBER_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // PIECEWISE_CONSTANT_WINDING_NUMBER Determine if a given mesh induces a
+      // piecewise constant winding number field: Is this mesh valid input to
+      // solid set operations.
+      // 
+      // Inputs:
+      //   V  #V by 3 list of mesh vertex positions
+      //   F  #F by 3 list of triangle indices into V
+      // Returns true if the mesh _combinatorially_ induces a piecewise
+      // constant winding number field.
+      template <
+        typename DerivedV,
+        typename DerivedF>
+      IGL_INLINE bool piecewise_constant_winding_number(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF>& F);
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "piecewise_constant_winding_number.cpp"
+#endif
+#endif
+

include/igl/copyleft/cgal/propagate_winding_numbers.cpp

@@ -11,8 +11,10 @@
 #include "../../extract_non_manifold_edge_curves.h"
 #include "../../facet_components.h"
 #include "../../unique_edge_map.h"
+#include "../../piecewise_constant_winding_number.h"
 #include "../../writeOBJ.h"
 #include "../../writePLY.h"
+#include "../../get_seconds.h"
 #include "order_facets_around_edge.h"
 #include "outer_facet.h"
 #include "closest_facet.h"
@@ -25,46 +27,7 @@
 #include <tuple>
 #include <queue>
 
-namespace propagate_winding_numbers_helper {
-  template<
-    typename DerivedF,
-    typename DeriveduE,
-    typename uE2EType >
-  bool is_orientable(
-      const Eigen::PlainObjectBase<DerivedF>& F,
-      const Eigen::PlainObjectBase<DeriveduE>& uE,
-      const std::vector<std::vector<uE2EType> >& uE2E) {
-    const size_t num_faces = F.rows();
-    const size_t num_edges = uE.rows();
-    auto edge_index_to_face_index = [&](size_t ei) {
-      return ei % num_faces;
-    };
-    auto is_consistent = [&](size_t fid, size_t s, size_t d) {
-      if ((size_t)F(fid, 0) == s && (size_t)F(fid, 1) == d) return true;
-      if ((size_t)F(fid, 1) == s && (size_t)F(fid, 2) == d) return true;
-      if ((size_t)F(fid, 2) == s && (size_t)F(fid, 0) == d) return true;
-
-      if ((size_t)F(fid, 0) == d && (size_t)F(fid, 1) == s) return false;
-      if ((size_t)F(fid, 1) == d && (size_t)F(fid, 2) == s) return false;
-      if ((size_t)F(fid, 2) == d && (size_t)F(fid, 0) == s) return false;
-      throw "Invalid face!!";
-    };
-    for (size_t i=0; i<num_edges; i++) {
-      const size_t s = uE(i,0);
-      const size_t d = uE(i,1);
-      int count=0;
-      for (const auto& ei : uE2E[i]) {
-        const size_t fid = edge_index_to_face_index(ei);
-        if (is_consistent(fid, s, d)) count++;
-        else count--;
-      }
-      if (count != 0) {
-        return false;
-      }
-    }
-    return true;
-  }
-}
+//#define PROPAGATE_WINDING_NUMBER_TIMING
 
 template<
   typename DerivedV,
@@ -76,6 +39,16 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     const Eigen::PlainObjectBase<DerivedF>& F,
     const Eigen::PlainObjectBase<DerivedL>& labels,
     Eigen::PlainObjectBase<DerivedW>& W) {
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  const auto & tictoc = []()
+  {
+    static double t_start = igl::get_seconds();
+    double diff = igl::get_seconds()-t_start;
+    t_start += diff;
+    return diff;
+  };
+  tictoc();
+#endif
   const size_t num_faces = F.rows();
   //typedef typename DerivedF::Scalar Index;
 
@@ -83,17 +56,24 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
   Eigen::VectorXi EMAP;
   std::vector<std::vector<size_t> > uE2E;
   igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-  if (!propagate_winding_numbers_helper::is_orientable(F, uE, uE2E)) {
-      std::cerr << "Input mesh is not orientable!" << std::endl;
+  if (!piecewise_constant_winding_number(F, uE, uE2E)) 
+  {
+    std::cerr << "Input mesh is not orientable!" << std::endl;
   }
 
   Eigen::VectorXi P;
   const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  std::cout << "extract manifold patches: " << tictoc() << std::endl;
+#endif
 
   DerivedW per_patch_cells;
   const size_t num_cells =
     igl::copyleft::cgal::extract_cells(
         V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  std::cout << "extract cells: " << tictoc() << std::endl;;
+#endif
 
   typedef std::tuple<size_t, bool, size_t> CellConnection;
   std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
@@ -103,6 +83,9 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
     cell_adjacency[positive_cell].emplace(negative_cell, false, i);
     cell_adjacency[negative_cell].emplace(positive_cell, true, i);
   }
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  std::cout << "cell connection: " << tictoc() << std::endl;
+#endif
 
   auto save_cell = [&](const std::string& filename, size_t cell_id) {
     std::vector<size_t> faces;
@@ -175,12 +158,19 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
                 }
                 std::cout << std::endl;
               }
-              assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
+              // Do not fail when odd cycle is detected because the resulting
+              // integer winding number field, although inconsistent, may still
+              // be used if the problem region is local and embedded within a
+              // valid volume.
+              //assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
             }
           }
         }
       }
     }
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+    std::cout << "check for odd cycle: " << tictoc() << std::endl;
+#endif
   }
 #endif
 
@@ -189,6 +179,9 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
   Eigen::VectorXi I;
   I.setLinSpaced(num_faces, 0, num_faces-1);
   igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  std::cout << "outer facet: " << tictoc() << std::endl;
+#endif
 
   const size_t outer_patch = P[outer_facet];
   const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
@@ -229,20 +222,30 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
         Q.push(neighbor_cell);
       } else {
 #ifndef NDEBUG
+        // Checking for winding number consistency.
+        // This check would inevitably fail for meshes that contain open
+        // boundary or non-orientable.  However, the inconsistent winding number
+        // field would still be useful in some cases such as when problem region
+        // is local and embedded within the volume.  This, unfortunately, is the
+        // best we can do because the problem of computing integer winding
+        // number is ill-defined for open and non-orientable surfaces.
         for (size_t i=0; i<num_labels; i++) {
           if ((int)i == patch_labels[patch_idx]) {
             int inc = direction ? -1:1;
-            assert(per_cell_W(neighbor_cell, i) ==
-                per_cell_W(curr_cell, i) + inc);
+            //assert(per_cell_W(neighbor_cell, i) ==
+            //    per_cell_W(curr_cell, i) + inc);
           } else {
-            assert(per_cell_W(neighbor_cell, i) ==
-                per_cell_W(curr_cell, i));
+            //assert(per_cell_W(neighbor_cell, i) ==
+            //    per_cell_W(curr_cell, i));
           }
         }
 #endif
       }
     }
   }
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  std::cout << "propagate winding number: " << tictoc() << std::endl;
+#endif
 
   W.resize(num_faces, num_labels*2);
   for (size_t i=0; i<num_faces; i++) {
@@ -254,6 +257,9 @@ IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
       W(i,j*2+1) = per_cell_W(negative_cell, j);
     }
   }
+#ifdef PROPAGATE_WINDING_NUMBER_TIMING
+  std::cout << "save result: " << tictoc() << std::endl;
+#endif
 }
 
 

include/igl/cross_field_missmatch.cpp

@@ -28,6 +28,7 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
   private:
@@ -35,6 +36,7 @@ namespace igl {
     std::vector<bool> V_border; // bool
     std::vector<std::vector<int> > VF;
     std::vector<std::vector<int> > VFi;
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 
@@ -88,7 +90,7 @@ public:
     igl::per_face_normals(V,F,N);
     V_border = igl::is_border_vertex(V,F);
     igl::vertex_triangle_adjacency(V,F,VF,VFi);
-    igl::triangle_triangle_adjacency(V,F,TT,TTi);
+    igl::triangle_triangle_adjacency(F,TT,TTi);
   }
 
   inline void calculateMissmatch(Eigen::PlainObjectBase<DerivedF> &Handle_MMatch)
@@ -116,8 +118,8 @@ IGL_INLINE void igl::cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV
                                            const bool isCombed,
                                            Eigen::PlainObjectBase<DerivedF> &missmatch)
 {
-  Eigen::PlainObjectBase<DerivedV> PD1_combed;
-  Eigen::PlainObjectBase<DerivedV> PD2_combed;
+  DerivedV PD1_combed;
+  DerivedV PD2_combed;
 
   if (!isCombed)
     igl::comb_cross_field(V,F,PD1,PD2,PD1_combed,PD2_combed);

include/igl/cut_mesh.cpp

@@ -29,6 +29,7 @@ namespace igl {
     int num_scalar_variables;
 
     // per face indexes of vertex in the solver
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> HandleS_Index;
 
     // per vertex variable indexes
@@ -308,7 +309,7 @@ IGL_INLINE void igl::cut_mesh(
   Eigen::MatrixXd Vt = V;
   Eigen::MatrixXi Ft = F;
   Eigen::MatrixXi TT, TTi;
-  igl::triangle_triangle_adjacency(Vt,Ft,TT,TTi);
+  igl::triangle_triangle_adjacency(Ft,TT,TTi);
 
   std::vector<bool> V_border = igl::is_border_vertex(V,F);
 

include/igl/cut_mesh_from_singularities.cpp

@@ -30,6 +30,7 @@ namespace igl {
     const Eigen::PlainObjectBase<DerivedM> &Handle_MMatch;
 
     Eigen::VectorXi F_visited;
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedF> TT;
     Eigen::PlainObjectBase<DerivedF> TTi;
 
@@ -147,7 +148,7 @@ namespace igl {
     F(F_),
     Handle_MMatch(Handle_MMatch_)
     {
-      triangle_triangle_adjacency(V,F,TT,TTi);
+      triangle_triangle_adjacency(F,TT,TTi);
       edge_topology(V,F,E,F2E,E2F);
     };
 

include/igl/doublearea.cpp

@@ -25,7 +25,7 @@ IGL_INLINE void igl::doublearea(
   assert(F.cols() == 3);
   const size_t m = F.rows();
   // Compute edge lengths
-  Eigen::PlainObjectBase<DerivedV> l;
+  Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3> l;
   // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
   // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
 
@@ -140,7 +140,7 @@ IGL_INLINE void igl::doublearea(
   assert(ul.cols() == 3);
   // Number of triangles
   const Index m = ul.rows();
-  Eigen::PlainObjectBase<Derivedl> l;
+  Eigen::Matrix<typename Derivedl::Scalar, Eigen::Dynamic, 3> l;
   MatrixXi _;
   sort(ul,2,false,l,_);
   // semiperimeters

include/igl/embree/EmbreeIntersector.h

@@ -16,7 +16,7 @@
 #ifndef IGL_EMBREE_EMBREE_INTERSECTOR_H
 #define IGL_EMBREE_EMBREE_INTERSECTOR_H
 
-#include "Hit.h"
+#include "../Hit.h"
 #include <Eigen/Geometry>
 #include <Eigen/Core>
 #include <Eigen/Geometry>
@@ -57,11 +57,13 @@ namespace igl
       // Inputs:
       //   V  #V by 3 list of vertex positions
       //   F  #F by 3 list of Oriented triangles
+      //   isStatic  scene is optimized for static geometry
       // Side effects:
       //   The first time this is ever called the embree engine is initialized.
       inline void init(
         const PointMatrixType& V,
-        const FaceMatrixType& F);
+        const FaceMatrixType& F,
+        bool isStatic = false);
 
       // Initialize with a given mesh.
       //
@@ -69,12 +71,14 @@ namespace igl
       //   V  vector of #V by 3 list of vertex positions for each geometry
       //   F  vector of #F by 3 list of Oriented triangles for each geometry
       //   masks  a 32 bit mask to identify active geometries.
+      //   isStatic  scene is optimized for static geometry
       // Side effects:
       //   The first time this is ever called the embree engine is initialized.
       inline void init(
         const std::vector<const PointMatrixType*>& V,
         const std::vector<const FaceMatrixType*>& F,
-        const std::vector<int>& masks);
+        const std::vector<int>& masks,
+        bool isStatic = false);
 
       // Deinitialize embree datasctructures for current mesh.  Also called on
       // destruction: no need to call if you just want to init() once and
@@ -252,7 +256,8 @@ inline igl::embree::EmbreeIntersector & igl::embree::EmbreeIntersector::operator
 
 inline void igl::embree::EmbreeIntersector::init(
   const PointMatrixType& V,
-  const FaceMatrixType& F)
+  const FaceMatrixType& F,
+  bool isStatic)
 {
   std::vector<const PointMatrixType*> Vtemp;
   std::vector<const FaceMatrixType*> Ftemp;
@@ -260,13 +265,14 @@ inline void igl::embree::EmbreeIntersector::init(
   Vtemp.push_back(&V);
   Ftemp.push_back(&F);
   masks.push_back(0xFFFFFFFF);
-  init(Vtemp,Ftemp,masks);
+  init(Vtemp,Ftemp,masks,isStatic);
 }
 
 inline void igl::embree::EmbreeIntersector::init(
   const std::vector<const PointMatrixType*>& V,
   const std::vector<const FaceMatrixType*>& F,
-  const std::vector<int>& masks)
+  const std::vector<int>& masks,
+  bool isStatic)
 {
 
   if(initialized)
@@ -282,7 +288,10 @@ inline void igl::embree::EmbreeIntersector::init(
   }
 
   // create a scene
-  scene = rtcNewScene(RTC_SCENE_ROBUST | RTC_SCENE_HIGH_QUALITY,RTC_INTERSECT1);
+  RTCSceneFlags flags = RTC_SCENE_ROBUST | RTC_SCENE_HIGH_QUALITY;
+  if(isStatic)
+    flags = flags | RTC_SCENE_STATIC;
+  scene = rtcNewScene(flags,RTC_INTERSECT1);
 
   for(int g=0;g<(int)V.size();g++)
   {

include/igl/embree/Hit.h

@@ -1,26 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-//               2014 Christian Schüller <schuellchr@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_EMBREE_HIT_H
-#define IGL_EMBREE_HIT_H
-
-namespace igl
-{
-  namespace embree
-  {
-    // Reimplementation of the embree::Hit struct from embree1.0
-    struct Hit
-    {
-      int id; // primitive id
-      int gid; // geometry id
-      float u,v; // barycentric coordinates
-      float t; // distance = direction*t to intersection
-    };
-  }
-}
-#endif 

include/igl/embree/ambient_occlusion.cpp

@@ -6,9 +6,9 @@
 // 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 "ambient_occlusion.h"
+#include "../ambient_occlusion.h"
 #include "EmbreeIntersector.h"
-#include <igl/random_dir.h>
-#include <igl/EPS.h>
+#include "../Hit.h"
 
 template <
   typename DerivedP,
@@ -21,38 +21,15 @@ IGL_INLINE void igl::embree::ambient_occlusion(
   const int num_samples,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
-  using namespace Eigen;
-  using namespace igl;
-  const int n = P.rows();
-  // Resize output
-  S.resize(n,1);
-  // Embree seems to be parallel when constructing but not when tracing rays
-#pragma omp parallel for
-  // loop over mesh vertices
-  for(int p = 0;p<n;p++)
+  const auto & shoot_ray = [&ei](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir)->bool
   {
-    const Vector3f origin = P.row(p).template cast<float>();
-    const Vector3f normal = N.row(p).template cast<float>();
-    int num_hits = 0;
-    MatrixXf D = random_dir_stratified(num_samples).cast<float>();
-    for(int s = 0;s<num_samples;s++)
-    {
-      //Vector3d d = random_dir();
-      Vector3f d = D.row(s);
-      if(d.dot(normal) < 0)
-      {
-        // reverse ray
-        d *= -1;
-      }
-      igl::embree::Hit hit;
-      const float tnear = 1e-4f;
-      if(ei.intersectRay(origin,d,hit,tnear))
-      {
-        num_hits++;
-      }
-    }
-    S(p) = (double)num_hits/(double)num_samples;
-  }
+    igl::Hit hit;
+    const float tnear = 1e-4f;
+    return ei.intersectRay(s,dir,hit,tnear);
+  };
+  return igl::ambient_occlusion(shoot_ray,P,N,num_samples,S);
 }
 
 template <

include/igl/embree/bone_visible.cpp

@@ -6,9 +6,10 @@
 // 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 "bone_visible.h"
-#include <igl/project_to_line.h>
-#include <igl/EPS.h>
-#include <igl/Timer.h>
+#include "../project_to_line.h"
+#include "../EPS.h"
+#include "../Hit.h"
+#include "../Timer.h"
 #include <iostream>
 
 template <
@@ -86,7 +87,7 @@ IGL_INLINE void igl::embree::bone_visible(
         projv = d;
       }
     }
-    igl::embree::Hit hit;
+    igl::Hit hit;
     // perhaps 1.0 should be 1.0-epsilon, or actually since we checking the
     // incident face, perhaps 1.0 should be 1.0+eps
     const Vector3d dir = (Vv-projv)*1.0;

include/igl/embree/line_mesh_intersection.cpp

@@ -6,6 +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 "line_mesh_intersection.h"
+#include "../Hit.h"
 
 // For error printing
 #include <cstdio>
@@ -40,7 +41,7 @@ IGL_INLINE ScalarMatrix igl::embree::line_mesh_intersection
   // Shoot rays from the source to the target
   for (unsigned i=0; i<ray_pos.rows(); ++i)
   {
-    igl::embree::Hit A,B;
+    igl::Hit A,B;
     // Shoot ray A
     Eigen::RowVector3d A_pos = ray_pos.row(i) + tol * ray_dir.row(i);
     Eigen::RowVector3d A_dir = -ray_dir.row(i);

include/igl/embree/unproject_in_mesh.cpp

@@ -7,7 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "unproject_in_mesh.h"
 #include "EmbreeIntersector.h"
-#include "../unproject.h"
+#include "../unproject_ray.h"
+#include "../unproject_in_mesh.h"
 #include <vector>
 
 template <typename Derivedobj>
@@ -18,41 +19,20 @@ IGL_INLINE int igl::embree::unproject_in_mesh(
   const Eigen::Vector4f& viewport,
   const EmbreeIntersector & ei,
   Eigen::PlainObjectBase<Derivedobj> & obj,
-  std::vector<igl::embree::Hit > & hits)
+  std::vector<igl::Hit > & hits)
 {
   using namespace igl;
   using namespace std;
   using namespace Eigen;
-  // Source and direction on screen
-  Vector3f win_s(pos(0),pos(1),0);
-  Vector3f win_d(pos(0),pos(1),1);
-  // Source, destination and direction in world
-  Vector3f s,d,dir;
-  s = igl::unproject(win_s,model,proj,viewport);
-  d = igl::unproject(win_d,model,proj,viewport);
-  dir = d-s;
-
-  // Shoot ray, collect all hits (could just collect first two)
-  int num_rays_shot;
-  hits.clear();
-  ei.intersectRay(s,dir,hits,num_rays_shot);
-  switch(hits.size())
+  const auto & shoot_ray = [&ei](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir,
+    std::vector<igl::Hit> & hits)
   {
-    case 0:
-      break;
-    case 1:
-    {
-      obj = (s + dir*hits[0].t).cast<typename Derivedobj::Scalar>();
-      break;
-    }
-    case 2:
-    default:
-    {
-      obj = 0.5*((s + dir*hits[0].t) + (s + dir*hits[1].t)).cast<typename Derivedobj::Scalar>();
-      break;
-    }
-  }
-  return hits.size();
+    int num_rays_shot;
+    ei.intersectRay(s,dir,hits,num_rays_shot);
+  };
+  return igl::unproject_in_mesh(pos,model,proj,viewport,shoot_ray,obj,hits);
 }
 
 template <typename Derivedobj>
@@ -64,13 +44,13 @@ IGL_INLINE int igl::embree::unproject_in_mesh(
     const EmbreeIntersector & ei,
     Eigen::PlainObjectBase<Derivedobj> & obj)
 {
-  std::vector<igl::embree::Hit> hits;
+  std::vector<igl::Hit> hits;
   return unproject_in_mesh(pos,model,proj,viewport,ei,obj,hits);
 }
 
 
 #ifdef IGL_STATIC_LIBRARY
-template int igl::embree::unproject_in_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<igl::embree::Hit, std::allocator<igl::embree::Hit> >&);
-template int igl::embree::unproject_in_mesh<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, std::vector<igl::embree::Hit, std::allocator<igl::embree::Hit> >&);
+template int igl::embree::unproject_in_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+template int igl::embree::unproject_in_mesh<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
 template int igl::embree::unproject_in_mesh<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
 #endif

include/igl/embree/unproject_in_mesh.h

@@ -11,7 +11,7 @@
 #include <Eigen/Core>
 
 #include <vector>
-#include "Hit.h"
+#include "../Hit.h"
 
 namespace igl
 {
@@ -54,7 +54,7 @@ namespace igl
       const Eigen::Vector4f& viewport,
       const EmbreeIntersector & ei,
       Eigen::PlainObjectBase<Derivedobj> & obj,
-      std::vector<igl::embree::Hit > & hits);
+      std::vector<igl::Hit > & hits);
     template < typename Derivedobj>
     IGL_INLINE int unproject_in_mesh(
       const Eigen::Vector2f& pos,

include/igl/embree/unproject_onto_mesh.cpp

@@ -24,7 +24,7 @@ IGL_INLINE bool igl::embree::unproject_onto_mesh(
   using namespace std;
   using namespace Eigen;
   MatrixXd obj;
-  vector<igl::embree::Hit> hits;
+  vector<igl::Hit> hits;
 
   // This is lazy, it will find more than just the first hit
   unproject_in_mesh(pos,model,proj,viewport,ei,obj,hits);

include/igl/embree/unproject_onto_mesh.h

@@ -11,7 +11,6 @@
 #include <Eigen/Core>
 
 #include <vector>
-#include "Hit.h"
 
 namespace igl
 {

include/igl/find_cross_field_singularities.cpp

@@ -71,7 +71,7 @@ IGL_INLINE void igl::find_cross_field_singularities(const Eigen::PlainObjectBase
   // Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > Handle_MMatch;
   // Eigen::PlainObjectBase<Eigen::Matrix<int,  Eigen::Dynamic, 3> > Handle_MMatch;
   // Eigen::Matrix<int, Eigen::Dynamic, 3> Handle_MMatch;
-  Eigen::PlainObjectBase<DerivedF> Handle_MMatch;
+  Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, 3> Handle_MMatch;
   // Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, 3> Handle_MMatch;
   // Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, Eigen::Dynamic> Handle_MMatch;
   igl::cross_field_missmatch<DerivedV, DerivedF>(V, F, PD1, PD2, isCombed, Handle_MMatch);

include/igl/fit_rigid.cpp

@@ -1,21 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#include "fit_rigid.h"
-#include "procrustes.h"
-
-IGL_INLINE void igl::fit_rigid(
-  const Eigen::MatrixXd & A,
-  const Eigen::MatrixXd & B,
-  Eigen::Rotation2Dd & R,
-  Eigen::RowVector2d & t)
-{
-  using namespace Eigen;
-  Matrix2d Rmat;
-  procrustes(A,B,Rmat,t);
-  R.fromRotationMatrix(Rmat);
-}

include/igl/fit_rigid.h

@@ -1,37 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-#ifndef IGL_FIT_RIGID_H
-#define IGL_FIT_RIGID_H
-
-#if defined(_WIN32)
-  #pragma message("Deprecated. Use igl/procrustes.h instead")
-#else
-  #warning "Deprecated. Use igl/procrustes.h instead"
-#endif
-
-#include "igl_inline.h"
-#include <Eigen/Core>
-#include <Eigen/Geometry>
-
-namespace igl
-{
-  // Deprecated: please use procrustes(...) instead.
-  // Fit a rigid 
-  IGL_INLINE void fit_rigid(
-    const Eigen::MatrixXd & A,
-    const Eigen::MatrixXd & B,
-    Eigen::Rotation2Dd & R,
-    Eigen::RowVector2d & t);
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "fit_rigid.cpp"
-#endif
-
-#endif
-

include/igl/full.cpp

@@ -1,44 +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 "full.h"
-
-template <typename T,typename DerivedB>
-IGL_INLINE void igl::full(
-  const Eigen::SparseMatrix<T> & A,
-  Eigen::PlainObjectBase<DerivedB>& B)
-{
-  assert(false && "Obsolete. Just call B = Matrix(A)");
-  using namespace Eigen;
-  B = PlainObjectBase<DerivedB >::Zero(A.rows(),A.cols());
-  // Iterate over outside
-  for(int k=0; k<A.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename SparseMatrix<T>::InnerIterator it (A,k); it; ++it)
-    {
-      B(it.row(),it.col()) = it.value();
-    }
-  }
-}
-
-template <typename DerivedA,typename DerivedB>
-IGL_INLINE void igl::full(
-  const Eigen::PlainObjectBase<DerivedA>& A,
-  Eigen::PlainObjectBase<DerivedB>& B)
-{
-  assert(false && "Obsolete. Just call B = Matrix(A)");
-  B = A;
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-// generated by autoexplicit.sh
-template void igl::full<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::full<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-#endif

include/igl/full.h

@@ -1,41 +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_FULL_H
-#define IGL_FULL_H
-#include "igl_inline.h"
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-namespace igl
-{
-  // This is totally unnecessary. You can just call MatrixXd B = MatrixXd(A);
-  //
-  // Convert a sparsematrix into a full one
-  //
-  // Templates:
-  //   T  should be a eigen sparse matrix primitive type like int or double
-  // Input:
-  //   A  m by n sparse matrix
-  // Output:
-  //   B  m by n dense/full matrix
-  template <typename T,typename DerivedB>
-  IGL_INLINE void full(
-    const Eigen::SparseMatrix<T> & A,
-    Eigen::PlainObjectBase<DerivedB> & B);
-  // If already full then this will just be a copy by assignment
-  template <typename DerivedA,typename DerivedB>
-  IGL_INLINE void full(
-    const Eigen::PlainObjectBase<DerivedA>& A,
-    Eigen::PlainObjectBase<DerivedB>& B);
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "full.cpp"
-#endif
-
-#endif

include/igl/grad.cpp

@@ -14,9 +14,8 @@ IGL_INLINE void igl::grad(const Eigen::PlainObjectBase<DerivedV>&V,
                      const Eigen::PlainObjectBase<DerivedF>&F,
                     Eigen::SparseMatrix<typename DerivedV::Scalar> &G)
 {
-  Eigen::PlainObjectBase<DerivedV > eperp21, eperp13;
-  eperp21.resize(F.rows(),3);
-  eperp13.resize(F.rows(),3);
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> 
+    eperp21(F.rows(),3), eperp13(F.rows(),3);
 
   for (int i=0;i<F.rows();++i)
   {

include/igl/is_border_vertex.cpp

@@ -12,11 +12,11 @@
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE std::vector<bool> igl::is_border_vertex(
-    const Eigen::PlainObjectBase<DerivedV> &V, 
+    const Eigen::PlainObjectBase<DerivedV> &V,
     const Eigen::PlainObjectBase<DerivedF> &F)
 {
-  Eigen::PlainObjectBase<DerivedF> FF;
-  igl::triangle_triangle_adjacency(V,F,FF);
+  DerivedF FF;
+  igl::triangle_triangle_adjacency(F,FF);
   std::vector<bool> ret(V.rows());
   for(unsigned i=0; i<ret.size();++i)
     ret[i] = false;
@@ -35,4 +35,5 @@ IGL_INLINE std::vector<bool> igl::is_border_vertex(
 // Explicit template specialization
 template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&);
 template std::vector<bool, std::allocator<bool> > igl::is_border_vertex<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 std::vector<bool, std::allocator<bool> > igl::is_border_vertex<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 #endif

include/igl/lim/lim.cpp

@@ -8,37 +8,39 @@
 #include "lim.h"
 #include <LIMSolverInterface.h>
 
-IGL_INLINE int igl::lim::lim(
+using namespace igl::lim;
+
+IGL_INLINE State lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
   const Eigen::SparseMatrix<double>& constraintMatrix,
   const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-  int energyType,
+  Energy energyType,
   double tolerance,
   int maxIteration,
   bool findLocalMinima)
 {
-  return ComputeLIM(
+  return (State)ComputeLIM(
     vertices,
     initialVertices,
     elements,
     constraintMatrix,
     constraintTargets,
-    energyType,
+    (EnergyType)energyType,
     tolerance,
     maxIteration,
     findLocalMinima
     );
 }
 
-IGL_INLINE int igl::lim(
+IGL_INLINE State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
   const Eigen::SparseMatrix<double>& constraintMatrix,
   const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-  int energyType,
+  Energy energyType,
   double tolerance,
   int maxIteration,
   bool findLocalMinima,
@@ -48,13 +50,13 @@ IGL_INLINE int igl::lim(
   double beta,
   double eps)
 {
-  return ComputeLIM(
+  return (State)ComputeLIM(
     vertices,
     initialVertices,
     elements,
     constraintMatrix,
     constraintTargets,
-    energyType,
+    (EnergyType)energyType,
     tolerance,
     maxIteration,
     findLocalMinima,
@@ -66,7 +68,7 @@ IGL_INLINE int igl::lim(
     );
 }
 
-IGL_INLINE int igl::lim(
+IGL_INLINE State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -74,12 +76,12 @@ IGL_INLINE int igl::lim(
   const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
   const Eigen::SparseMatrix<double>& constraintMatrix,
   const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-  int energyType,
+  Energy energyType,
   double tolerance,
   int maxIteration,
   bool findLocalMinima)
 {
-  return ComputeLIM(
+  return (State)ComputeLIM(
     vertices,
     initialVertices,
     elements,
@@ -87,14 +89,14 @@ IGL_INLINE int igl::lim(
     gradients,
     constraintMatrix,
     constraintTargets,
-    energyType,
+    (EnergyType)energyType,
     tolerance,
     maxIteration,
     findLocalMinima
     );
 }
 
-IGL_INLINE int igl::lim(
+IGL_INLINE State igl::lim::lim(
   Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
   const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
   const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -102,7 +104,7 @@ IGL_INLINE int igl::lim(
   const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
   const Eigen::SparseMatrix<double>& constraintMatrix,
   const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-  int energyType,
+  Energy energyType,
   double tolerance,
   int maxIteration,
   bool findLocalMinima,
@@ -112,7 +114,7 @@ IGL_INLINE int igl::lim(
   double beta,
   double eps)
 {
-  return ComputeLIM(
+  return (State)ComputeLIM(
     vertices,
     initialVertices,
     elements,
@@ -120,7 +122,7 @@ IGL_INLINE int igl::lim(
     gradients,
     constraintMatrix,
     constraintTargets,
-    energyType,
+    (EnergyType)energyType,
     tolerance,
     maxIteration,
     findLocalMinima,

include/igl/lim/lim.h

@@ -15,9 +15,6 @@ namespace igl
 {
   namespace lim
   {
-    // Known issues: energy type should be a readable enum rather than magic
-    // ints.
-    //
     // Computes a locally injective mapping of a triangle or tet-mesh based on
     // a deformation energy subject to some provided linear positional
     // constraints Cv-d.
@@ -39,7 +36,7 @@ namespace igl
     //                     y_2, ..., x_v,y_v])
     //   constraintTargets d: c vector target positions
     //   energyType        type of used energy:
-    //                     0=Dirichlet,1=Laplacian,2=Green,3=ARAP,4=LSCM
+    //                     Dirichlet, Laplacian, Green, ARAP, LSCM, Poisson (only 2D), UniformLaplacian, Identity
     //   tolerance         max squared positional constraints error
     //   maxIteration      max number of iterations
     //   findLocalMinima   iterating until a local minima is found. If not
@@ -61,30 +58,32 @@ namespace igl
     //                   mesh
     //--------------------------------------------------------------------------
     // Return values:
-    //  1 : Optimization deemed successful because either (a) it stagnated
-    //    (very step size) or (b) positional constraints were satisfied. (re:
-    //    https://github.com/libigl/libigl/issues/79 )
-    // -1 : Max iteration reached before tolerance was fulfilled
-    // -2 : not feasible -> has inverted elements (may want to decrease eps?)
+    //  Succeeded : Successful optimization with fulfilled tolerance
+    //  LocalMinima : Convergenged to a local minima / tolerance not fullfilled
+    //  IterationLimit : Max iteration reached before tolerance was fulfilled
+    //  Infeasible : not feasible -> has inverted elements (decrease eps?)
   
-    int lim(
+    enum Energy { Dirichlet = 0, Laplacian=1, Green=2, ARAP=3, LSCM=4, Poisson=5, UniformLaplacian=6, Identity=7 };
+    enum State { Uninitialized = -4, Infeasible = -3, IterationLimit = -2, LocalMinima = -1, Running = 0, Succeeded = 1 };
+
+    State lim(
       Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
       const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
       const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
       const Eigen::SparseMatrix<double>& constraintMatrix,
       const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-      int energyType,
+      Energy energyType,
       double tolerance,
       int maxIteration,
       bool findLocalMinima);
   
-    int lim(
+    State lim(
       Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
       const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
       const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
       const Eigen::SparseMatrix<double>& constraintMatrix,
       const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-      int energyType,
+      Energy energyType,
       double tolerance,
       int maxIteration,
       bool findLocalMinima,
@@ -94,7 +93,7 @@ namespace igl
       double beta,
       double eps);
   
-    int lim(
+    State lim(
       Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
       const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
       const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -102,12 +101,12 @@ namespace igl
       const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
       const Eigen::SparseMatrix<double>& constraintMatrix,
       const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-      int energyType,
+      Energy energyType,
       double tolerance,
       int maxIteration,
       bool findLocalMinima);
   
-    int lim(
+    State lim(
       Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
       const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
       const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
@@ -115,7 +114,7 @@ namespace igl
       const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
       const Eigen::SparseMatrix<double>& constraintMatrix,
       const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
-      int energyType,
+      Energy energyType,
       double tolerance,
       int maxIteration,
       bool findLocalMinima,

include/igl/line_field_missmatch.cpp

@@ -31,6 +31,7 @@ public:
     const Eigen::PlainObjectBase<DerivedF> &F;
     const Eigen::PlainObjectBase<DerivedV> &PD1;
     const Eigen::PlainObjectBase<DerivedV> &PD2;
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> N;
 
 private:
@@ -93,7 +94,7 @@ public:
         igl::per_face_normals(V,F,N);
         V_border = igl::is_border_vertex(V,F);
         igl::vertex_triangle_adjacency(V,F,VF,VFi);
-        igl::triangle_triangle_adjacency(V,F,TT,TTi);
+        igl::triangle_triangle_adjacency(F,TT,TTi);
     }
 
     inline void calculateMissmatchLine(Eigen::PlainObjectBase<DerivedO> &Handle_MMatch)

include/igl/list_to_matrix.cpp

@@ -15,8 +15,8 @@
 #include "max_size.h"
 #include "min_size.h"
 
-template <typename T, class Mat>
-IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Mat & M)
+template <typename T, typename Derived>
+IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Eigen::PlainObjectBase<Derived>& M)
 {
   // number of rows
   int m = V.size();
@@ -48,12 +48,12 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Mat
   return true;
 }
 
-template <typename T, class Mat>
+template <typename T, typename Derived>
 IGL_INLINE bool igl::list_to_matrix(
   const std::vector<std::vector<T > > & V,
   const int n,
   const T & padding,
-  Mat & M)
+  Eigen::PlainObjectBase<Derived>& M)
 {
   const int m = V.size();
   M.resize(m,n);
@@ -77,8 +77,8 @@ IGL_INLINE bool igl::list_to_matrix(
   return true;
 }
 
-template <typename T, class Mat>
-IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Mat & M)
+template <typename T, typename Derived>
+IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Eigen::PlainObjectBase<Derived>& M)
 {
   // number of rows
   int m = V.size();
@@ -86,10 +86,10 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Mat & M)
   {
     //fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
     //return false;
-    if(Mat::ColsAtCompileTime == 1)
+    if(Eigen::PlainObjectBase<Derived>::ColsAtCompileTime == 1)
     {
       M.resize(0,1);
-    }else if(Mat::RowsAtCompileTime == 1)
+    }else if(Eigen::PlainObjectBase<Derived>::RowsAtCompileTime == 1)
     {
       M.resize(1,0);
     }else
@@ -99,7 +99,7 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Mat & M)
     return true;
   }
   // Resize output
-  if(Mat::RowsAtCompileTime == 1)
+  if(Eigen::PlainObjectBase<Derived>::RowsAtCompileTime == 1)
   {
     M.resize(1,m);
   }else
@@ -118,50 +118,21 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Mat & M)
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<unsigned int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > >(std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<float, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::vector<double, std::allocator<double> > const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-// generated by autoexplicit.sh
-template bool igl::list_to_matrix<double, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1>&);
-template bool igl::list_to_matrix<bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template bool igl::list_to_matrix<bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > >(std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template bool igl::list_to_matrix<int, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template bool igl::list_to_matrix<int, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<int, std::allocator<int> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2> >&);
-template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >&);
-template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<int, std::allocator<int> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
-template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
-template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template bool igl::list_to_matrix<double, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>&);
-template bool igl::list_to_matrix<unsigned long, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<unsigned long, std::allocator<unsigned long> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template bool igl::list_to_matrix<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<double, std::allocator<double> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
-template bool igl::list_to_matrix<float, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
-template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::Matrix<double, -1, 2, 0, -1, 2>&);
-template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::Matrix<double, -1, 3, 0, -1, 3>&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> > >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 2, 1, 0, 2, 1> >&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> > >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 4, 1, 0, 4, 1> >&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > >(std::vector<double, std::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
-template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::Matrix<int, -1, 2, 0, -1, 2>&);
-template bool igl::list_to_matrix<unsigned long long int,Eigen::PlainObjectBase<Eigen::Matrix<int,-1,1,0,-1,1> > >(std::vector<unsigned long long int,std::allocator<unsigned long long int> > const &,Eigen::PlainObjectBase<Eigen::Matrix<int,-1,1,0,-1,1> > &);
-template bool igl::list_to_matrix<bool, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> > >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
-template bool igl::list_to_matrix<bool, Eigen::PlainObjectBase<Eigen::Array<bool, -1, -1, 0, -1, -1> > >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<int, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >&);
+template bool igl::list_to_matrix<double, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<float, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<unsigned int, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >(std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> >&);
+template bool igl::list_to_matrix<int, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<int, std::allocator<int> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template bool igl::list_to_matrix<unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(std::vector<unsigned long, std::allocator<unsigned long> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template bool igl::list_to_matrix<float, Eigen::Matrix<float, -1, -1, 0, -1, -1> >(std::vector<std::vector<float, std::allocator<float> >, std::allocator<std::vector<float, std::allocator<float> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<bool, Eigen::Array<bool, -1, 1, 0, -1, 1> >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
+template bool igl::list_to_matrix<bool, Eigen::Array<bool, -1, -1, 0, -1, -1> >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<bool, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<bool, std::allocator<bool> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<bool, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<bool, std::allocator<bool> >, std::allocator<std::vector<bool, std::allocator<bool> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template bool igl::list_to_matrix<double, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::__1::vector<double, std::__1::allocator<double> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/list_to_matrix.h

@@ -1,18 +1,20 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_LIST_TO_MATRIX_H
 #define IGL_LIST_TO_MATRIX_H
 #include "igl_inline.h"
 #include <vector>
+#include <Eigen/Core>
+
 namespace igl
 {
   // Convert a list (std::vector) of row vectors of the same length to a matrix
-  // Template: 
+  // Template:
   //   T  type that can be safely cast to type in Mat via '='
   //   Mat  Matrix type, must implement:
   //     .resize(m,n)
@@ -22,10 +24,10 @@ namespace igl
   // Outputs:
   //   M  an m by n matrix
   // Returns true on success, false on errors
-  template <typename T, class Mat>
+  template <typename T, typename Derived>
   IGL_INLINE bool list_to_matrix(
     const std::vector<std::vector<T > > & V,
-    Mat & M);
+    Eigen::PlainObjectBase<Derived>& M);
   // Convert a list of row vectors of `n` or less to a matrix and pad on
   // the right with `padding`:
   //
@@ -35,15 +37,15 @@ namespace igl
   //   padding  value to fill in from right for short rows
   // Outputs:
   //   M  an m by n matrix
-  template <typename T, class Mat>
+  template <typename T, typename Derived>
   IGL_INLINE bool list_to_matrix(
     const std::vector<std::vector<T > > & V,
     const int n,
     const T & padding,
-    Mat & M);
+    Eigen::PlainObjectBase<Derived>& M);
   // Vector wrapper
-  template <typename T, class Mat>
-  IGL_INLINE bool list_to_matrix(const std::vector<T > & V,Mat & M);
+  template <typename T, typename Derived>
+  IGL_INLINE bool list_to_matrix(const std::vector<T > & V,Eigen::PlainObjectBase<Derived>& M);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/lu_lagrange.cpp

@@ -1,159 +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 "lu_lagrange.h"
-
-// Cholesky LLT decomposition for symmetric positive definite
-//#include <Eigen/SparseExtra>
-// Bug in unsupported/Eigen/SparseExtra needs iostream first
-#include <iostream>
-#include <unsupported/Eigen/SparseExtra>
-#include <cassert>
-#include <cstdio>
-#include "find.h"
-#include "sparse.h"
-
-template <typename T>
-IGL_INLINE bool igl::lu_lagrange(
-  const Eigen::SparseMatrix<T> & ATA,
-  const Eigen::SparseMatrix<T> & C,
-  Eigen::SparseMatrix<T> & L,
-  Eigen::SparseMatrix<T> & U)
-{
-#if EIGEN_VERSION_AT_LEAST(3,0,92)
-#if defined(_WIN32)
-  #pragma message("lu_lagrange has not yet been implemented for your Eigen Version")
-#else
-  #warning lu_lagrange has not yet been implemented for your Eigen Version
-#endif
-
-  return false;
-#else
-  // number of unknowns
-  int n = ATA.rows();
-  // number of lagrange multipliers
-  int m = C.cols();
-
-  assert(ATA.cols() == n);
-  if(m != 0)
-  {
-    assert(C.rows() == n);
-    if(C.nonZeros() == 0)
-    {
-      // See note above about empty columns in C
-      fprintf(stderr,"Error: lu_lagrange() C has columns but no entries\n");
-      return false;
-    }
-  }
-
-  // Check that each column of C has at least one entry
-  std::vector<bool> has_entry; has_entry.resize(C.cols(),false);
-  // Iterate over outside
-  for(int k=0; k<C.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (C,k); it; ++it)
-    {
-      has_entry[it.col()] = true;
-    }
-  }
-  for(int i=0;i<(int)has_entry.size();i++)
-  {
-    if(!has_entry[i])
-    {
-      // See note above about empty columns in C
-      fprintf(stderr,"Error: lu_lagrange() C(:,%d) has no entries\n",i);
-      return false;
-    }
-  }
-
-
-
-  // Cholesky factorization of ATA
-  //// Eigen fails if you give a full view of the matrix like this:
-  //Eigen::SparseLLT<SparseMatrix<T> > ATA_LLT(ATA);
-  Eigen::SparseMatrix<T> ATA_LT = ATA.template triangularView<Eigen::Lower>();
-  Eigen::SparseLLT<Eigen::SparseMatrix<T> > ATA_LLT(ATA_LT);
-
-  Eigen::SparseMatrix<T> J = ATA_LLT.matrixL();
-
-  //if(!ATA_LLT.succeeded())
-  if(!((J*0).eval().nonZeros() == 0))
-  {
-    fprintf(stderr,"Error: lu_lagrange() failed to factor ATA\n");
-    return false;
-  }
-
-  if(m == 0)
-  {
-    // If there are no constraints (C is empty) then LU decomposition is just L
-    // and L' from cholesky decomposition
-    L = J;
-    U = J.transpose();
-  }else
-  {
-    // Construct helper matrix M
-    Eigen::SparseMatrix<T> M = C;
-    J.template triangularView<Eigen::Lower>().solveInPlace(M);
-
-    // Compute cholesky factorizaiton of M'*M
-    Eigen::SparseMatrix<T> MTM = M.transpose() * M;
-
-    Eigen::SparseLLT<Eigen::SparseMatrix<T> > MTM_LLT(MTM.template triangularView<Eigen::Lower>());
-
-    Eigen::SparseMatrix<T> K = MTM_LLT.matrixL();
-
-    //if(!MTM_LLT.succeeded())
-    if(!((K*0).eval().nonZeros() == 0))
-    {
-      fprintf(stderr,"Error: lu_lagrange() failed to factor MTM\n");
-      return false;
-    }
-
-    // assemble LU decomposition of Q
-    Eigen::Matrix<int,Eigen::Dynamic,1> MI;
-    Eigen::Matrix<int,Eigen::Dynamic,1> MJ;
-    Eigen::Matrix<T,Eigen::Dynamic,1> MV;
-    igl::find(M,MI,MJ,MV);
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> KI;
-    Eigen::Matrix<int,Eigen::Dynamic,1> KJ;
-    Eigen::Matrix<T,Eigen::Dynamic,1> KV;
-    igl::find(K,KI,KJ,KV);
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> JI;
-    Eigen::Matrix<int,Eigen::Dynamic,1> JJ;
-    Eigen::Matrix<T,Eigen::Dynamic,1> JV;
-    igl::find(J,JI,JJ,JV);
-
-    int nnz = JV.size()  + MV.size() + KV.size();
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> UI(nnz);
-    Eigen::Matrix<int,Eigen::Dynamic,1> UJ(nnz);
-    Eigen::Matrix<T,Eigen::Dynamic,1> UV(nnz);
-    UI << JJ,                        MI, (KJ.array() + n).matrix();
-    UJ << JI, (MJ.array() + n).matrix(), (KI.array() + n).matrix(); 
-    UV << JV,                        MV,                     KV*-1;
-    igl::sparse(UI,UJ,UV,U);
-
-    Eigen::Matrix<int,Eigen::Dynamic,1> LI(nnz);
-    Eigen::Matrix<int,Eigen::Dynamic,1> LJ(nnz);
-    Eigen::Matrix<T,Eigen::Dynamic,1> LV(nnz);
-    LI << JI, (MJ.array() + n).matrix(), (KI.array() + n).matrix();
-    LJ << JJ,                        MI, (KJ.array() + n).matrix(); 
-    LV << JV,                        MV,                        KV;
-    igl::sparse(LI,LJ,LV,L);
-  }
-
-  return true;
-  #endif
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template specialization
-template bool igl::lu_lagrange<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::SparseMatrix<double, 0, int>&);
-#endif

include/igl/lu_lagrange.h

@@ -1,67 +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_LU_LAGRANGE_H
-#define IGL_LU_LAGRANGE_H
-#include "igl_inline.h"
-
-#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-namespace igl
-{
-  // KNOWN BUGS: This does not seem to be correct for non-empty C
-  //
-  // LU_LAGRANGE Compute a LU decomposition for a special type of
-  // matrix Q that is symmetric but not positive-definite:
-  // Q = [A'*A C
-  //      C'   0];
-  // where A'*A, or ATA, is given as a symmetric positive definite matrix and C
-  // has full column-rank(?)
-  //
-  // [J] = lu_lagrange(ATA,C)
-  //
-  // Templates:
-  //   T  should be a eigen matrix primitive type like int or double
-  // Inputs:
-  //   ATA   n by n square, symmetric, positive-definite system matrix, usually
-  //     the quadratic coefficients corresponding to the original unknowns in a
-  //     system
-  //   C  n by m rectangular matrix corresponding the quadratic coefficients of
-  //     the original unknowns times the lagrange multipliers enforcing linear
-  //     equality constraints
-  // Outputs:
-  //   L  lower triangular matrix such that Q = L*U
-  //   U  upper triangular matrix such that Q = L*U
-  // Returns true on success, false on error
-  //
-  // Note: C should *not* have any empty columns. Typically C is the slice of
-  // the linear constraints matrix Aeq concerning the unknown variables of a
-  // quadratic optimization. Generally constraints may deal with unknowns as
-  // well as knowns. Each linear constraint corresponds to a column of Aeq. As
-  // long as each constraint concerns at least one unknown then the
-  // corresponding column in C will have at least one non zero entry. If a
-  // constraint concerns *no* unknowns, you should double check that this is a
-  // valid constraint. How can you constrain known values to each other? This
-  // is either a contradiction to the knowns' values or redundent. In either
-  // case, it's not this functions responsiblilty to handle empty constraints
-  // so you will get an error.
-  //
-  template <typename T>
-  IGL_INLINE bool lu_lagrange(
-    const Eigen::SparseMatrix<T> & ATA,
-    const Eigen::SparseMatrix<T> & C,
-    Eigen::SparseMatrix<T> & L,
-    Eigen::SparseMatrix<T> & U);
-
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "lu_lagrange.cpp"
-#endif
-
-#endif

include/igl/matlab/prepare_lhs.cpp

@@ -14,22 +14,17 @@ IGL_INLINE void igl::matlab::prepare_lhs_double(
 {
   using namespace std;
   using namespace Eigen;
-  plhs[0] = mxCreateDoubleMatrix(V.rows(),V.cols(), mxREAL);
+  const int m = V.rows();
+  const int n = V.cols();
+  plhs[0] = mxCreateDoubleMatrix(m,n, mxREAL);
   double * Vp = mxGetPr(plhs[0]);
-
-  typedef typename DerivedV::Scalar Scalar;
-  const Scalar * V_data;
-  Matrix<Scalar, DerivedV::ColsAtCompileTime, DerivedV::RowsAtCompileTime, RowMajor> VT;
-  if(DerivedV::IsRowMajor)
-  {
-    VT = V.transpose();
-    V_data = VT.data();
-  }else
+  for(int i = 0;i<m;i++)
   {
-    V_data = V.data();
+    for(int j = 0;j<n;j++)
+    {
+      Vp[i+m*j] = V(i,j);
+    }
   }
-
-  copy(V_data,V_data+V.size(),Vp);
 }
 
 template <typename DerivedV>
@@ -39,22 +34,17 @@ IGL_INLINE void igl::matlab::prepare_lhs_logical(
 {
   using namespace std;
   using namespace Eigen;
-  plhs[0] = mxCreateLogicalMatrix(V.rows(),V.cols());
+  const int m = V.rows();
+  const int n = V.cols();
+  plhs[0] = mxCreateLogicalMatrix(m,n);
   mxLogical * Vp = static_cast<mxLogical*>(mxGetData(plhs[0]));
-
-  typedef typename DerivedV::Scalar Scalar;
-  const Scalar * V_data;
-  Matrix<Scalar, DerivedV::ColsAtCompileTime, DerivedV::RowsAtCompileTime, RowMajor> VT;
-  if(DerivedV::IsRowMajor)
-  {
-    VT = V.transpose();
-    V_data = VT.data();
-  }else
+  for(int i = 0;i<m;i++)
   {
-    V_data = V.data();
+    for(int j = 0;j<n;j++)
+    {
+      Vp[i+m*j] = V(i,j);
+    }
   }
-
-  copy(V_data,V_data+V.size(),Vp);
 }
 
 template <typename DerivedV>
@@ -63,7 +53,6 @@ IGL_INLINE void igl::matlab::prepare_lhs_index(
   mxArray *plhs[])
 {
   // Treat indices as reals
-  
   const auto Vd = (V.template cast<double>().array()+1).eval();
   return prepare_lhs_double(Vd,plhs);
 }

include/igl/min_quad_with_fixed.cpp

@@ -1,6 +1,6 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 //
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 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
@@ -12,8 +12,6 @@
 #include "find.h"
 #include "sparse.h"
 #include "repmat.h"
-//#include "lu_lagrange.h"
-#include "full.h"
 #include "matlab_format.h"
 #include "EPS.h"
 #include "cat.h"
@@ -536,7 +534,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   const Eigen::PlainObjectBase<DerivedBeq> & Beq,
   Eigen::PlainObjectBase<DerivedZ> & Z)
 {
-  Eigen::PlainObjectBase<DerivedZ> sol;
+  Eigen::Matrix<typename DerivedZ::Scalar, Eigen::Dynamic, Eigen::Dynamic> sol;
   return min_quad_with_fixed_solve(data,B,Y,Beq,Z,sol);
 }
 

include/igl/mod.cpp

@@ -14,16 +14,19 @@ IGL_INLINE void igl::mod(
   Eigen::PlainObjectBase<DerivedB> & B)
 {
   B.resize(A.rows(),A.cols());
-  for(int i = 0;i<B.size();i++)
+  for(int i = 0;i<A.rows();i++)
   {
-    *(B.data()+i) = (*(A.data()+i))%base;
+    for(int j = 0;j<A.cols();j++)
+    {
+      B(i,j) = A(i,j)%base;
+    }
   }
 }
 template <typename DerivedA>
-IGL_INLINE Eigen::PlainObjectBase<DerivedA> igl::mod(
+IGL_INLINE DerivedA igl::mod(
   const Eigen::PlainObjectBase<DerivedA> & A, const int base)
 {
-  Eigen::PlainObjectBase<DerivedA> B;
+  DerivedA B;
   mod(A,base,B);
   return B;
 }

include/igl/mod.h

@@ -24,7 +24,7 @@ namespace igl
     const int base,
     Eigen::PlainObjectBase<DerivedB> & B);
   template <typename DerivedA>
-  IGL_INLINE Eigen::PlainObjectBase<DerivedA> mod(
+  IGL_INLINE DerivedA mod(
     const Eigen::PlainObjectBase<DerivedA> & A, const int base);
 }
 #ifndef IGL_STATIC_LIBRARY

include/igl/n_polyvector.cpp

@@ -39,6 +39,7 @@ namespace igl {
     Eigen::VectorXi indInteriorToFull;
     Eigen::VectorXi indFullToInterior;
 
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
 
     IGL_INLINE void computek();

include/igl/n_polyvector_general.cpp

@@ -36,6 +36,7 @@ namespace igl {
     Eigen::VectorXi indInteriorToFull;
     Eigen::VectorXi indFullToInterior;
 
+#warning "Constructing Eigen::PlainObjectBase directly is deprecated"
     Eigen::PlainObjectBase<DerivedV> B1, B2, FN;
 
     IGL_INLINE void computek();

include/igl/per_corner_normals.cpp

@@ -11,23 +11,27 @@
 #include "per_face_normals.h"
 #include "PI.h"
 
-template <typename DerivedV, typename DerivedF>
+template <typename DerivedV, typename DerivedF, typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   const double corner_threshold,
-  Eigen::PlainObjectBase<DerivedV> & CN)
+  Eigen::PlainObjectBase<DerivedCN> & CN)
 {
   using namespace Eigen;
   using namespace std;
-  Eigen::PlainObjectBase<DerivedV> FN;
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> FN;
   per_face_normals(V,F,FN);
   vector<vector<int> > VF,VFi;
   vertex_triangle_adjacency(V,F,VF,VFi);
   return per_corner_normals(V,F,FN,VF,corner_threshold,CN);
 }
 
-template <typename DerivedV, typename DerivedF, typename DerivedFN, typename DerivedCN>
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedFN, 
+  typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
@@ -42,14 +46,19 @@ IGL_INLINE void igl::per_corner_normals(
   return per_corner_normals(V,F,FN,VF,corner_threshold,CN);
 }
 
-template <typename DerivedV, typename DerivedF, typename IndexType>
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedFN, 
+  typename IndexType,
+  typename DerivedCN>
 IGL_INLINE void igl::per_corner_normals(
   const Eigen::PlainObjectBase<DerivedV>& /*V*/,
   const Eigen::PlainObjectBase<DerivedF>& F,
-  const Eigen::PlainObjectBase<DerivedV>& FN,
+  const Eigen::PlainObjectBase<DerivedFN>& FN,
   const std::vector<std::vector<IndexType> >& VF,
   const double corner_threshold,
-  Eigen::PlainObjectBase<DerivedV> & CN)
+  Eigen::PlainObjectBase<DerivedCN> & CN)
 {
   using namespace Eigen;
   using namespace std;
@@ -96,9 +105,7 @@ IGL_INLINE void igl::per_corner_normals(
 // Explicit template specialization
 // generated by autoexplicit.sh
 template void igl::per_corner_normals<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&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::per_corner_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, unsigned int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, std::vector<std::vector<unsigned int, std::allocator<unsigned int> >, std::allocator<std::vector<unsigned int, std::allocator<unsigned int> > > > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
 template void igl::per_corner_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::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&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::per_corner_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_corner_normals<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&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_corner_normals<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<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&, double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/per_corner_normals.h

@@ -21,12 +21,12 @@ namespace igl
   // Output:
   //   CN  #F*3 by 3 eigen Matrix of mesh vertex 3D normals, where the normal
   //     for corner F(i,j) is at CN(i*3+j,:) 
-  template <typename DerivedV, typename DerivedF>
+  template <typename DerivedV, typename DerivedF, typename DerivedCN>
   IGL_INLINE void per_corner_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
     const double corner_threshold,
-    Eigen::PlainObjectBase<DerivedV> & CN);
+    Eigen::PlainObjectBase<DerivedCN> & CN);
   // Other Inputs:
   //   FN  #F by 3 eigen Matrix of face normals
   template <
@@ -42,14 +42,19 @@ namespace igl
     Eigen::PlainObjectBase<DerivedCN> & CN);
   // Other Inputs:
   //   VF  map from vertices to list of incident faces
-  template <typename DerivedV, typename DerivedF, typename IndexType>
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedFN, 
+    typename IndexType,
+    typename DerivedCN>
   IGL_INLINE void per_corner_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
-    const Eigen::PlainObjectBase<DerivedV>& FN,
+    const Eigen::PlainObjectBase<DerivedFN>& FN,
     const std::vector<std::vector<IndexType> >& VF,
     const double corner_threshold,
-    Eigen::PlainObjectBase<DerivedV> & CN);
+    Eigen::PlainObjectBase<DerivedCN> & CN);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/per_edge_normals.cpp

@@ -90,7 +90,7 @@ IGL_INLINE void igl::per_edge_normals(
   Eigen::PlainObjectBase<DerivedE> & E,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
 {
-  Eigen::PlainObjectBase<DerivedN> FN;
+  Eigen::Matrix<typename DerivedN::Scalar,Eigen::Dynamic,3> FN;
   per_face_normals(V,F,FN);
   return per_edge_normals(V,F,weighting,FN,N,E,EMAP);
 }

include/igl/per_face_normals.cpp

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

include/igl/per_vertex_normals.cpp

@@ -19,7 +19,7 @@ IGL_INLINE void igl::per_vertex_normals(
   const igl::PerVertexNormalsWeightingType weighting,
   Eigen::PlainObjectBase<DerivedV> & N)
 {
-  Eigen::PlainObjectBase<DerivedV> PFN;
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> PFN;
   igl::per_face_normals(V,F,PFN);
   return per_vertex_normals(V,F,weighting,PFN,N);
 }

include/igl/piecewise_constant_winding_number.cpp

@@ -0,0 +1,70 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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 "piecewise_constant_winding_number.h"
+#include "unique_edge_map.h"
+
+template <
+  typename DerivedF,
+  typename DeriveduE,
+  typename uE2EType>
+IGL_INLINE bool igl::piecewise_constant_winding_number(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DeriveduE>& uE,
+  const std::vector<std::vector<uE2EType> >& uE2E)
+{
+  const size_t num_faces = F.rows();
+  const size_t num_edges = uE.rows();
+  const auto edge_index_to_face_index = [&](size_t ei)
+  {
+    return ei % num_faces;
+  };
+  const auto is_consistent = [&](size_t fid, size_t s, size_t d)
+  {
+    if ((size_t)F(fid, 0) == s && (size_t)F(fid, 1) == d) return true;
+    if ((size_t)F(fid, 1) == s && (size_t)F(fid, 2) == d) return true;
+    if ((size_t)F(fid, 2) == s && (size_t)F(fid, 0) == d) return true;
+
+    if ((size_t)F(fid, 0) == d && (size_t)F(fid, 1) == s) return false;
+    if ((size_t)F(fid, 1) == d && (size_t)F(fid, 2) == s) return false;
+    if ((size_t)F(fid, 2) == d && (size_t)F(fid, 0) == s) return false;
+    throw "Invalid face!!";
+  };
+  for (size_t i=0; i<num_edges; i++)
+  {
+    const size_t s = uE(i,0);
+    const size_t d = uE(i,1);
+    int count=0;
+    for (const auto& ei : uE2E[i])
+    {
+      const size_t fid = edge_index_to_face_index(ei);
+      if (is_consistent(fid, s, d))
+      {
+        count++;
+      }
+      else
+      {
+        count--;
+      }
+    }
+    if (count != 0)
+    {
+      return false;
+    }
+  }
+  return true;
+}
+template <typename DerivedF>
+IGL_INLINE bool igl::piecewise_constant_winding_number(
+  const Eigen::PlainObjectBase<DerivedF>& F)
+{
+  Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,2> E, uE;
+  Eigen::Matrix<typename DerivedF::Scalar,Eigen::Dynamic,1> EMAP;
+  std::vector<std::vector<size_t> > uE2E;
+  unique_edge_map(F, E, uE, EMAP, uE2E);
+  return piecewise_constant_winding_number(F,uE,uE2E);
+}

include/igl/piecewise_constant_winding_number.h

@@ -0,0 +1,49 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_PIECEWISE_CONSTANT_WINDING_NUMBER_H
+#define IGL_PIECEWISE_CONSTANT_WINDING_NUMBER_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // PIECEWISE_CONSTANT_WINDING_NUMBER Determine if a given mesh induces a
+  // piecewise constant winding number field: Is this mesh valid input to solid
+  // set operations.  **Assumes** that `(V,F)` contains no self-intersections
+  // (including degeneracies and co-incidences).  If there are co-planar and
+  // co-incident vertex placements, a mesh could _fail_ this combinatorial test
+  // but still induce a piecewise-constant winding number _geometrically_. For
+  // example, consider a hemisphere with boundary and then pinch the boundary
+  // "shut" along a line segment. The **_bullet-proof_** check is to first
+  // resolve all self-intersections in `(V,F) -> (SV,SF)` (i.e. what the
+  // `igl::copyleft::cgal::piecewise_constant_winding_number` overload does).
+  //
+  // Inputs:
+  //   F  #F by 3 list of triangle indices into some (abstract) list of
+  //     vertices V
+  //   uE  #uE by 2 list of unique edges indices into V
+  //   uE2E  #uE list of lists of indices into directed edges (#F * 3)
+  // Returns true if the mesh _combinatorially_ induces a piecewise constant
+  // winding number field.
+  //
+  template <
+    typename DerivedF,
+    typename DeriveduE,
+    typename uE2EType>
+  IGL_INLINE bool piecewise_constant_winding_number(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DeriveduE>& uE,
+    const std::vector<std::vector<uE2EType> >& uE2E);
+  template <typename DerivedF>
+  IGL_INLINE bool piecewise_constant_winding_number(
+    const Eigen::PlainObjectBase<DerivedF>& F);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "piecewise_constant_winding_number.cpp"
+#endif
+#endif

include/igl/ply.h.REMOVED.git-id

@@ -1 +1 @@
-31ac787262543763805c0a1c35a57998ebb1628c
+a11a359439a307b36f2320e726048ce3dd344fda

include/igl/polar_dec.cpp

@@ -84,8 +84,8 @@ IGL_INLINE void igl::polar_dec(
   Eigen::PlainObjectBase<DerivedR> & R,
   Eigen::PlainObjectBase<DerivedT> & T)
 {
-  Eigen::PlainObjectBase<DerivedA> U;
-  Eigen::PlainObjectBase<DerivedA> V;
+  DerivedA U;
+  DerivedA V;
   Eigen::Matrix<typename DerivedA::Scalar,DerivedA::RowsAtCompileTime,1> S;
   return igl::polar_dec(A,R,T,U,S,V);
 }

include/igl/polar_svd.cpp

@@ -20,8 +20,8 @@ IGL_INLINE void igl::polar_svd(
   Eigen::PlainObjectBase<DerivedR> & R,
   Eigen::PlainObjectBase<DerivedT> & T)
 {
-  Eigen::PlainObjectBase<DerivedA> U;
-  Eigen::PlainObjectBase<DerivedA> V;
+  DerivedA U;
+  DerivedA V;
   Eigen::Matrix<typename DerivedA::Scalar,DerivedA::RowsAtCompileTime,1> S;
   return igl::polar_svd(A,R,T,U,S,V);
 }

include/igl/polar_svd3x3.cpp

@@ -85,7 +85,7 @@ IGL_INLINE void igl::polar_svd3x3_avx(const Eigen::Matrix<T, 3*8, 3>& A, Eigen::
     Eigen::Matrix3f Rpart_SSE = R.block(3*k, 0, 3, 3);
     Eigen::Matrix3f diff = Rpart - Rpart_SSE;
     float diffNorm = diff.norm();
-    if (abs(diffNorm) > 0.001)
+    if (std::abs(diffNorm) > 0.001)
     {
       printf("Huh: diffNorm = %15f (k = %i)\n", diffNorm, k);
     }

include/igl/polyvector_field_comb_from_matchings_and_cuts.cpp

@@ -96,7 +96,7 @@ IGL_INLINE void igl::polyvector_field_comb_from_matchings_and_cuts(
   Eigen::PlainObjectBase<DerivedS> &sol3D_combed)
   {
     Eigen::MatrixXi TT, TTi;
-    igl::triangle_triangle_adjacency(V,F,TT,TTi);
+    igl::triangle_triangle_adjacency(F,TT,TTi);
 
     Eigen::MatrixXi E, E2F, F2E;
     igl::edge_topology(V,F,E,F2E,E2F);

include/igl/polyvector_field_cut_mesh_with_singularities.cpp

@@ -91,7 +91,7 @@ IGL_INLINE void igl::polyvector_field_cut_mesh_with_singularities(
   igl::adjacency_list(F, VV);
 
   Eigen::MatrixXi TT, TTi;
-  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  igl::triangle_triangle_adjacency(F,TT,TTi);
 
   igl::polyvector_field_cut_mesh_with_singularities(V, F, VF, VV, TT, TTi, singularities, cuts);
 

include/igl/polyvector_field_matchings.cpp

@@ -168,7 +168,8 @@ IGL_INLINE typename DerivedC::Scalar igl::polyvector_field_matchings(
   Eigen::MatrixXi E, E2F, F2E;
   igl::edge_topology(V,F,E,F2E,E2F);
 
-  Eigen::PlainObjectBase<DerivedV> FN;
+#warning "Poor templating of igl::polyvector_field_matchings forces FN to be DerivedV (this could cause issues if DerivedV has fixed number of rows)"
+  DerivedV FN;
   igl::per_face_normals(V,F,FN);
 
   return igl::polyvector_field_matchings(sol3D, V, F, E, FN, E2F, match_with_curl, match_ab, match_ba, curl);

include/igl/polyvector_field_singularities_from_matchings.cpp

@@ -80,7 +80,7 @@ IGL_INLINE void igl::polyvector_field_singularities_from_matchings(
   igl::vertex_triangle_adjacency(V,F,VF,VFi);
 
   Eigen::MatrixXi TT, TTi;
-  igl::triangle_triangle_adjacency(V,F,TT,TTi);
+  igl::triangle_triangle_adjacency(F,TT,TTi);
 
   Eigen::MatrixXi E, E2F, F2E;
   igl::edge_topology(V,F,E,F2E,E2F);

include/igl/project_to_line.cpp

@@ -1,18 +1,18 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "project_to_line.h"
 #include <cassert>
 #include <Eigen/Core>
 
 template <
-  typename MatP, 
-  typename MatL, 
-  typename Matt, 
+  typename MatP,
+  typename MatL,
+  typename Matt,
   typename MatsqrD>
 IGL_INLINE void igl::project_to_line(
   const MatP & P,
@@ -38,11 +38,11 @@ IGL_INLINE void igl::project_to_line(
   // resize output
   t.resize(np,1);
   sqrD.resize(np,1);
-  // loop over points 
+  // loop over points
 #pragma omp parallel for if (np>10000)
   for(int i = 0;i<np;i++)
   {
-    MatL Pi = P.row(i);
+    MatP Pi = P.row(i);
     // vector from point i to source
     MatL SmPi = S-Pi;
     t(i) = -(DmS.array()*SmPi.array()).sum() / v_sqrlen;

include/igl/project_to_line_segment.cpp

@@ -1,19 +1,19 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "project_to_line_segment.h"
 #include "project_to_line.h"
 #include <Eigen/Core>
 
 template <
-  typename DerivedP, 
-  typename DerivedS, 
-  typename DerivedD, 
-  typename Derivedt, 
+  typename DerivedP,
+  typename DerivedS,
+  typename DerivedD,
+  typename Derivedt,
   typename DerivedsqrD>
 IGL_INLINE void igl::project_to_line_segment(
   const Eigen::PlainObjectBase<DerivedP> & P,
@@ -28,7 +28,7 @@ IGL_INLINE void igl::project_to_line_segment(
 #pragma omp parallel for if (np>10000)
   for(int p = 0;p<np;p++)
   {
-    const Eigen::PlainObjectBase<DerivedS> Pp = P.row(p);
+    const Eigen::PlainObjectBase<DerivedP> Pp = P.row(p);
     if(t(p)<0)
     {
       sqrD(p) = (Pp-S).squaredNorm();

include/igl/randperm.cpp

@@ -20,16 +20,16 @@ IGL_INLINE void igl::randperm(
   std::random_shuffle(I.data(),I.data()+n);
 }
 
-template <typename DerivedI>
-IGL_INLINE Eigen::PlainObjectBase<DerivedI> igl::randperm( const int n)
-{
-  Eigen::PlainObjectBase<DerivedI> I;
-  randperm(n,I);
-  return I;
-}
+//template <typename DerivedI>
+//IGL_INLINE Eigen::PlainObjectBase<DerivedI> igl::randperm( const int n)
+//{
+//  Eigen::PlainObjectBase<DerivedI> I;
+//  randperm(n,I);
+//  return I;
+//}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int);
+//template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > igl::randperm<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int);
 #endif

include/igl/randperm.h

@@ -21,8 +21,8 @@ namespace igl
   IGL_INLINE void randperm(
     const int n,
     Eigen::PlainObjectBase<DerivedI> & I);
-  template <typename DerivedI>
-  IGL_INLINE Eigen::PlainObjectBase<DerivedI> randperm( const int n);
+  //template <typename DerivedI>
+  //IGL_INLINE DerivedI randperm( const int n);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "randperm.cpp"

include/igl/ray_box_intersect.cpp

@@ -0,0 +1,138 @@
+#include "ray_box_intersect.h"
+#include <vector>
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename Scalar>
+IGL_INLINE bool igl::ray_box_intersect(
+  const Eigen::PlainObjectBase<Derivedsource> & origin,
+  const Eigen::PlainObjectBase<Deriveddir> & dir,
+  const Eigen::AlignedBox<Scalar,3> & box,
+  const Scalar & t0,
+  const Scalar & t1,
+  Scalar & tmin,
+  Scalar & tmax)
+{
+#ifdef false
+  // https://github.com/RMonica/basic_next_best_view/blob/master/src/RayTracer.cpp
+  const auto & intersectRayBox = [](
+    const Eigen::Vector3f& rayo,
+    const Eigen::Vector3f& rayd,
+    const Eigen::Vector3f& bmin,
+    const Eigen::Vector3f& bmax, 
+    float & tnear,
+    float & tfar
+    )->bool
+  {
+    Eigen::Vector3f bnear;
+    Eigen::Vector3f bfar;
+    // Checks for intersection testing on each direction coordinate
+    // Computes 
+    float t1, t2;
+    tnear = -1e+6f, tfar = 1e+6f; //, tCube;
+    bool intersectFlag = true;
+    for (int i = 0; i < 3; ++i) {
+  //    std::cout << "coordinate " << i << ": bmin " << bmin(i) << ", bmax " << bmax(i) << std::endl; 
+      assert(bmin(i) <= bmax(i));
+      if (::fabs(rayd(i)) < 1e-6) {   // Ray parallel to axis i-th
+        if (rayo(i) < bmin(i) || rayo(i) > bmax(i)) {
+          intersectFlag = false;
+        }
+      }
+      else {
+        // Finds the nearest and the farthest vertices of the box from the ray origin
+        if (::fabs(bmin(i) - rayo(i)) < ::fabs(bmax(i) - rayo(i))) {
+          bnear(i) = bmin(i);
+          bfar(i) = bmax(i);
+        }
+        else {
+          bnear(i) = bmax(i);
+          bfar(i) = bmin(i);
+        }
+  //      std::cout << "  bnear " << bnear(i) << ", bfar " << bfar(i) << std::endl;
+        // Finds the distance parameters t1 and t2 of the two ray-box intersections:
+        // t1 must be the closest to the ray origin rayo. 
+        t1 = (bnear(i) - rayo(i)) / rayd(i);
+        t2 = (bfar(i) - rayo(i)) / rayd(i);
+        if (t1 > t2) {
+          std::swap(t1,t2);
+        } 
+        // The two intersection values are used to saturate tnear and tfar
+        if (t1 > tnear) {
+          tnear = t1;
+        }
+        if (t2 < tfar) {
+          tfar = t2;
+        }
+  //      std::cout << "  t1 " << t1 << ", t2 " << t2 << ", tnear " << tnear << ", tfar " << tfar 
+  //        << "  tnear > tfar? " << (tnear > tfar) << ", tfar < 0? " << (tfar < 0) << std::endl;
+        if(tnear > tfar) {
+          intersectFlag = false;
+        }
+        if(tfar < 0) {
+        intersectFlag = false;
+        }
+      }
+    }
+    // Checks whether intersection occurs or not
+    return intersectFlag;
+  };
+  float tmin_f, tmax_f;
+  bool ret = intersectRayBox(
+      origin.   template cast<float>(),
+      dir.      template cast<float>(),
+      box.min().template cast<float>(),
+      box.max().template cast<float>(),
+      tmin_f,
+      tmax_f);
+  tmin = tmin_f;
+  tmax = tmax_f;
+  return ret;
+#else
+  using namespace Eigen;
+  // This should be precomputed and provided as input
+  typedef Matrix<Scalar,1,3>  RowVector3S;
+  const RowVector3S inv_dir( 1./dir(0),1./dir(1),1./dir(2));
+  const std::vector<bool> sign = { inv_dir(0)<0, inv_dir(1)<0, inv_dir(2)<0};
+  // http://people.csail.mit.edu/amy/papers/box-jgt.pdf
+  // "An Efficient and Robust Ray–Box Intersection Algorithm"
+  Scalar tymin, tymax, tzmin, tzmax;
+  std::vector<RowVector3S> bounds = {box.min(),box.max()};
+  tmin = ( bounds[sign[0]](0)   - origin(0)) * inv_dir(0);
+  tmax = ( bounds[1-sign[0]](0) - origin(0)) * inv_dir(0);
+  tymin = (bounds[sign[1]](1)   - origin(1)) * inv_dir(1);
+  tymax = (bounds[1-sign[1]](1) - origin(1)) * inv_dir(1);
+  if ( (tmin > tymax) || (tymin > tmax) )
+  {
+    return false;
+  }
+  if (tymin > tmin)
+  {
+    tmin = tymin;
+  }
+  if (tymax < tmax)
+  {
+    tmax = tymax;
+  }
+  tzmin = (bounds[sign[2]](2) - origin(2))   * inv_dir(2);
+  tzmax = (bounds[1-sign[2]](2) - origin(2)) * inv_dir(2);
+  if ( (tmin > tzmax) || (tzmin > tmax) )
+  {
+    return false;
+  }
+  if (tzmin > tmin)
+  {
+    tmin = tzmin;
+  }
+  if (tzmax < tmax)
+  {
+    tmax = tzmax;
+  }
+  if(!( (tmin < t1) && (tmax > t0) ))
+  {
+    return false;
+  }
+  return true;
+#endif
+}

include/igl/ray_box_intersect.h

@@ -0,0 +1,44 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_RAY_BOX_INTERSECT_H
+#define IGL_RAY_BOX_INTERSECT_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+namespace igl
+{
+  // Determine whether a ray origin+t*dir and box intersect within the ray's parameterized
+  // range (t0,t1)
+  //
+  // Inputs:
+  //   source  3-vector origin of ray
+  //   dir  3-vector direction of ray
+  //   box  axis aligned box
+  //   t0  hit only if hit.t less than t0
+  //   t1  hit only if hit.t greater than t1
+  // Outputs:
+  //   tmin  minimum of interval of overlap within [t0,t1]
+  //   tmax  maximum of interval of overlap within [t0,t1]
+  // Returns true if hit
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename Scalar>
+  IGL_INLINE bool ray_box_intersect(
+    const Eigen::PlainObjectBase<Derivedsource> & source,
+    const Eigen::PlainObjectBase<Deriveddir> & dir,
+    const Eigen::AlignedBox<Scalar,3> & box,
+    const Scalar & t0,
+    const Scalar & t1,
+    Scalar & tmin,
+    Scalar & tmax);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "ray_box_intersect.cpp"
+#endif
+#endif

include/igl/ray_mesh_intersect.cpp

@@ -0,0 +1,73 @@
+#include "ray_mesh_intersect.h"
+
+extern "C"
+{
+#include "raytri.c"
+}
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename DerivedV, 
+  typename DerivedF> 
+IGL_INLINE bool igl::ray_mesh_intersect(
+  const Eigen::PlainObjectBase<Derivedsource> & s,
+  const Eigen::PlainObjectBase<Deriveddir> & dir,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  std::vector<igl::Hit> & hits)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  // Should be but can't be const 
+  Vector3d s_d = s.template cast<double>();
+  Vector3d dir_d = dir.template cast<double>();
+  hits.clear();
+  // loop over all triangles
+  for(int f = 0;f<F.rows();f++)
+  {
+    // Should be but can't be const 
+    RowVector3d v0 = V.row(F(f,0)).template cast<double>();
+    RowVector3d v1 = V.row(F(f,1)).template cast<double>();
+    RowVector3d v2 = V.row(F(f,2)).template cast<double>();
+    // shoot ray, record hit
+    double t,u,v;
+    if(intersect_triangle1(
+      s_d.data(), dir_d.data(), v0.data(), v1.data(), v2.data(), &t, &u, &v) &&
+      t>0)
+    {
+      hits.push_back({(int)f,(int)-1,(float)u,(float)v,(float)t});
+    }
+  }
+  // Sort hits based on distance
+  std::sort(
+    hits.begin(),
+    hits.end(),
+    [](const Hit & a, const Hit & b)->bool{ return a.t < b.t;});
+  return hits.size() > 0;
+}
+
+template <
+  typename Derivedsource,
+  typename Deriveddir,
+  typename DerivedV, 
+  typename DerivedF> 
+IGL_INLINE bool igl::ray_mesh_intersect(
+  const Eigen::PlainObjectBase<Derivedsource> & source,
+  const Eigen::PlainObjectBase<Deriveddir> & dir,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  igl::Hit & hit)
+{
+  std::vector<igl::Hit> hits;
+  ray_mesh_intersect(source,dir,V,F,hits);
+  if(hits.size() > 0)
+  {
+    hit = hits.front();
+    return true;
+  }else
+  {
+    return false;
+  }
+}

include/igl/ray_mesh_intersect.h

@@ -0,0 +1,49 @@
+#ifndef IGL_RAY_MESH_INTERSECT_H
+#define IGL_RAY_MESH_INTERSECT_H
+#include "igl_inline.h"
+#include "Hit.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Shoot a ray against a mesh (V,F) and collect all hits.
+  //
+  // Inputs:
+  //   source  3-vector origin of ray
+  //   dir  3-vector direction of ray
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of mesh face indices into V
+  // Outputs:
+  //    hits  **sorted** list of hits
+  // Returns true if there were any hits (hits.size() > 0)
+  //
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename DerivedV, 
+    typename DerivedF> 
+  IGL_INLINE bool ray_mesh_intersect(
+    const Eigen::PlainObjectBase<Derivedsource> & source,
+    const Eigen::PlainObjectBase<Deriveddir> & dir,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    std::vector<igl::Hit> & hits);
+  // Outputs:
+  //   hit  first hit, set only if it exists
+  // Returns true if there was a hit
+  template <
+    typename Derivedsource,
+    typename Deriveddir,
+    typename DerivedV, 
+    typename DerivedF> 
+  IGL_INLINE bool ray_mesh_intersect(
+    const Eigen::PlainObjectBase<Derivedsource> & source,
+    const Eigen::PlainObjectBase<Deriveddir> & dir,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    igl::Hit & hit);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "ray_mesh_intersect.cpp"
+#endif
+#endif

include/igl/raytri.c

@@ -0,0 +1,264 @@
+/* Ray-Triangle Intersection Test Routines          */
+/* Different optimizations of my and Ben Trumbore's */
+/* code from journals of graphics tools (JGT)       */
+/* http://www.acm.org/jgt/                          */
+/* by Tomas Moller, May 2000                        */
+
+
+// Alec: I've added an include guard, made all functions inline and added
+// IGL_RAY_TRI_ to #define macros
+#ifndef IGL_RAY_TRI_C
+#define IGL_RAY_TRI_C
+
+#include <math.h>
+
+#define IGL_RAY_TRI_EPSILON 0.000001
+#define IGL_RAY_TRI_CROSS(dest,v1,v2) \
+          dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
+          dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
+          dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
+#define IGL_RAY_TRI_DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
+#define IGL_RAY_TRI_SUB(dest,v1,v2) \
+          dest[0]=v1[0]-v2[0]; \
+          dest[1]=v1[1]-v2[1]; \
+          dest[2]=v1[2]-v2[2]; 
+
+/* the original jgt code */
+inline int intersect_triangle(double orig[3], double dir[3],
+		       double vert0[3], double vert1[3], double vert2[3],
+		       double *t, double *u, double *v)
+{
+   double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
+   double det,inv_det;
+
+   /* find vectors for two edges sharing vert0 */
+   IGL_RAY_TRI_SUB(edge1, vert1, vert0);
+   IGL_RAY_TRI_SUB(edge2, vert2, vert0);
+
+   /* begin calculating determinant - also used to calculate U parameter */
+   IGL_RAY_TRI_CROSS(pvec, dir, edge2);
+
+   /* if determinant is near zero, ray lies in plane of triangle */
+   det = IGL_RAY_TRI_DOT(edge1, pvec);
+
+   if (det > -IGL_RAY_TRI_EPSILON && det < IGL_RAY_TRI_EPSILON)
+     return 0;
+   inv_det = 1.0 / det;
+
+   /* calculate distance from vert0 to ray origin */
+   IGL_RAY_TRI_SUB(tvec, orig, vert0);
+
+   /* calculate U parameter and test bounds */
+   *u = IGL_RAY_TRI_DOT(tvec, pvec) * inv_det;
+   if (*u < 0.0 || *u > 1.0)
+     return 0;
+
+   /* prepare to test V parameter */
+   IGL_RAY_TRI_CROSS(qvec, tvec, edge1);
+
+   /* calculate V parameter and test bounds */
+   *v = IGL_RAY_TRI_DOT(dir, qvec) * inv_det;
+   if (*v < 0.0 || *u + *v > 1.0)
+     return 0;
+
+   /* calculate t, ray intersects triangle */
+   *t = IGL_RAY_TRI_DOT(edge2, qvec) * inv_det;
+
+   return 1;
+}
+
+
+/* code rewritten to do tests on the sign of the determinant */
+/* the division is at the end in the code                    */
+inline int intersect_triangle1(double orig[3], double dir[3],
+			double vert0[3], double vert1[3], double vert2[3],
+			double *t, double *u, double *v)
+{
+   double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
+   double det,inv_det;
+
+   /* find vectors for two edges sharing vert0 */
+   IGL_RAY_TRI_SUB(edge1, vert1, vert0);
+   IGL_RAY_TRI_SUB(edge2, vert2, vert0);
+
+   /* begin calculating determinant - also used to calculate U parameter */
+   IGL_RAY_TRI_CROSS(pvec, dir, edge2);
+
+   /* if determinant is near zero, ray lies in plane of triangle */
+   det = IGL_RAY_TRI_DOT(edge1, pvec);
+
+   if (det > IGL_RAY_TRI_EPSILON)
+   {
+      /* calculate distance from vert0 to ray origin */
+      IGL_RAY_TRI_SUB(tvec, orig, vert0);
+      
+      /* calculate U parameter and test bounds */
+      *u = IGL_RAY_TRI_DOT(tvec, pvec);
+      if (*u < 0.0 || *u > det)
+	 return 0;
+      
+      /* prepare to test V parameter */
+      IGL_RAY_TRI_CROSS(qvec, tvec, edge1);
+      
+      /* calculate V parameter and test bounds */
+      *v = IGL_RAY_TRI_DOT(dir, qvec);
+      if (*v < 0.0 || *u + *v > det)
+	 return 0;
+      
+   }
+   else if(det < -IGL_RAY_TRI_EPSILON)
+   {
+      /* calculate distance from vert0 to ray origin */
+      IGL_RAY_TRI_SUB(tvec, orig, vert0);
+      
+      /* calculate U parameter and test bounds */
+      *u = IGL_RAY_TRI_DOT(tvec, pvec);
+/*      printf("*u=%f\n",(float)*u); */
+/*      printf("det=%f\n",det); */
+      if (*u > 0.0 || *u < det)
+	 return 0;
+      
+      /* prepare to test V parameter */
+      IGL_RAY_TRI_CROSS(qvec, tvec, edge1);
+      
+      /* calculate V parameter and test bounds */
+      *v = IGL_RAY_TRI_DOT(dir, qvec) ;
+      if (*v > 0.0 || *u + *v < det)
+	 return 0;
+   }
+   else return 0;  /* ray is parallell to the plane of the triangle */
+
+
+   inv_det = 1.0 / det;
+
+   /* calculate t, ray intersects triangle */
+   *t = IGL_RAY_TRI_DOT(edge2, qvec) * inv_det;
+   (*u) *= inv_det;
+   (*v) *= inv_det;
+
+   return 1;
+}
+
+/* code rewritten to do tests on the sign of the determinant */
+/* the division is before the test of the sign of the det    */
+inline int intersect_triangle2(double orig[3], double dir[3],
+			double vert0[3], double vert1[3], double vert2[3],
+			double *t, double *u, double *v)
+{
+   double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
+   double det,inv_det;
+
+   /* find vectors for two edges sharing vert0 */
+   IGL_RAY_TRI_SUB(edge1, vert1, vert0);
+   IGL_RAY_TRI_SUB(edge2, vert2, vert0);
+
+   /* begin calculating determinant - also used to calculate U parameter */
+   IGL_RAY_TRI_CROSS(pvec, dir, edge2);
+
+   /* if determinant is near zero, ray lies in plane of triangle */
+   det = IGL_RAY_TRI_DOT(edge1, pvec);
+
+   /* calculate distance from vert0 to ray origin */
+   IGL_RAY_TRI_SUB(tvec, orig, vert0);
+   inv_det = 1.0 / det;
+   
+   if (det > IGL_RAY_TRI_EPSILON)
+   {
+      /* calculate U parameter and test bounds */
+      *u = IGL_RAY_TRI_DOT(tvec, pvec);
+      if (*u < 0.0 || *u > det)
+	 return 0;
+      
+      /* prepare to test V parameter */
+      IGL_RAY_TRI_CROSS(qvec, tvec, edge1);
+      
+      /* calculate V parameter and test bounds */
+      *v = IGL_RAY_TRI_DOT(dir, qvec);
+      if (*v < 0.0 || *u + *v > det)
+	 return 0;
+      
+   }
+   else if(det < -IGL_RAY_TRI_EPSILON)
+   {
+      /* calculate U parameter and test bounds */
+      *u = IGL_RAY_TRI_DOT(tvec, pvec);
+      if (*u > 0.0 || *u < det)
+	 return 0;
+      
+      /* prepare to test V parameter */
+      IGL_RAY_TRI_CROSS(qvec, tvec, edge1);
+      
+      /* calculate V parameter and test bounds */
+      *v = IGL_RAY_TRI_DOT(dir, qvec) ;
+      if (*v > 0.0 || *u + *v < det)
+	 return 0;
+   }
+   else return 0;  /* ray is parallell to the plane of the triangle */
+
+   /* calculate t, ray intersects triangle */
+   *t = IGL_RAY_TRI_DOT(edge2, qvec) * inv_det;
+   (*u) *= inv_det;
+   (*v) *= inv_det;
+
+   return 1;
+}
+
+/* code rewritten to do tests on the sign of the determinant */
+/* the division is before the test of the sign of the det    */
+/* and one IGL_RAY_TRI_CROSS has been moved out from the if-else if-else */
+inline int intersect_triangle3(double orig[3], double dir[3],
+			double vert0[3], double vert1[3], double vert2[3],
+			double *t, double *u, double *v)
+{
+   double edge1[3], edge2[3], tvec[3], pvec[3], qvec[3];
+   double det,inv_det;
+
+   /* find vectors for two edges sharing vert0 */
+   IGL_RAY_TRI_SUB(edge1, vert1, vert0);
+   IGL_RAY_TRI_SUB(edge2, vert2, vert0);
+
+   /* begin calculating determinant - also used to calculate U parameter */
+   IGL_RAY_TRI_CROSS(pvec, dir, edge2);
+
+   /* if determinant is near zero, ray lies in plane of triangle */
+   det = IGL_RAY_TRI_DOT(edge1, pvec);
+
+   /* calculate distance from vert0 to ray origin */
+   IGL_RAY_TRI_SUB(tvec, orig, vert0);
+   inv_det = 1.0 / det;
+   
+   IGL_RAY_TRI_CROSS(qvec, tvec, edge1);
+      
+   if (det > IGL_RAY_TRI_EPSILON)
+   {
+      *u = IGL_RAY_TRI_DOT(tvec, pvec);
+      if (*u < 0.0 || *u > det)
+	 return 0;
+            
+      /* calculate V parameter and test bounds */
+      *v = IGL_RAY_TRI_DOT(dir, qvec);
+      if (*v < 0.0 || *u + *v > det)
+	 return 0;
+      
+   }
+   else if(det < -IGL_RAY_TRI_EPSILON)
+   {
+      /* calculate U parameter and test bounds */
+      *u = IGL_RAY_TRI_DOT(tvec, pvec);
+      if (*u > 0.0 || *u < det)
+	 return 0;
+      
+      /* calculate V parameter and test bounds */
+      *v = IGL_RAY_TRI_DOT(dir, qvec) ;
+      if (*v > 0.0 || *u + *v < det)
+	 return 0;
+   }
+   else return 0;  /* ray is parallell to the plane of the triangle */
+
+   *t = IGL_RAY_TRI_DOT(edge2, qvec) * inv_det;
+   (*u) *= inv_det;
+   (*v) *= inv_det;
+
+   return 1;
+}
+#endif

include/igl/remove_duplicate_vertices.cpp

@@ -26,7 +26,7 @@ IGL_INLINE void igl::remove_duplicate_vertices(
 {
   if(epsilon > 0)
   {
-    Eigen::PlainObjectBase<DerivedV> rV,rSV;
+    DerivedV rV,rSV;
     round((V/(10.0*epsilon)).eval(),rV);
     unique_rows(rV,rSV,SVI,SVJ);
     slice(V,SVI,colon<int>(0,V.cols()-1),SV);

include/igl/signed_distance.cpp

@@ -154,7 +154,7 @@ IGL_INLINE void igl::signed_distance(
           signed_distance_pseudonormal(tree3,V,F,FN,VN,EN,EMAP,q3,s,sqrd,i,c3,n3):
           signed_distance_pseudonormal(tree2,V,F,FN,VN,q2,s,sqrd,i,c2,n2);
         Eigen::RowVectorXd n;
-        (dim==3 ? n = n3 : n = n3);
+        (dim==3 ? n = n3 : n = n2);
         N.row(p) = n;
         break;
       }

include/igl/slice.cpp

@@ -225,22 +225,22 @@ IGL_INLINE void igl::slice(
 }
 
 template <typename DerivedX>
-IGL_INLINE Eigen::PlainObjectBase<DerivedX> igl::slice(
+IGL_INLINE DerivedX igl::slice(
   const Eigen::PlainObjectBase<DerivedX> & X,
   const Eigen::Matrix<int,Eigen::Dynamic,1> & R)
 {
-  Eigen::PlainObjectBase<DerivedX> Y;
+  DerivedX Y;
   igl::slice(X,R,Y);
   return Y;
 }
 
 template <typename DerivedX>
-IGL_INLINE Eigen::PlainObjectBase<DerivedX> igl::slice(
+IGL_INLINE DerivedX igl::slice(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
   const int dim)
 {
-  Eigen::PlainObjectBase<DerivedX> Y;
+  DerivedX Y;
   igl::slice(X,R,dim,Y);
   return Y;
 }
@@ -248,14 +248,14 @@ IGL_INLINE Eigen::PlainObjectBase<DerivedX> igl::slice(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
-template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
-template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
-template Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
-template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > igl::slice<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
-template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > igl::slice<Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::Matrix<double, -1, 1, 0, -1, 1>  igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::Matrix<double, -1, -1, 0, -1, -1>  igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::Matrix<double, -1, 1, 0, -1, 1>  igl::slice<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
+template Eigen::Matrix<double, -1, 3, 0, -1, 3>  igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::Matrix<double, 1, -1, 1, 1, -1>  igl::slice<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::Matrix<int, -1, 3, 0, -1, 3>  igl::slice<Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
+template Eigen::Matrix<int, -1, -1, 0, -1, -1 > igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template void igl::slice<std::complex<double>, std::complex<double> >(Eigen::SparseMatrix<std::complex<double>, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<std::complex<double>, 0, int>&);
-template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int);
 template void igl::slice<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::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::slice<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);

include/igl/slice.h

@@ -63,12 +63,17 @@ namespace igl
     const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
     Eigen::PlainObjectBase<DerivedY> & Y);
   // VectorXi Y = slice(X,R);
+  //
+  // This templating is bad because the return type might not have the same
+  // size as `DerivedX`. This will probably only work if DerivedX has Dynamic
+  // as it's non-trivial sizes or if the number of rows in R happens to equal
+  // the number of rows in `DerivedX`.
   template <typename DerivedX>
-  IGL_INLINE Eigen::PlainObjectBase<DerivedX> slice(
+  IGL_INLINE DerivedX slice(
     const Eigen::PlainObjectBase<DerivedX> & X,
     const Eigen::Matrix<int,Eigen::Dynamic,1> & R);
   template <typename DerivedX>
-  IGL_INLINE Eigen::PlainObjectBase<DerivedX> slice(
+  IGL_INLINE DerivedX slice(
     const Eigen::PlainObjectBase<DerivedX>& X,
     const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
     const int dim);

include/igl/slice_mask.cpp

@@ -94,23 +94,23 @@ IGL_INLINE void igl::slice_mask(
 }
 
 template <typename DerivedX>
-IGL_INLINE Eigen::PlainObjectBase<DerivedX> igl::slice_mask(
+IGL_INLINE DerivedX igl::slice_mask(
   const Eigen::PlainObjectBase<DerivedX> & X,
   const Eigen::Array<bool,Eigen::Dynamic,1> & R,
   const Eigen::Array<bool,Eigen::Dynamic,1> & C)
 {
-  Eigen::PlainObjectBase<DerivedX> Y;
+  DerivedX Y;
   igl::slice_mask(X,R,C,Y);
   return Y;
 }
 
 template <typename DerivedX>
-IGL_INLINE Eigen::PlainObjectBase<DerivedX> igl::slice_mask(
+IGL_INLINE DerivedX igl::slice_mask(
   const Eigen::PlainObjectBase<DerivedX>& X,
   const Eigen::Array<bool,Eigen::Dynamic,1> & R,
   const int dim)
 {
-  Eigen::PlainObjectBase<DerivedX> Y;
+  DerivedX Y;
   igl::slice_mask(X,R,dim,Y);
   return Y;
 }

include/igl/slice_mask.h

@@ -37,14 +37,18 @@ namespace igl
     const Eigen::Array<bool,Eigen::Dynamic,1> & R,
     const int dim,
     Eigen::PlainObjectBase<DerivedX> & Y);
-
+  //
+  // This templating is bad because the return type might not have the same
+  // size as `DerivedX`. This will probably only work if DerivedX has Dynamic
+  // as it's non-trivial sizes or if the number of rows in R happens to equal
+  // the number of rows in `DerivedX`.
   template <typename DerivedX>
-  IGL_INLINE Eigen::PlainObjectBase<DerivedX> slice_mask(
+  IGL_INLINE DerivedX slice_mask(
     const Eigen::PlainObjectBase<DerivedX> & X,
     const Eigen::Array<bool,Eigen::Dynamic,1> & R,
     const Eigen::Array<bool,Eigen::Dynamic,1> & C);
   template <typename DerivedX>
-  IGL_INLINE Eigen::PlainObjectBase<DerivedX> slice_mask(
+  IGL_INLINE DerivedX slice_mask(
     const Eigen::PlainObjectBase<DerivedX> & X,
     const Eigen::Array<bool,Eigen::Dynamic,1> & R,
     const int dim);

include/igl/sort.cpp

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

include/igl/triangle_triangle_adjacency.cpp

@@ -13,19 +13,40 @@
 #include <algorithm>
 #include <iostream>
 
-template <typename Scalar, typename Index>
-IGL_INLINE void igl::triangle_triangle_adjacency_preprocess(
-    const Eigen::PlainObjectBase<Scalar>& /*V*/,
-    const Eigen::PlainObjectBase<Index>& F,
-    std::vector<std::vector<int> >& TTT)
+// Extract the face adjacencies
+template <typename DerivedF, typename TTT_type, typename DerivedTT>
+IGL_INLINE void igl::triangle_triangle_adjacency_extractTT(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  std::vector<std::vector<TTT_type> >& TTT,
+  Eigen::PlainObjectBase<DerivedTT>& TT)
+{
+  TT.setConstant((int)(F.rows()),F.cols(),-1);
+
+  for(int i=1;i<(int)TTT.size();++i)
+  {
+    std::vector<int>& r1 = TTT[i-1];
+    std::vector<int>& r2 = TTT[i];
+    if ((r1[0] == r2[0]) && (r1[1] == r2[1]))
+    {
+      TT(r1[2],r1[3]) = r2[2];
+      TT(r2[2],r2[3]) = r1[2];
+    }
+  }
+}
+
+template <typename DerivedF, typename DerivedTT>
+IGL_INLINE void igl::triangle_triangle_adjacency(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedTT>& TT)
 {
-  return triangle_triangle_adjacency_preprocess(F,TTT);
+  DerivedTT TTi;
+  return triangle_triangle_adjacency(F,TT,TTi);
 }
 
-template <typename Index>
+template <typename DerivedF, typename TTT_type>
 IGL_INLINE void igl::triangle_triangle_adjacency_preprocess(
-    const Eigen::PlainObjectBase<Index>& F,
-    std::vector<std::vector<int> >& TTT)
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  std::vector<std::vector<TTT_type> >& TTT)
 {
   for(int f=0;f<F.rows();++f)
     for (int i=0;i<F.cols();++i)
@@ -42,32 +63,11 @@ IGL_INLINE void igl::triangle_triangle_adjacency_preprocess(
   std::sort(TTT.begin(),TTT.end());
 }
 
-// Extract the face adjacencies
-template <typename DerivedF, typename DerivedTT>
-IGL_INLINE void igl::triangle_triangle_adjacency_extractTT(
-  const Eigen::PlainObjectBase<DerivedF>& F,
-  std::vector<std::vector<int> >& TTT,
-  Eigen::PlainObjectBase<DerivedTT>& TT)
-{
-  TT.setConstant((int)(F.rows()),F.cols(),-1);
-
-  for(int i=1;i<(int)TTT.size();++i)
-  {
-    std::vector<int>& r1 = TTT[i-1];
-    std::vector<int>& r2 = TTT[i];
-    if ((r1[0] == r2[0]) && (r1[1] == r2[1]))
-    {
-      TT(r1[2],r1[3]) = r2[2];
-      TT(r2[2],r2[3]) = r1[2];
-    }
-  }
-}
-
 // Extract the face adjacencies indices (needed for fast traversal)
-template <typename DerivedF, typename DerivedTTi>
+template <typename DerivedF, typename TTT_type, typename DerivedTTi>
 IGL_INLINE void igl::triangle_triangle_adjacency_extractTTi(
   const Eigen::PlainObjectBase<DerivedF>& F,
-  std::vector<std::vector<int> >& TTT,
+  std::vector<std::vector<TTT_type> >& TTT,
   Eigen::PlainObjectBase<DerivedTTi>& TTi)
 {
   TTi.setConstant((int)(F.rows()),F.cols(),-1);
@@ -84,34 +84,6 @@ IGL_INLINE void igl::triangle_triangle_adjacency_extractTTi(
   }
 }
 
-// Compute triangle-triangle adjacency
-template <typename Scalar, typename Index>
-IGL_INLINE void igl::triangle_triangle_adjacency(const Eigen::PlainObjectBase<Scalar>& V,
-                        const Eigen::PlainObjectBase<Index>& F,
-                        Eigen::PlainObjectBase<Index>& TT)
-{
-  //assert(igl::is_edge_manifold(V,F));
-  std::vector<std::vector<int> > TTT;
-
-  triangle_triangle_adjacency_preprocess(V,F,TTT);
-  triangle_triangle_adjacency_extractTT(F,TTT,TT);
-}
-
-// Compute triangle-triangle adjacency with indices
-template <typename Scalar, typename Index>
-IGL_INLINE void igl::triangle_triangle_adjacency(
-  const Eigen::PlainObjectBase<Scalar>& /*V*/,
-  const Eigen::PlainObjectBase<Index>& F,
-  Eigen::PlainObjectBase<Index>& TT,
-  Eigen::PlainObjectBase<Index>& TTi)
-{
-  //assert(igl::is_edge_manifold(V,F));
-  std::vector<std::vector<int> > TTT;
-  triangle_triangle_adjacency_preprocess(F,TTT);
-  triangle_triangle_adjacency_extractTT(F,TTT,TT);
-  triangle_triangle_adjacency_extractTTi(F,TTT,TTi);
-}
-
 // Compute triangle-triangle adjacency with indices
 template <typename DerivedF, typename DerivedTT, typename DerivedTTi>
 IGL_INLINE void igl::triangle_triangle_adjacency(
@@ -119,7 +91,6 @@ IGL_INLINE void igl::triangle_triangle_adjacency(
   Eigen::PlainObjectBase<DerivedTT>& TT,
   Eigen::PlainObjectBase<DerivedTTi>& TTi)
 {
-  //assert(igl::is_edge_manifold(V,F));
   std::vector<std::vector<int> > TTT;
   triangle_triangle_adjacency_preprocess(F,TTT);
   triangle_triangle_adjacency_extractTT(F,TTT,TT);
@@ -235,10 +206,10 @@ template <
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::triangle_triangle_adjacency<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<int, -1, 3, 0, -1, 3> >&);
-template void igl::triangle_triangle_adjacency<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> >&);
-template void igl::triangle_triangle_adjacency<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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-template void igl::triangle_triangle_adjacency<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<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::triangle_triangle_adjacency<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> >&);
+// generated by autoexplicit.sh
+template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::triangle_triangle_adjacency<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> >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, long, long, long>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > const&, bool, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&, std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > >&);
 template void igl::triangle_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, int, int>(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, bool, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&, std::vector<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >, std::allocator<std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > > >&);

include/igl/triangle_triangle_adjacency.h

@@ -22,7 +22,6 @@ namespace igl
   //   Index  derived type of eigen matrix for F (e.g. derived from
   //     MatrixXi)
   // Inputs:
-  //   V  #V by dim list of mesh vertex positions
   //   F  #F by simplex_size list of mesh faces (must be triangles)
   // Outputs:
   //   TT   #F by #3 adjacent matrix, the element i,j is the id of the triangle adjacent to the j edge of triangle i
@@ -30,56 +29,40 @@ namespace igl
   // NOTE: the first edge of a triangle is [0,1] the second [1,2] and the third [2,3].
   //       this convention is DIFFERENT from cotmatrix_entries.h
   // Known bug: this should not need to take V as input.
-
-  template <typename Scalar, typename Index>
-  IGL_INLINE void triangle_triangle_adjacency(
-    const Eigen::PlainObjectBase<Scalar>& V,
-    const Eigen::PlainObjectBase<Index>& F,
-    Eigen::PlainObjectBase<Index>& TT);
-  // Compute triangle-triangle adjacency with indices
-  template <typename Scalar, typename Index>
-  IGL_INLINE void triangle_triangle_adjacency(
-    const Eigen::PlainObjectBase<Scalar>& V,
-    const Eigen::PlainObjectBase<Index>& F,
-    Eigen::PlainObjectBase<Index>& TT,
-    Eigen::PlainObjectBase<Index>& TTi);
-
   template <typename DerivedF, typename DerivedTT, typename DerivedTTi>
   IGL_INLINE void triangle_triangle_adjacency(
     const Eigen::PlainObjectBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedTT>& TT,
     Eigen::PlainObjectBase<DerivedTTi>& TTi);
-
-
+  template <typename DerivedF, typename DerivedTT>
+  IGL_INLINE void triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedTT>& TT);
   // Preprocessing
-  template <typename Scalar, typename Index>
-  IGL_INLINE void triangle_triangle_adjacency_preprocess(
-    const Eigen::PlainObjectBase<Scalar>& V,
-    const Eigen::PlainObjectBase<Index>& F,
-    std::vector<std::vector<int> >& TTT);
-  template <typename DerivedF>
+  template <typename DerivedF, typename TTT_type>
   IGL_INLINE void triangle_triangle_adjacency_preprocess(
     const Eigen::PlainObjectBase<DerivedF>& F,
-    std::vector<std::vector<int> >& TTT);
+    std::vector<std::vector<TTT_type> >& TTT);
   // Extract the face adjacencies
-  template <typename DerivedF, typename DerivedTT>
+  template <typename DerivedF, typename TTT_type, typename DerivedTT>
   IGL_INLINE void triangle_triangle_adjacency_extractTT(
     const Eigen::PlainObjectBase<DerivedF>& F,
-    std::vector<std::vector<int> >& TTT,
+    std::vector<std::vector<TTT_type> >& TTT,
     Eigen::PlainObjectBase<DerivedTT>& TT);
   // Extract the face adjacencies indices (needed for fast traversal)
-  template <typename DerivedF, typename DerivedTTi>
+  template <typename DerivedF, typename TTT_type, typename DerivedTTi>
   IGL_INLINE void triangle_triangle_adjacency_extractTTi(
     const Eigen::PlainObjectBase<DerivedF>& F,
-    std::vector<std::vector<int> >& TTT,
+    std::vector<std::vector<TTT_type> >& TTT,
     Eigen::PlainObjectBase<DerivedTTi>& TTi);
   // Adjacency list version, which works with non-manifold meshes
   //
   // Inputs:
   //   F  #F by 3 list of triangle indices
   // Outputs:
-  //   TT  #F by 3 list of lists so that TT[i][c] --> {j,k,...} means that faces j and
-  //     k etc. are edge-neighbors of face i on face i's edge opposite corner c
+  //   TT  #F by 3 list of lists so that TT[i][c] --> {j,k,...} means that
+  //     faces j and k etc. are edge-neighbors of face i on face i's edge
+  //     opposite corner c
   //   TTj  #F list of lists so that TTj[i][c] --> {j,k,...} means that face
   //     TT[i][c][0] is an edge-neighbor of face i incident on the edge of face
   //     TT[i][c][0] opposite corner j, and TT[i][c][1] " corner k, etc.

include/igl/unique.cpp

@@ -209,7 +209,7 @@ IGL_INLINE void igl::unique_rows(
   using namespace std;
   using namespace Eigen;
   VectorXi IM;
-  Eigen::PlainObjectBase<DerivedA> sortA;
+  DerivedA sortA;
   sortrows(A,true,sortA,IM);
 
 

include/igl/unproject.cpp

@@ -10,14 +10,22 @@
 #include <Eigen/Dense>
 #include <Eigen/LU>
 
-template <typename Scalar>
-IGL_INLINE Eigen::Matrix<Scalar,3,1> igl::unproject(
-  const    Eigen::Matrix<Scalar,3,1>&  win,
-  const    Eigen::Matrix<Scalar,4,4>& model,
-  const    Eigen::Matrix<Scalar,4,4>& proj,
-  const    Eigen::Matrix<Scalar,4,1>&  viewport)
+template <
+  typename Derivedwin,
+  typename Derivedmodel,
+  typename Derivedproj,
+  typename Derivedviewport,
+  typename Derivedscene>
+IGL_INLINE void igl::unproject(
+  const Eigen::PlainObjectBase<Derivedwin>&  win,
+  const Eigen::PlainObjectBase<Derivedmodel>& model,
+  const Eigen::PlainObjectBase<Derivedproj>& proj,
+  const Eigen::PlainObjectBase<Derivedviewport>&  viewport,
+  Eigen::PlainObjectBase<Derivedscene> & scene)
 {
-  Eigen::Matrix<Scalar,4,4> Inverse = (proj * model).inverse();
+  typedef typename Derivedscene::Scalar Scalar;
+  Eigen::Matrix<Scalar,4,4> Inverse = 
+    (proj.template cast<Scalar>() * model.template cast<Scalar>()).inverse();
 
   Eigen::Matrix<Scalar,4,1> tmp;
   tmp << win, 1;
@@ -28,7 +36,19 @@ IGL_INLINE Eigen::Matrix<Scalar,3,1> igl::unproject(
   Eigen::Matrix<Scalar,4,1> obj = Inverse * tmp;
   obj /= obj(3);
 
-  return obj.head(3);
+  scene = obj.head(3);
+}
+
+template <typename Scalar>
+IGL_INLINE Eigen::Matrix<Scalar,3,1> igl::unproject(
+  const    Eigen::Matrix<Scalar,3,1>&  win,
+  const    Eigen::Matrix<Scalar,4,4>& model,
+  const    Eigen::Matrix<Scalar,4,4>& proj,
+  const    Eigen::Matrix<Scalar,4,1>&  viewport)
+{
+  Eigen::Matrix<Scalar,3,1> scene;
+  unproject(win,model,proj,viewport,scene);
+  return scene;
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/unproject.h

@@ -12,17 +12,32 @@
 namespace igl
 {
   // Eigen reimplementation of gluUnproject
+  //
   // Inputs:
   //   win  screen space x, y, and z coordinates
-  // Returns:
-  //   the unprojected x, y, and z coordinates
-  // Returns return value of gluUnProject call
+  //   model  4x4 model-view matrix
+  //   proj  4x4 projection matrix
+  //   viewport  4-long viewport vector
+  // Outputs:
+  //   scene  the unprojected x, y, and z coordinates
+  template <
+    typename Derivedwin,
+    typename Derivedmodel,
+    typename Derivedproj,
+    typename Derivedviewport,
+    typename Derivedscene>
+  IGL_INLINE void unproject(
+    const Eigen::PlainObjectBase<Derivedwin>&  win,
+    const Eigen::PlainObjectBase<Derivedmodel>& model,
+    const Eigen::PlainObjectBase<Derivedproj>& proj,
+    const Eigen::PlainObjectBase<Derivedviewport>&  viewport,
+    Eigen::PlainObjectBase<Derivedscene> & scene);
   template <typename Scalar>
   IGL_INLINE Eigen::Matrix<Scalar,3,1> unproject(
-    const    Eigen::Matrix<Scalar,3,1>&  win,
-    const    Eigen::Matrix<Scalar,4,4>& model,
-    const    Eigen::Matrix<Scalar,4,4>& proj,
-    const    Eigen::Matrix<Scalar,4,1>&  viewport);
+    const Eigen::Matrix<Scalar,3,1>&  win,
+    const Eigen::Matrix<Scalar,4,4>& model,
+    const Eigen::Matrix<Scalar,4,4>& proj,
+    const Eigen::Matrix<Scalar,4,1>&  viewport);
 }
 
 

include/igl/unproject_in_mesh.cpp

@@ -0,0 +1,98 @@
+// 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
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "unproject_in_mesh.h"
+#include "unproject_ray.h"
+#include "ray_mesh_intersect.h"
+
+template < typename Derivedobj>
+  IGL_INLINE int igl::unproject_in_mesh(
+      const Eigen::Vector2f& pos,
+      const Eigen::Matrix4f& model,
+      const Eigen::Matrix4f& proj,
+      const Eigen::Vector4f& viewport,
+      const std::function<
+        void(
+          const Eigen::Vector3f&,
+          const Eigen::Vector3f&,
+          std::vector<igl::Hit> &)
+          > & shoot_ray,
+      Eigen::PlainObjectBase<Derivedobj> & obj,
+      std::vector<igl::Hit > & hits)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  Vector3f s,dir;
+  unproject_ray(pos,model,proj,viewport,s,dir);
+  shoot_ray(s,dir,hits);
+  switch(hits.size())
+  {
+    case 0:
+      break;
+    case 1:
+    {
+      obj = (s + dir*hits[0].t).cast<typename Derivedobj::Scalar>();
+      break;
+    }
+    case 2:
+    default:
+    {
+      obj = 0.5*((s + dir*hits[0].t) + (s + dir*hits[1].t)).cast<typename Derivedobj::Scalar>();
+      break;
+    }
+  }
+  return hits.size();
+}
+
+extern "C"
+{
+#include "raytri.c"
+}
+
+template < typename DerivedV, typename DerivedF, typename Derivedobj>
+  IGL_INLINE int igl::unproject_in_mesh(
+      const Eigen::Vector2f& pos,
+      const Eigen::Matrix4f& model,
+      const Eigen::Matrix4f& proj,
+      const Eigen::Vector4f& viewport,
+      const Eigen::PlainObjectBase<DerivedV> & V,
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      Eigen::PlainObjectBase<Derivedobj> & obj,
+      std::vector<igl::Hit > & hits)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  const auto & shoot_ray = [&V,&F](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir,
+    std::vector<igl::Hit> & hits)
+  {
+    ray_mesh_intersect(s,dir,V,F,hits);
+  };
+  return unproject_in_mesh(pos,model,proj,viewport,shoot_ray,obj,hits);
+}
+
+template < typename DerivedV, typename DerivedF, typename Derivedobj>
+  IGL_INLINE int igl::unproject_in_mesh(
+      const Eigen::Vector2f& pos,
+      const Eigen::Matrix4f& model,
+      const Eigen::Matrix4f& proj,
+      const Eigen::Vector4f& viewport,
+      const Eigen::PlainObjectBase<DerivedV> & V,
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      Eigen::PlainObjectBase<Derivedobj> & obj)
+{
+  std::vector<igl::Hit> hits;
+  return unproject_in_mesh(pos,model,proj,viewport,V,F,obj,hits);
+}
+#ifdef IGL_STATIC_LIBRARY
+template int igl::unproject_in_mesh<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, std::function<void (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+template int igl::unproject_in_mesh<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, std::function<void (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+template int igl::unproject_in_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, std::function<void (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+#endif

include/igl/unproject_in_mesh.h

@@ -0,0 +1,87 @@
+// 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
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_UNPROJECT_IN_MESH
+#define IGL_UNPROJECT_IN_MESH
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+#include <vector>
+#include "Hit.h"
+
+namespace igl
+{
+  // Unproject a screen location (using current opengl viewport, projection, and
+  // model view) to a 3D position _inside_ a given mesh. If the ray through the
+  // given screen location (x,y) _hits_ the mesh more than twice then the 3D
+  // midpoint between the first two hits is return. If it hits once, then that
+  // point is return. If it does not hit the mesh then obj is not set.
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    V   #V by 3 list of mesh vertex positions
+  //    F   #F by 3 list of mesh triangle indices into V
+  // Outputs:
+  //    obj        3d unprojected mouse point in mesh
+  //    hits       vector of hits
+  // Returns number of hits
+  //
+  template < typename DerivedV, typename DerivedF, typename Derivedobj>
+    IGL_INLINE int unproject_in_mesh(
+        const Eigen::Vector2f& pos,
+        const Eigen::Matrix4f& model,
+        const Eigen::Matrix4f& proj,
+        const Eigen::Vector4f& viewport,
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<Derivedobj> & obj,
+        std::vector<igl::Hit > & hits);
+  //
+  // Inputs:
+  //    pos        screen space coordinates
+  //    model      model matrix
+  //    proj       projection matrix
+  //    viewport   vieweport vector
+  //    shoot_ray  function handle that outputs hits of a given ray against a
+  //      mesh (embedded in function handles as captured variable/data)
+  // Outputs:
+  //    obj        3d unprojected mouse point in mesh
+  //    hits       vector of hits
+  // Returns number of hits
+  //
+  template < typename Derivedobj>
+    IGL_INLINE int unproject_in_mesh(
+        const Eigen::Vector2f& pos,
+        const Eigen::Matrix4f& model,
+        const Eigen::Matrix4f& proj,
+        const Eigen::Vector4f& viewport,
+        const std::function<
+          void(
+            const Eigen::Vector3f&,
+            const Eigen::Vector3f&,
+            std::vector<igl::Hit> &)
+            > & shoot_ray,
+        Eigen::PlainObjectBase<Derivedobj> & obj,
+        std::vector<igl::Hit > & hits);
+  template < typename DerivedV, typename DerivedF, typename Derivedobj>
+    IGL_INLINE int unproject_in_mesh(
+        const Eigen::Vector2f& pos,
+        const Eigen::Matrix4f& model,
+        const Eigen::Matrix4f& proj,
+        const Eigen::Vector4f& viewport,
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<Derivedobj> & obj);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "unproject_in_mesh.cpp"
+#endif
+#endif
+

include/igl/unproject_ray.cpp

@@ -0,0 +1,42 @@
+// 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
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "unproject_ray.h"
+#include "unproject.h"
+
+template <
+  typename Derivedpos,
+  typename Derivedmodel,
+  typename Derivedproj,
+  typename Derivedviewport,
+  typename Deriveds,
+  typename Deriveddir>
+IGL_INLINE void igl::unproject_ray(
+  const Eigen::PlainObjectBase<Derivedpos> & pos,
+  const Eigen::PlainObjectBase<Derivedmodel> & model,
+  const Eigen::PlainObjectBase<Derivedproj> & proj,
+  const Eigen::PlainObjectBase<Derivedviewport> & viewport,
+  Eigen::PlainObjectBase<Deriveds> & s,
+  Eigen::PlainObjectBase<Deriveddir> & dir)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  // Source and direction on screen
+  typedef Eigen::Matrix<typename Deriveds::Scalar,3,1> Vec3;
+  Vec3 win_s(pos(0),pos(1),0);
+  Vec3 win_d(pos(0),pos(1),1);
+  // Source, destination and direction in world
+  Vec3 d;
+  igl::unproject(win_s,model,proj,viewport,s);
+  igl::unproject(win_d,model,proj,viewport,d);
+  dir = d-s;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::unproject_ray<Eigen::Matrix<float, 2, 1, 0, 2, 1>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 4, 0, 4, 4>, Eigen::Matrix<float, 4, 1, 0, 4, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 2, 1, 0, 2, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 4, 0, 4, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 4, 1, 0, 4, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
+#endif