Merge branch 'master' of https://github.com/libigl/libigl

Former-commit-id: 1ae1ca38c38922adffbb562d906f0c40c20ddfa0

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+1.0.3  Bone heat method
 1.0.2  Bug fix in winding number code
 1.0.1  Bug fixes and more CGAL support
 1.0.0  Major beta release: many renames, tutorial, triangle wrapper, org. build

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-1.0.2
+1.0.3

include/igl/WindingNumberAABB.h

@@ -31,6 +31,7 @@ namespace igl
         NUM_SPLIT_METHODS = 3
       } split_method;
     public:
+      inline WindingNumberAABB(){}
       inline WindingNumberAABB(
         const Eigen::MatrixXd & V,
         const Eigen::MatrixXi & F);
@@ -38,6 +39,9 @@ namespace igl
         const WindingNumberTree<Point> & parent,
         const Eigen::MatrixXi & F);
       // Initialize some things
+      inline void set_mesh(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
       inline void init();
       inline bool inside(const Point & p) const;
       inline virtual void grow();
@@ -67,6 +71,15 @@ namespace igl
 #  define WindingNumberAABB_MIN_F 100
 #endif
 
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::set_mesh(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F)
+{
+  igl::WindingNumberTree<Point>::set_mesh(V,F);
+  init();
+}
+
 template <typename Point>
 inline void igl::WindingNumberAABB<Point>::init()
 {
@@ -142,7 +155,6 @@ inline void igl::WindingNumberAABB<Point>::grow()
 
   // Blerg, why is selecting rows so difficult
 
-  vector<int> id( this->getF().rows());
   double split_value;
   // Split in longest direction
   switch(split_method)
@@ -160,6 +172,7 @@ inline void igl::WindingNumberAABB<Point>::grow()
   //cout<<"c: "<<0.5*(max_corner[max_d] + min_corner[max_d])<<" "<<
   //  "m: "<<split_value<<endl;;
 
+  vector<int> id( this->getF().rows());
   for(int i = 0;i<this->getF().rows();i++)
   {
     if(BC(i,max_d) <= split_value)
@@ -216,7 +229,7 @@ inline bool igl::WindingNumberAABB<Point>::inside(const Point & p) const
   {
     //// Perfect matching is **not** robust
     //if( p(i) < min_corner(i) || p(i) >= max_corner(i))
-    // **MUST** be conservative!!
+    // **MUST** be conservative
     if( p(i) < min_corner(i) || p(i) > max_corner(i))
     {
       return false;

include/igl/WindingNumberTree.h

@@ -12,6 +12,7 @@
 #include <Eigen/Dense>
 #include "WindingNumberMethod.h"
 
+static Eigen::MatrixXd dummyV;
 namespace igl
 {
   // Templates:
@@ -44,6 +45,7 @@ namespace igl
       // (Approximate) center (of mass)
       Point center;
     public:
+      inline WindingNumberTree():V(dummyV){}
       // For root
       inline WindingNumberTree(
         const Eigen::MatrixXd & V,
@@ -53,6 +55,9 @@ namespace igl
         const WindingNumberTree<Point> & parent,
         const Eigen::MatrixXi & F);
       inline virtual ~WindingNumberTree();
+      inline virtual void set_mesh(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
       // Set method
       inline void set_method( const WindingNumberMethod & m);
     public:
@@ -137,7 +142,6 @@ template <typename Point>
 std::map< std::pair<const igl::WindingNumberTree<Point>*,const igl::WindingNumberTree<Point>*>, double>
   igl::WindingNumberTree<Point>::cached;
 
-static Eigen::MatrixXd dummyV;
 template <typename Point>
 inline igl::WindingNumberTree<Point>::WindingNumberTree(
   const Eigen::MatrixXd & _V,
@@ -152,6 +156,15 @@ inline igl::WindingNumberTree<Point>::WindingNumberTree(
   radius(std::numeric_limits<double>::infinity()),
   center(0,0,0)
 {
+  set_mesh(_V,_F);
+}
+
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::set_mesh(
+    const Eigen::MatrixXd & _V,
+    const Eigen::MatrixXi & _F)
+{
+  using namespace std;
   // Remove any exactly duplicate vertices
   // Q: Can this ever increase the complexity of the boundary?
   // Q: Would we gain even more by remove almost exactly duplicate vertices?

include/igl/bounding_box.cpp

@@ -0,0 +1,85 @@
+// 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 "bounding_box.h"
+#include <iostream>
+
+template <typename DerivedV, typename DerivedBV, typename DerivedBF>
+IGL_INLINE void igl::bounding_box(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedBV>& BV,
+  Eigen::PlainObjectBase<DerivedBF>& BF)
+{
+  using namespace std;
+
+  const int dim = V.cols();
+  const auto & minV = V.colwise().minCoeff();
+  const auto & maxV = V.colwise().maxCoeff();
+  // 2^n vertices
+  BV.resize((1<<dim),dim);
+
+  // Recursive lambda to generate all 2^n combinations
+  const std::function<void(const int,const int,int*,int)> combos =
+  [&BV,&minV,&maxV,&combos](
+    const int dim,
+    const int i,
+    int * X,
+    const int pre_index)
+  {
+    for(X[i] = 0;X[i]<2;X[i]++)
+    {
+      int index = pre_index*2+X[i];
+      if((i+1)<dim)
+      {
+        combos(dim,i+1,X,index);
+      }else
+      {
+        for(int d = 0;d<dim;d++)
+        {
+          BV(index,d) = (X[d]?minV[d]:maxV[d]);
+        }
+      }
+    }
+  };
+
+  Eigen::VectorXi X(dim);
+  combos(dim,0,X.data(),0);
+  switch(dim)
+  {
+    case 2:
+      BF.resize(4,2);
+      BF<<
+        3,1,
+        1,0,
+        0,2,
+        2,3;
+      break;
+    case 3:
+      BF.resize(12,3);
+      BF<<
+        2,0,6,
+        0,4,6,
+        5,4,0,
+        5,0,1,
+        6,4,5,
+        5,7,6,
+        3,0,2,
+        1,0,3,
+        3,2,6,
+        6,7,3,
+        5,1,3,
+        3,7,5;
+      break;
+    default:
+      assert(false && "Unsupported dimension.");
+      break;
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/bounding_box.h

@@ -0,0 +1,33 @@
+// 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_BOUNDING_BOX_H
+#define IGL_BOUNDING_BOX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Build a triangle mesh of the bounding box of a given list of vertices
+  // 
+  // Inputs:
+  //   V  #V by dim list of rest domain positions
+  // Outputs:
+  //   BV  2^dim by dim list of bounding box corners positions
+  //   BF  #BF by dim list of simplex facets 
+  template <typename DerivedV, typename DerivedBV, typename DerivedBF>
+  IGL_INLINE void bounding_box(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedBV>& BV,
+    Eigen::PlainObjectBase<DerivedBF>& BF);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bounding_box.cpp"
+#endif
+
+#endif
+

include/igl/cgal/signed_distance.cpp

@@ -6,6 +6,100 @@
 // 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 "signed_distance.h"
+#include "../per_vertex_normals.h"
+#include "../per_edge_normals.h"
+#include "../per_face_normals.h"
+#include "point_mesh_squared_distance.h"
+
+
+template <typename Kernel>
+IGL_INLINE void igl::signed_distance(
+  const Eigen::MatrixXd & P,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const SignedDistanceType sign_type,
+  Eigen::VectorXd & S,
+  Eigen::VectorXi & I,
+  Eigen::MatrixXd & C,
+  Eigen::MatrixXd & N)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(V.cols() == 3 && "V should have 3d positions");
+  assert(P.cols() == 3 && "P should have 3d positions");
+  assert(F.cols() == 3 && "F should have triangles");
+
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  // Prepare distance computation
+  Tree tree;
+  vector<Triangle_3 > T;
+  point_mesh_squared_distance_precompute(V,F,tree,T);
+
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  WindingNumberAABB<Eigen::Vector3d> hier;
+  switch(sign_type)
+  {
+    default:
+      assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_DEFAULT:
+    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+      hier.set_mesh(V,F);
+      hier.grow();
+      break;
+    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+      // [Bærentzen & Aanæs 2005]
+      per_face_normals(V,F,FN);
+      per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
+      per_edge_normals(
+        V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,EN,E,EMAP);
+      N.resize(P.rows(),3);
+      break;
+  }
+
+  S.resize(P.rows(),1);
+  I.resize(P.rows(),1);
+  C.resize(P.rows(),3);
+  for(int p = 0;p<P.rows();p++)
+  {
+    const Point_3 q(P(p,0),P(p,1),P(p,2));
+    double s,sqrd;
+    Point_and_primitive_id pp;
+    switch(sign_type)
+    {
+      default:
+        assert(false && "Unknown SignedDistanceType");
+      case SIGNED_DISTANCE_TYPE_DEFAULT:
+      case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+        signed_distance_winding_number(tree,hier,q,s,sqrd,pp);
+        break;
+      case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      {
+        Vector3d n(0,0,0);
+        signed_distance_pseudonormal(tree,T,F,FN,VN,EN,EMAP,q,s,sqrd,pp,n);
+        N.row(p) = n;
+        break;
+      }
+    }
+    I(p) = pp.second - T.begin();
+    S(p) = s*sqrt(sqrd);
+    C(p,0) = pp.first.x();
+    C(p,1) = pp.first.y();
+    C(p,2) = pp.first.z();
+  }
+}
+
+
 template <typename Kernel>
 IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   const CGAL::AABB_tree<
@@ -22,6 +116,44 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   const Eigen::MatrixXd & EN,
   const Eigen::VectorXi & EMAP,
   const typename Kernel::Point_3 & q)
+{
+  typename CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
+        > 
+      >::Point_and_primitive_id pp;
+  typename Kernel::FT s,sqrd;
+  Eigen::Vector3d n(0,0,0);
+  signed_distance_pseudonormal<Kernel>(tree,T,F,FN,VN,EN,EMAP,q,s,sqrd,pp,n);
+  return s*sqrt(sqrd);
+}
+
+template <typename Kernel>
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::MatrixXd & EN,
+  const Eigen::VectorXi & EMAP,
+  const typename Kernel::Point_3 & q,
+  typename Kernel::FT & s,
+  typename Kernel::FT & sqrd,
+  typename CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
+        > 
+      >::Point_and_primitive_id & pp,
+   Eigen::Vector3d & n)
 {
   using namespace Eigen;
   using namespace std;
@@ -34,10 +166,10 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
   typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
 
-  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  pp = tree.closest_point_and_primitive(q);
   Point_3 & p = pp.first;
   const auto & qp = q-p;
-  const FT sqrd = qp.squared_length();
+  sqrd = qp.squared_length();
   Vector3d v(qp.x(),qp.y(),qp.z());
   const int f = pp.second - T.begin();
   const Triangle_3 & t = *pp.second;
@@ -49,7 +181,6 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
   Vector3d b(area(1,2),area(2,0),area(0,1));
   b /= b.sum();
   // Determine which normal to use
-  Vector3d n;
   const double epsilon = 1e-12;
   const int type = (b.array()<=epsilon).cast<int>().sum();
   switch(type)
@@ -82,8 +213,7 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
       n = FN.row(f);
       break;
   }
-  const double s = (v.dot(n) >= 0 ? 1. : -1.);
-  return s*sqrt(sqrd);
+  s = (v.dot(n) >= 0 ? 1. : -1.);
 }
 
 template <typename Kernel>
@@ -97,6 +227,38 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
   > & tree,
   const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
   const typename Kernel::Point_3 & q)
+{
+  typename CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
+        > 
+      >::Point_and_primitive_id pp;
+  typename Kernel::FT s,sqrd;
+  signed_distance_winding_number<Kernel>(tree,hier,q,s,sqrd,pp);
+  return s*sqrt(sqrd);
+}
+
+
+template <typename Kernel>
+IGL_INLINE void igl::signed_distance_winding_number(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+  const typename Kernel::Point_3 & q,
+  typename Kernel::FT & s,
+  typename Kernel::FT & sqrd,
+  typename CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
+        > 
+      >::Point_and_primitive_id & pp)
 {
   typedef typename Kernel::FT FT;
   typedef typename Kernel::Point_3 Point_3;
@@ -109,14 +271,13 @@ IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
   using namespace Eigen;
   using namespace std;
   using namespace igl;
-  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  pp = tree.closest_point_and_primitive(q);
   Point_3 & p = pp.first;
   const auto & qp = q-p;
   const Vector3d eq(q.x(),q.y(),q.z()); 
-  const FT sqrd = qp.squared_length();
+  sqrd = qp.squared_length();
   const double w = hier.winding_number(eq);
-  const FT s = 1.-2.*w;
-  return s*sqrt(sqrd);
+  s = 1.-2.*w;
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/cgal/signed_distance.h

@@ -19,8 +19,38 @@ namespace igl
     SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
     SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
     SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
-    NUM_SIGNED_DISTANCE_TYPE            = 3
+    SIGNED_DISTANCE_TYPE_IGNORE         = 3,
+    NUM_SIGNED_DISTANCE_TYPE            = 4
   };
+  // Computes signed distance to a mesh
+  //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Simple_cartesian<double>)
+  // Inputs:
+  //   P  #P by 3 list of query point positions
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
+  // Outputs:
+  //   S  #P list of smallest signed distances
+  //   I  #P list of facet indices corresponding to smallest distances
+  //   C  #P by 3 list of closest points
+  //   N  #P by 3 list of closest normals (only set if
+  //     sign_type=SIGNED_DISTANCE_TYPE_PSEUDONORMAL)
+  //
+  // Known bugs: This only computes distances to triangles. So unreferenced
+  // vertices are ignored.
+  template <typename Kernel>
+  IGL_INLINE void signed_distance(
+    const Eigen::MatrixXd & P,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const SignedDistanceType sign_type,
+    Eigen::VectorXd & S,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C,
+    Eigen::MatrixXd & N);
   // Computes signed distance to mesh
   //
   // Templates:
@@ -53,6 +83,37 @@ namespace igl
     const Eigen::MatrixXd & EN,
     const Eigen::VectorXi & EMAP,
     const typename Kernel::Point_3 & q);
+  // Outputs:
+  //   s  sign
+  //   sqrd  squared distance
+  //   pp  closest point and primitve
+  //   n  normal
+  template <typename Kernel>
+  IGL_INLINE void signed_distance_pseudonormal(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const std::vector<CGAL::Triangle_3<Kernel> > & T,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    const typename Kernel::Point_3 & q,
+    typename Kernel::FT & s,
+    typename Kernel::FT & sqrd,
+    typename CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
+          > 
+        >::Point_and_primitive_id & pp,
+   Eigen::Vector3d & n);
+
   // Inputs:
   //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
   //   hier  Winding number evaluation hierarchy
@@ -69,6 +130,29 @@ namespace igl
     > & tree,
     const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
     const typename Kernel::Point_3 & q);
+  // Outputs:
+  //   s  sign
+  //   sqrd  squared distance
+  //   pp  closest point and primitve
+  template <typename Kernel>
+  IGL_INLINE void signed_distance_winding_number(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+    const typename Kernel::Point_3 & q,
+    typename Kernel::FT & s,
+    typename Kernel::FT & sqrd,
+    typename CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator> 
+          > 
+        >::Point_and_primitive_id & pp);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/cgal/signed_distance_isosurface.cpp

@@ -62,20 +62,33 @@ IGL_INLINE bool igl::signed_distance_isosurface(
   typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
   typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
 
-  Eigen::MatrixXd FN,VN,EN;
-  Eigen::MatrixXi E;
-  Eigen::VectorXi EMAP;
-  // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
-  // [Bærentzen & Aanæs 2005]
-  per_face_normals(IV,IF,FN);
-  per_vertex_normals(IV,IF,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,VN);
-  per_edge_normals(IV,IF,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,EN,E,EMAP);
-  // Prepare distance computation
   Tree tree;
   vector<Triangle_3 > T;
   point_mesh_squared_distance_precompute(IV,IF,tree,T);
-  WindingNumberAABB<Eigen::Vector3d> hier(IV,IF);
-  hier.grow();
+
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  WindingNumberAABB<Eigen::Vector3d> hier;
+  switch(sign_type)
+  {
+    default:
+      assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_DEFAULT:
+    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+      hier.set_mesh(IV,IF);
+      hier.grow();
+      break;
+    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+      // [Bærentzen & Aanæs 2005]
+      per_face_normals(IV,IF,FN);
+      per_vertex_normals(IV,IF,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
+      per_edge_normals(
+        IV,IF,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
+      break;
+  }
+
   Tr tr;            // 3D-Delaunay triangulation
   C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
   // defining the surface

include/igl/cotmatrix_entries.h

@@ -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) 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 

include/igl/doublearea.cpp

@@ -25,18 +25,30 @@ IGL_INLINE void igl::doublearea(
   Eigen::PlainObjectBase<DerivedV> l;
   // "Lecture Notes on Geometric Robustness" Shewchuck 09, Section 3.1
   // http://www.cs.berkeley.edu/~jrs/meshpapers/robnotes.pdf
+
+  // Projected area helper
+  const auto & proj_doublearea = 
+    [&V,&F](const int x, const int y, const int f)->double
+  {
+    auto rx = V(F(f,0),x)-V(F(f,2),x);
+    auto sx = V(F(f,1),x)-V(F(f,2),x);
+    auto ry = V(F(f,0),y)-V(F(f,2),y);
+    auto sy = V(F(f,1),y)-V(F(f,2),y);
+    return rx*sy - ry*sx;
+  };
+
   switch(dim)
   {
     case 3:
     {
       dblA = Eigen::PlainObjectBase<DeriveddblA>::Zero(m,1);
-      for(int d = 0;d<3;d++)
+      for(int f = 0;f<F.rows();f++)
       {
-        Eigen::Matrix<typename DerivedV::Scalar, DerivedV::RowsAtCompileTime, 2> Vd(V.rows(),2);
-        Vd << V.col(d), V.col((d+1)%3);
-        Eigen::PlainObjectBase<DeriveddblA> dblAd;
-        doublearea(Vd,F,dblAd);
-        dblA.array() += dblAd.array().square();
+        for(int d = 0;d<3;d++)
+        {
+          double dblAd = proj_doublearea(d,(d+1)%3,f);
+          dblA[f] += dblAd*dblAd;
+        }
       }
       dblA = dblA.array().sqrt().eval();
       break;
@@ -46,9 +58,7 @@ IGL_INLINE void igl::doublearea(
       dblA.resize(m,1);
       for(int f = 0;f<m;f++)
       {
-        auto r = V.row(F(f,0))-V.row(F(f,2));
-        auto s = V.row(F(f,1))-V.row(F(f,2));
-        dblA(f) = r(0)*s(1) - r(1)*s(0);
+        dblA(f) = proj_doublearea(0,1,f);
       }
       break;
     }

include/igl/doublearea.h

@@ -26,6 +26,8 @@ namespace igl
   // Outputs:
   //   dblA  #F list of triangle double areas (SIGNED only for 2D input)
   //
+  // Known bug: For dim==3 complexity is O(#V + #F)!! Not just O(#F). This is a big deal
+  // if you have 1million unreferenced vertices and 1 face
   template <typename DerivedV, typename DerivedF, typename DeriveddblA>
   IGL_INLINE void doublearea( 
     const Eigen::PlainObjectBase<DerivedV> & V, 

include/igl/embree/bone_heat.cpp

@@ -0,0 +1,116 @@
+// 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 "bone_heat.h"
+#include "EmbreeIntersector.h"
+#include "bone_visible.h"
+#include "../project_to_line_segment.h"
+#include "../cotmatrix.h"
+#include "../massmatrix.h"
+#include "../mat_min.h"
+#include <Eigen/Sparse>
+
+bool igl::bone_heat(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & C,
+  const Eigen::VectorXi & P,
+  const Eigen::MatrixXi & BE,
+  const Eigen::MatrixXi & CE,
+  Eigen::MatrixXd & W)
+{
+  using namespace std;
+  using namespace Eigen;
+  assert(CE.rows() == 0 && "Cage edges not supported.");
+  assert(C.cols() == V.cols() && "V and C should have same #cols");
+  assert(BE.cols() == 2 && "BE should have #cols=2");
+  assert(F.cols() == 3 && "F should contain triangles.");
+  assert(V.cols() == 3 && "V should contain 3D positions.");
+
+  const int n = V.rows();
+  const int np = P.rows();
+  const int nb = BE.rows();
+  const int m = np + nb;
+
+  // "double sided lighting"
+  MatrixXi FF;
+  FF.resize(F.rows()*2,F.cols());
+  FF << F, F.rowwise().reverse();
+  // Initialize intersector
+  EmbreeIntersector ei;
+  ei.init(V.cast<float>(),F.cast<int>());
+
+  typedef Matrix<bool,Dynamic,1> VectorXb;
+  typedef Matrix<bool,Dynamic,Dynamic> MatrixXb;
+  MatrixXb vis_mask(n,m);
+  // Distances
+  MatrixXd D(n,m);
+  // loop over points
+  for(int j = 0;j<np;j++)
+  {
+    const Vector3d p = C.row(P(j));
+    D.col(j) = (V.rowwise()-p.transpose()).rowwise().norm();
+    VectorXb vj;
+    bone_visible(V,F,ei,p,p,vj);
+    vis_mask.col(j) = vj;
+  }
+
+  // loop over bones
+  for(int j = 0;j<nb;j++)
+  {
+    const Vector3d s = C.row(BE(j,0));
+    const Vector3d d = C.row(BE(j,1));
+    VectorXd t,sqrD;
+    project_to_line_segment(V,s,d,t,sqrD);
+    D.col(np+j) = sqrD.array().sqrt();
+    VectorXb vj;
+    bone_visible(V,F,ei,s,d,vj);
+    vis_mask.col(np+j) = vj;
+  }
+
+  if(CE.rows() > 0)
+  {
+    cerr<<"Error: Cage edges are not supported. Ignored."<<endl;
+  }
+
+  MatrixXd PP = MatrixXd::Zero(n,m);
+  VectorXd min_D;
+  VectorXd Hdiag = VectorXd::Zero(n);
+  VectorXi J;
+  mat_min(D,2,min_D,J);
+  for(int i = 0;i<n;i++)
+  {
+    PP(i,J(i)) = 1;
+    if(vis_mask(i,J(i)))
+    {
+      double hii = pow(min_D(i),-2.); 
+      Hdiag(i) = (hii>1e10?1e10:hii);
+    }
+  }
+  SparseMatrix<double> Q,L,M;
+  cotmatrix(V,F,L);
+  massmatrix(V,F,MASSMATRIX_TYPE_DEFAULT,M);
+  const auto & H = Hdiag.asDiagonal();
+  Q = (-L+M*H);
+  SimplicialLLT <SparseMatrix<double > > llt;
+  llt.compute(Q);
+  switch(llt.info())
+  {
+    case Eigen::Success:
+      break;
+    case Eigen::NumericalIssue:
+      cerr<<"Error: Numerical issue."<<endl;
+      return false;
+    default:
+      cerr<<"Error: Other."<<endl;
+      return false;
+  }
+
+  const auto & rhs = M*H*PP;
+  W = llt.solve(rhs);
+  return true;
+}

include/igl/embree/bone_heat.h

@@ -0,0 +1,44 @@
+// 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_BONE_HEAT_H
+#define IGL_BONE_HEAT_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // BONE_HEAT  Compute skinning weights W given a surface mesh (V,F) and an
+  // internal skeleton (C,BE) according to "Automatic Rigging" [Baran and
+  // Popovic 2007].
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of mesh corner indices into V
+  //   C  #C by 3 list of joint locations
+  //   P  #P list of point handle indices into C
+  //   BE  #BE by 2 list of bone edge indices into C
+  //   CE  #CE by 2 list of cage edge indices into **P**
+  // Outputs:
+  //   W  #V by #P+#BE matrix of weights.
+  // Returns true only on success.
+  //
+  IGL_INLINE bool bone_heat(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & C,
+    const Eigen::VectorXi & P,
+    const Eigen::MatrixXi & BE,
+    const Eigen::MatrixXi & CE,
+    Eigen::MatrixXd & W);
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bone_heat.cpp"
+#endif
+
+#endif

include/igl/embree/bone_visible.cpp

@@ -6,7 +6,6 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "bone_visible.h"
-#include "EmbreeIntersector.h"
 #include <igl/project_to_line.h>
 #include <igl/EPS.h>
 #include <igl/Timer.h>
@@ -24,10 +23,6 @@ IGL_INLINE void igl::bone_visible(
   const Eigen::PlainObjectBase<DerivedSD> & d,
   Eigen::PlainObjectBase<Derivedflag>  & flag)
 {
-  using namespace igl;
-  using namespace std;
-  using namespace Eigen;
-  flag.resize(V.rows());
   // "double sided lighting"
   Eigen::PlainObjectBase<DerivedF> FF;
   FF.resize(F.rows()*2,F.cols());
@@ -35,13 +30,33 @@ IGL_INLINE void igl::bone_visible(
   // Initialize intersector
   EmbreeIntersector ei;
   ei.init(V.template cast<float>(),FF.template cast<int>());
+  return bone_visible(V,F,ei,s,d,flag);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedSD,
+  typename Derivedflag>
+IGL_INLINE void igl::bone_visible(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const igl::EmbreeIntersector & ei,
+  const Eigen::PlainObjectBase<DerivedSD> & s,
+  const Eigen::PlainObjectBase<DerivedSD> & d,
+  Eigen::PlainObjectBase<Derivedflag>  & flag)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  flag.resize(V.rows());
   const double sd_norm = (s-d).norm();
   // Embree seems to be parallel when constructing but not when tracing rays
 #pragma omp parallel for
   // loop over mesh vertices
   for(int v = 0;v<V.rows();v++)
   {
-    Vector3d Vv = V.row(v);
+    const Vector3d Vv = V.row(v);
     // Project vertex v onto line segment sd
     //embree.intersectSegment
     double t,sqrd;

include/igl/embree/bone_visible.h

@@ -9,29 +9,30 @@
 #define IGL_BONE_VISIBLE_H
 #include <igl/igl_inline.h>
 #include <Eigen/Core>
-//
-// BONE_VISIBLE  test whether vertices of mesh are "visible" to a given bone,
-// where "visible" is defined as in [Baran & Popovic 07]. Instead of checking
-// whether each point can see *any* of the bone, we just check if each point
-// can see its own projection onto the bone segment. In other words, we project
-// each vertex v onto the bone, projv. Then we check if there are any
-// intersections between the line segment (projv-->v) and the mesh.
-//
-// [flag] = bone_visible(V,F,s,d);
-//
-// Input:
-//    s  row vector of position of start end point of bone
-//    d  row vector of position of dest end point of bone
-//    V  #V by 3 list of vertex positions
-//    F  #F by 3 list of triangle indices
-// Output:
-//    flag  #V by 1 list of bools (true) visible, (false) obstructed
-//
-// Note: This checks for hits along the segment which are facing in *any*
-// direction from the ray.
-//
+#include "EmbreeIntersector.h"
 namespace igl
 {
+  //
+  // BONE_VISIBLE  test whether vertices of mesh are "visible" to a given bone,
+  // where "visible" is defined as in [Baran & Popovic 07]. Instead of checking
+  // whether each point can see *any* of the bone, we just check if each point
+  // can see its own projection onto the bone segment. In other words, we project
+  // each vertex v onto the bone, projv. Then we check if there are any
+  // intersections between the line segment (projv-->v) and the mesh.
+  //
+  // [flag] = bone_visible(V,F,s,d);
+  //
+  // Input:
+  //    V  #V by 3 list of vertex positions
+  //    F  #F by 3 list of triangle indices
+  //    s  row vector of position of start end point of bone
+  //    d  row vector of position of dest end point of bone
+  // Output:
+  //    flag  #V by 1 list of bools (true) visible, (false) obstructed
+  //
+  // Note: This checks for hits along the segment which are facing in *any*
+  // direction from the ray.
+  //
   template <
     typename DerivedV, 
     typename DerivedF, 
@@ -43,6 +44,20 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedSD> & s,
     const Eigen::PlainObjectBase<DerivedSD> & d,
     Eigen::PlainObjectBase<Derivedflag>  & flag);
+  // Inputs:
+  //  ei  EmbreeIntersector for mesh (V,F) should be double sided
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedSD,
+    typename Derivedflag>
+  IGL_INLINE void bone_visible(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const igl::EmbreeIntersector & ei,
+    const Eigen::PlainObjectBase<DerivedSD> & s,
+    const Eigen::PlainObjectBase<DerivedSD> & d,
+    Eigen::PlainObjectBase<Derivedflag>  & flag);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "bone_visible.cpp"

include/igl/exterior_edges.cpp

@@ -2,7 +2,7 @@
 #include "all_edges.h"
 
 #include <cassert>
-#include <map>
+#include <unordered_map>
 #include <utility>
 
 //template <typename T> inline int sgn(T val) {
@@ -40,9 +40,21 @@ IGL_INLINE void igl::exterior_edges(
   assert(F.cols() == 3);
   MatrixXi all_E;
   all_edges(F,all_E);
+  long int n = F.maxCoeff()+1;
+  int m = F.minCoeff();
+  const auto & compress = [&n,&m](const int i, const int j)->long int
+  {
+    return n*(i-m)+(j-m);
+  };
+  const auto & decompress = [&n,&m](const long int l,int & i, int & j)
+  {
+    i = (l / n) + m;
+    j = (l % n) + m;
+  };
+
   // Count occurances looking only at pairs (i,j) where i<j, so we count and
   // edge i-->j as +1 and j-->i as -1
-  map<pair<const int,const int>,int> C;
+  unordered_map<long int,int> C;
   // Loop over edges
   for(int e = 0;e<all_E.rows();e++)
   {
@@ -51,11 +63,11 @@ IGL_INLINE void igl::exterior_edges(
     if(i<j)
     {
       // Forward direction --> +1
-      C[pair<const int,const int>(i,j)]++;
+      C[compress(i,j)]++;
     }else
     {
       // Backward direction --> -1
-      C[pair<const int,const int>(j,i)]--;
+      C[compress(j,i)]--;
     }
   }
   // Q: Why mod off factors of 2? +1/-1 already takes care of interior edges?
@@ -66,25 +78,20 @@ IGL_INLINE void igl::exterior_edges(
   E.resize(all_E.rows(),all_E.cols());
   int e = 0;
   // Find all edges with -1 or 1 occurances, flipping direction for -1
-  for(
-    map<pair<const int, const int>,int>::const_iterator cit=C.begin();
-    cit!=C.end();
-    cit++)
+  for(const auto & cit : C)
   {
     int i,j;
-    if(cit->second > 0)
+    if(cit.second > 0)
     {
-      i = cit->first.first;
-      j = cit->first.second;
-    } else if(cit->second < 0)
+      decompress(cit.first,i,j);
+    } else if(cit.second < 0)
     {
-      i = cit->first.second;
-      j = cit->first.first;
-    } else if(cit->second == 0)
+      decompress(cit.first,j,i);
+    } else if(cit.second == 0)
     {
       continue;
     }
-    for(int k = 0;k<abs(cit->second);k++)
+    for(int k = 0;k<abs(cit.second);k++)
     {
       E(e,0) = i;
       E(e,1) = j;

include/igl/min_quad_with_fixed.h

@@ -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) 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

include/igl/mosek/mosek_quadprog.cpp

@@ -16,11 +16,6 @@
 #include "../harwell_boeing.h"
 #include "../EPS.h"
 
-// Little function mosek needs in order to know how to print to std out
-static void MSKAPI printstr(void *handle, char str[])
-{
-  printf("%s",str);
-}
 
 igl::MosekData::MosekData()
 {
@@ -109,6 +104,11 @@ IGL_INLINE bool igl::mosek_quadprog(
   mosek_guarded(MSK_makeenv(&env,NULL,NULL,NULL,NULL));
 #endif
   ///* Directs the log stream to the 'printstr' function. */
+  //// Little function mosek needs in order to know how to print to std out
+  //const auto & printstr = [](void *handle, char str[])
+  //{
+  //  printf("%s",str);
+  //}
   //mosek_guarded(MSK_linkfunctoenvstream(env,MSK_STREAM_LOG,NULL,printstr));
   // initialize mosek environment
   mosek_guarded(MSK_initenv(env));

include/igl/per_edge_normals.cpp

@@ -8,6 +8,7 @@
 template <
   typename DerivedV, 
   typename DerivedF, 
+  typename DerivedFN, 
   typename DerivedN,
   typename DerivedE,
   typename DerivedEMAP>
@@ -15,9 +16,11 @@ IGL_INLINE void igl::per_edge_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   const PerEdgeNormalsWeightingType weighting,
+  const Eigen::PlainObjectBase<DerivedFN>& FN,
   Eigen::PlainObjectBase<DerivedN> & N,
   Eigen::PlainObjectBase<DerivedE> & E,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+
 {
   using namespace Eigen;
   using namespace std;
@@ -32,8 +35,6 @@ IGL_INLINE void igl::per_edge_normals(
   unique_simplices(allE,E,_,EMAP);
   // now sort(allE,2) == E(EMAP,:), that is, if EMAP(i) = j, then E.row(j) is
   // the undirected edge corresponding to the directed edge allE.row(i).
-  MatrixXd FN;
-  per_face_normals(V,F,FN);
 
   Eigen::VectorXd W(F.rows());
   switch(weighting)
@@ -60,7 +61,25 @@ IGL_INLINE void igl::per_edge_normals(
     }
   }
   N.rowwise().normalize();
-  
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedN,
+  typename DerivedE,
+  typename DerivedEMAP>
+IGL_INLINE void igl::per_edge_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const PerEdgeNormalsWeightingType weighting,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  Eigen::PlainObjectBase<DerivedN> FN;
+  per_face_normals(V,F,FN);
+  return per_edge_normals(V,F,weighting,FN,N,E,EMAP);
 }
 
 template <

include/igl/per_edge_normals.h

@@ -31,6 +31,21 @@ namespace igl
   //   E  #E by 2 matrix of edge indices per row
   //   EMAP  #E by 1 matrix of indices from all edges to E
   //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedFN,
+    typename DerivedN,
+    typename DerivedE,
+    typename DerivedEMAP>
+  IGL_INLINE void per_edge_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const PerEdgeNormalsWeightingType weight,
+    const Eigen::PlainObjectBase<DerivedFN>& FN,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
   template <
     typename DerivedV, 
     typename DerivedF, 
@@ -63,4 +78,3 @@ namespace igl
 #endif
 
 #endif
-

include/igl/project_to_line_segment.cpp

@@ -28,13 +28,14 @@ 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);
     if(t(p)<0)
     {
-      sqrD(p) = (P.row(p)-S).squaredNorm();
+      sqrD(p) = (Pp-S).squaredNorm();
       t(p) = 0;
     }else if(t(p)>1)
     {
-      sqrD(p) = (P.row(p)-D).squaredNorm();
+      sqrD(p) = (Pp-D).squaredNorm();
       t(p) = 1;
     }
   }

include/igl/readOBJ.cpp

@@ -203,6 +203,17 @@ IGL_INLINE bool igl::readOBJ(
   return true;
 }
 
+template <typename Scalar, typename Index>
+IGL_INLINE bool igl::readOBJ(
+  const std::string obj_file_name, 
+  std::vector<std::vector<Scalar > > & V,
+  std::vector<std::vector<Index > > & F)
+{
+  std::vector<std::vector<Scalar > > TC,N;
+  std::vector<std::vector<Index > > FTC,FN;
+  return readOBJ(obj_file_name,V,TC,N,F,FTC,FN);
+}
+
 template <typename DerivedV, typename DerivedF, typename DerivedT>
 IGL_INLINE bool igl::readOBJ(
   const std::string str,

include/igl/readOBJ.h

@@ -47,6 +47,12 @@ namespace igl
     std::vector<std::vector<Index > > & F,
     std::vector<std::vector<Index > > & FTC,
     std::vector<std::vector<Index > > & FN);
+  // Just V and F
+  template <typename Scalar, typename Index>
+  IGL_INLINE bool readOBJ(
+    const std::string obj_file_name, 
+    std::vector<std::vector<Scalar > > & V,
+    std::vector<std::vector<Index > > & F);
 
 #ifndef IGL_NO_EIGEN
   //! Read a mesh from an ascii obj file

include/igl/slice_into.cpp

@@ -20,11 +20,11 @@ IGL_INLINE void igl::slice_into(
   Eigen::SparseMatrix<T>& Y)
 {
 
+#ifndef NDEBUG
   int xm = X.rows();
   int xn = X.cols();
   assert(R.size() == xm);
   assert(C.size() == xn);
-#ifndef NDEBUG
   int ym = Y.size();
   int yn = Y.size();
   assert(R.minCoeff() >= 0);
@@ -57,9 +57,9 @@ IGL_INLINE void igl::slice_into(
 
   int xm = X.rows();
   int xn = X.cols();
+#ifndef NDEBUG
   assert(R.size() == xm);
   assert(C.size() == xn);
-#ifndef NDEBUG
   int ym = Y.size();
   int yn = Y.size();
   assert(R.minCoeff() >= 0);

include/igl/tetgen/cdt.cpp

@@ -0,0 +1,66 @@
+// 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 "cdt.h"
+#include "../bounding_box.h"
+#include "../writeOBJ.h"
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedTV, 
+  typename DerivedTT, 
+  typename DerivedTF>
+IGL_INLINE bool igl::cdt(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const igl::CDTParam & param,
+  Eigen::PlainObjectBase<DerivedTV>& TV,
+  Eigen::PlainObjectBase<DerivedTT>& TT,
+  Eigen::PlainObjectBase<DerivedTF>& TF)
+{
+  using namespace Eigen;
+  using namespace std;
+  typedef Eigen::PlainObjectBase<DerivedV> MatrixXS;
+  typedef Eigen::PlainObjectBase<DerivedF> MatrixXI;
+  // Effective input mesh
+  MatrixXS U;
+  MatrixXI G;
+  if(param.use_bounding_box)
+  {
+    // Construct bounding box mesh
+    MatrixXS BV;
+    MatrixXI BF;
+    bounding_box(V,BV,BF);
+    // scale bounding box
+    const RowVector3d mid = 
+     (BV.colwise().minCoeff() + BV.colwise().maxCoeff()).eval()*0.5;
+    BV.rowwise() -= mid;
+    assert(param.bounding_box_scale >= 1.);
+    BV.array() *= param.bounding_box_scale;
+    BV.rowwise() += mid;
+    // Append bounding box to mesh
+    U.resize(V.rows()+BV.rows(),V.cols());
+    U<<V,BV;
+    BF.array() += V.rows();
+    G.resize(F.rows()+BF.rows(),F.cols());
+    G<<F,BF;
+  }else
+  {
+    // needless copies
+    U = V;
+    G = F;
+  }
+  // effective flags;
+  string flags = param.flags + (param.use_bounding_box ? "" : "c");
+  writeOBJ("UG.obj",U,G);
+  return tetrahedralize(U,G,flags,TV,TT,TF);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/tetgen/cdt.h

@@ -0,0 +1,66 @@
+// 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_CDT_H
+#define IGL_CDT_H
+#include "../igl_inline.h"
+
+#include <Eigen/Core>
+#include <string>
+#ifndef TETLIBRARY
+#define TETLIBRARY 
+#endif
+#include "tetgen.h" // Defined REAL
+
+namespace igl
+{
+  struct CDTParam
+  {
+    // Tetgen can compute mesh of convex hull of input (i.e. "c") but often
+    // chokes. One workaround is to force it to mesh the entire bounding box.
+    // {false}
+    bool use_bounding_box = false;
+    // Scale the bounding box a bit so that vertices near it do not give tetgen
+    // problems. {1.01}
+    double bounding_box_scale = 1.01;
+    // Flags to tetgen. Do not include the "c" flag here! {"Y"}
+    std::string flags = "Y";
+  };
+  // Create a constrained delaunay tesselation containing convex hull of the
+  // given **non-selfintersecting** mesh.
+  //
+  // Inputs:
+  //    V  #V by 3 list of input mesh vertices
+  //    F  #F by 3 list of input mesh facets
+  //    param  see above
+  //    TV  #TV by 3 list of output mesh vertices (V come first)
+  //    TT  #TT by 3 list of tetrahedra indices into TV.
+  //    TF  #TF by 3 list of facets from F potentially subdivided.
+  // 
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedTV, 
+    typename DerivedTT, 
+    typename DerivedTF>
+  IGL_INLINE bool cdt(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const igl::CDTParam & param,
+    Eigen::PlainObjectBase<DerivedTV>& TV,
+    Eigen::PlainObjectBase<DerivedTT>& TT,
+    Eigen::PlainObjectBase<DerivedTF>& TF);
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cdt.cpp"
+#endif
+
+#endif
+
+

include/igl/tetgen/mesh_to_tetgenio.cpp

@@ -35,11 +35,12 @@ IGL_INLINE bool igl::mesh_to_tetgenio(
 
   in.numberoffacets = F.size();
   in.facetlist = new tetgenio::facet[in.numberoffacets];
-  in.facetmarkerlist = NULL;
+  in.facetmarkerlist = new int[in.numberoffacets];
 
   // loop over face
   for(int i = 0;i < (int)F.size(); i++)
   {
+    in.facetmarkerlist[i] = i;
     tetgenio::facet * f = &in.facetlist[i];
     f->numberofpolygons = 1;
     f->polygonlist = new tetgenio::polygon[f->numberofpolygons];

include/igl/tetgen/tetgenio_to_tetmesh.cpp

@@ -16,7 +16,8 @@
 IGL_INLINE bool igl::tetgenio_to_tetmesh(
   const tetgenio & out,
   std::vector<std::vector<REAL > > & V, 
-  std::vector<std::vector<int> > & T)
+  std::vector<std::vector<int> > & T,
+  std::vector<std::vector<int> > & F)
 {
   using namespace std;
   // process points
@@ -62,14 +63,42 @@ IGL_INLINE bool igl::tetgenio_to_tetmesh(
   assert(max_index >= 0);
   assert(max_index < (int)V.size());
 
+  cout<<out.numberoftrifaces<<endl;
+
+  // When would this not be 4?
+  F.clear();
+  // loop over tetrahedra
+  for(int i = 0; i < out.numberoftrifaces; i++)
+  {
+    if(out.trifacemarkerlist[i]>=0)
+    {
+      vector<int> face(3);
+      for(int j = 0; j<3; j++)
+      {
+        face[j] = out.trifacelist[i * 3 + j];
+      }
+      F.push_back(face);
+    }
+  }
+
   return true;
 }
 
-template <typename DerivedV, typename DerivedT>
+IGL_INLINE bool igl::tetgenio_to_tetmesh(
+  const tetgenio & out,
+  std::vector<std::vector<REAL > > & V, 
+  std::vector<std::vector<int> > & T)
+{
+  std::vector<std::vector<int> > F;
+  return igl::tetgenio_to_tetmesh(out,V,T,F);
+}
+
+template <typename DerivedV, typename DerivedT, typename DerivedF>
 IGL_INLINE bool igl::tetgenio_to_tetmesh(
   const tetgenio & out,
   Eigen::PlainObjectBase<DerivedV>& V,
-  Eigen::PlainObjectBase<DerivedT>& T)
+  Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<DerivedF>& F)
 {
   using namespace igl;
   using namespace std;
@@ -95,6 +124,16 @@ IGL_INLINE bool igl::tetgenio_to_tetmesh(
   return true;
 }
 
+template <typename DerivedV, typename DerivedT>
+IGL_INLINE bool igl::tetgenio_to_tetmesh(
+  const tetgenio & out,
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedT>& T)
+{
+  Eigen::PlainObjectBase<DerivedT> F;
+  return tetgenio_to_tetmesh(out,V,T,F);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template bool igl::tetgenio_to_tetmesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(tetgenio const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/tetgen/tetgenio_to_tetmesh.h

@@ -21,7 +21,13 @@ namespace igl
   // Outputs:
   //   V  #V by 3 vertex position list
   //   T  #T by 4 list of tetrahedra indices into V
+  //   F  #F by 3 list of marked facets
   // Returns true on success, false on error
+  IGL_INLINE bool tetgenio_to_tetmesh(
+    const tetgenio & out,
+    std::vector<std::vector<REAL > > & V, 
+    std::vector<std::vector<int> > & T,
+    std::vector<std::vector<int> > & F);
   IGL_INLINE bool tetgenio_to_tetmesh(
     const tetgenio & out,
     std::vector<std::vector<REAL > > & V, 
@@ -31,6 +37,12 @@ namespace igl
   // Templates:
   //   DerivedV  real-value: i.e. from MatrixXd
   //   DerivedT  integer-value: i.e. from MatrixXi
+  template <typename DerivedV, typename DerivedT, typename DerivedF>
+  IGL_INLINE bool tetgenio_to_tetmesh(
+    const tetgenio & out,
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<DerivedF>& F);
   template <typename DerivedV, typename DerivedT>
   IGL_INLINE bool tetgenio_to_tetmesh(
     const tetgenio & out,

include/igl/tetgen/tetrahedralize.cpp

@@ -50,12 +50,12 @@ IGL_INLINE int igl::tetrahedralize(
     cerr<<"^"<<__FUNCTION__<<": Tetgen failed to create tets"<<endl;
     return 2;
   }
-  success = tetgenio_to_tetmesh(out,TV,TT);
+  success = tetgenio_to_tetmesh(out,TV,TT,TF);
   if(!success)
   {
     return -1;
   }
-  boundary_facets(TT,TF);
+  //boundary_facets(TT,TF);
   return 0;
 }