Merge remote-tracking branch 'refs/remotes/origin/master'

Former-commit-id: 003bf677dd7d6c5c1787273014ef02f92c56dc9e

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/copyleft/boolean/CSGTree.h

@@ -26,13 +26,14 @@ namespace igl
       //template <typename DerivedF>
       class CSGTree
       {
-        private:
+        public:
           typedef CGAL::Epeck::FT ExactScalar;
-          typedef Eigen::Matrix<ExactScalar,Eigen::Dynamic,3> MatrixX3E;
           //typedef Eigen::PlainObjectBase<DerivedF> POBF;
           typedef Eigen::MatrixXi POBF;
           typedef POBF::Index FIndex;
+          typedef Eigen::Matrix<ExactScalar,Eigen::Dynamic,3> MatrixX3E;
           typedef Eigen::Matrix<FIndex,Eigen::Dynamic,1> VectorJ;
+        private:
           // Resulting mesh
           MatrixX3E m_V;
           POBF m_F;

include/igl/copyleft/boolean/mesh_boolean.cpp

@@ -65,6 +65,23 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
       DerivedFC FF(FA.rows() + FB.rows(), 3);
       VV << VA, VB;
       FF << FA, FB.array() + VA.rows();
+      //// Handle annoying empty cases
+      //if(VA.size()>0)
+      //{
+      //  VV<<VA;
+      //}
+      //if(VB.size()>0)
+      //{
+      //  VV<<VB;
+      //}
+      //if(FA.size()>0)
+      //{
+      //  FF<<FA;
+      //}
+      //if(FB.size()>0)
+      //{
+      //  FF<<FB.array()+VA.rows();
+      //}
       resolve_fun(VV, FF, V, F, CJ);
   }
 
@@ -165,13 +182,13 @@ IGL_INLINE void igl::copyleft::boolean::mesh_boolean(
     Eigen::PlainObjectBase<DerivedFC > & FC,
     Eigen::PlainObjectBase<DerivedJ > & J)
 {
-  typedef CGAL::Epeck Kernel;
-  typedef Kernel::FT ExactScalar;
-  typedef Eigen::Matrix<
-    ExactScalar,
-    Eigen::Dynamic,
-    Eigen::Dynamic,
-    DerivedVC::IsRowMajor> MatrixXES;
+  //typedef CGAL::Epeck Kernel;
+  //typedef Kernel::FT ExactScalar;
+  //typedef Eigen::Matrix<
+  //  ExactScalar,
+  //  Eigen::Dynamic,
+  //  Eigen::Dynamic,
+  //  DerivedVC::IsRowMajor> MatrixXES;
 
   switch (type) {
     case MESH_BOOLEAN_TYPE_UNION:

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -354,7 +354,10 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
     tit != T.end(); 
     ++tit)
   {
-    boxes.push_back(Box(tit->bbox(), tit));
+    if (!tit->is_degenerate())
+    {
+      boxes.push_back(Box(tit->bbox(), tit));
+    }
   }
   // Leapfrog callback
   std::function<void(const Box &a,const Box &b)> cb = 
@@ -831,10 +834,10 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
   Index fb = b.handle()-T.begin();
   const Triangle_3 & A = *a.handle();
   const Triangle_3 & B = *b.handle();
-  // I'm not going to deal with degenerate triangles, though at some point we
-  // should
-  assert(!a.handle()->is_degenerate());
-  assert(!b.handle()->is_degenerate());
+  //// I'm not going to deal with degenerate triangles, though at some point we
+  //// should
+  //assert(!a.handle()->is_degenerate());
+  //assert(!b.handle()->is_degenerate());
   // Number of combinatorially shared vertices
   Index comb_shared_vertices = 0;
   // Number of geometrically shared vertices (*not* including combinatorially

include/igl/copyleft/cgal/closest_facet.cpp

@@ -10,6 +10,7 @@
 
 #include <vector>
 #include <stdexcept>
+#include <unordered_map>
 
 #include <CGAL/AABB_tree.h>
 #include <CGAL/AABB_traits.h>
@@ -18,343 +19,412 @@
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
 
 #include "order_facets_around_edge.h"
+#include "../../vertex_triangle_adjacency.h"
+#include "../../writePLY.h"
 
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedI,
-    typename DerivedP,
-    typename DerivedR,
-    typename DerivedS >
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename DerivedP,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename DerivedR,
+  typename DerivedS >
 IGL_INLINE void igl::copyleft::cgal::closest_facet(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedI>& I,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        Eigen::PlainObjectBase<DerivedR>& R,
-        Eigen::PlainObjectBase<DerivedS>& S) {
-    typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
-    typedef Kernel::Point_3 Point_3;
-    typedef Kernel::Plane_3 Plane_3;
-    typedef Kernel::Segment_3 Segment_3;
-    typedef Kernel::Triangle_3 Triangle;
-    typedef std::vector<Triangle>::iterator Iterator;
-    typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
-    typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
-    typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
-
-    if (F.rows() <= 0 || I.rows() <= 0) {
-        throw std::runtime_error(
-                "Closest facet cannot be computed on empty mesh.");
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedI>& I,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedS>& S) {
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+  typedef Kernel::Point_3 Point_3;
+  typedef Kernel::Plane_3 Plane_3;
+  typedef Kernel::Segment_3 Segment_3;
+  typedef Kernel::Triangle_3 Triangle;
+  typedef std::vector<Triangle>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+
+  if (F.rows() <= 0 || I.rows() <= 0) {
+    throw std::runtime_error(
+        "Closest facet cannot be computed on empty mesh.");
+  }
+
+  std::vector<std::vector<size_t> > VF;
+  std::vector<std::vector<size_t> > VFi;
+  igl::vertex_triangle_adjacency(V.rows(), F, VF, VFi);
+
+  std::vector<bool> in_I(F.rows(), false);
+  const size_t num_faces = I.rows();
+  std::vector<Triangle> triangles;
+  for (size_t i=0; i<num_faces; i++) {
+    const Eigen::Vector3i f = F.row(I(i, 0));
+    in_I[I(i,0)] = true;
+    triangles.emplace_back(
+        Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
+        Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
+        Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
+    if (triangles.back().is_degenerate()) {
+      throw std::runtime_error(
+          "Input facet components contains degenerated triangles");
+    }
+  }
+  Tree tree(triangles.begin(), triangles.end());
+  tree.accelerate_distance_queries();
+
+  auto on_the_positive_side = [&](size_t fid, const Point_3& p) {
+    const auto& f = F.row(fid).eval();
+    Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+    Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+    Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+    auto ori = CGAL::orientation(v0, v1, v2, p);
+    switch (ori) {
+      case CGAL::POSITIVE:
+        return true;
+      case CGAL::NEGATIVE:
+        return false;
+      case CGAL::COPLANAR:
+        // Warning:
+        // This can only happen if fid contains a boundary edge.
+        // Catergorized this ambiguous case as negative side.
+        return false;
+      default:
+        throw std::runtime_error("Unknown CGAL state.");
+    }
+    return false;
+  };
+
+  auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
+    const auto& f = F.row(fid);
+    if      ((size_t)f[0] == s && (size_t)f[1] == d) return false;
+    else if ((size_t)f[1] == s && (size_t)f[2] == d) return false;
+    else if ((size_t)f[2] == s && (size_t)f[0] == d) return false;
+    else if ((size_t)f[0] == d && (size_t)f[1] == s) return true;
+    else if ((size_t)f[1] == d && (size_t)f[2] == s) return true;
+    else if ((size_t)f[2] == d && (size_t)f[0] == s) return true;
+    else {
+      throw std::runtime_error(
+          "Cannot compute orientation due to incorrect connectivity");
+      return false;
+    }
+  };
+  auto index_to_signed_index = [&](size_t index, bool ori) -> int{
+    return (index+1) * (ori? 1:-1);
+  };
+  //auto signed_index_to_index = [&](int signed_index) -> size_t {
+  //    return abs(signed_index) - 1;
+  //};
+
+  enum ElementType { VERTEX, EDGE, FACE };
+  auto determine_element_type = [&](const Point_3& p, const size_t fid,
+      size_t& element_index) {
+    const auto& tri = triangles[fid];
+    const Point_3 p0 = tri[0];
+    const Point_3 p1 = tri[1];
+    const Point_3 p2 = tri[2];
+
+    if (p == p0) { element_index = 0; return VERTEX; }
+    if (p == p1) { element_index = 1; return VERTEX; }
+    if (p == p2) { element_index = 2; return VERTEX; }
+    if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
+    if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
+    if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
+
+    element_index = 0;
+    return FACE;
+  };
+
+  auto process_edge_case = [&](
+      size_t query_idx,
+      const size_t s, const size_t d,
+      size_t preferred_facet,
+      bool& orientation) {
+    Point_3 query_point(
+        P(query_idx, 0),
+        P(query_idx, 1),
+        P(query_idx, 2));
+
+    size_t corner_idx = std::numeric_limits<size_t>::max();
+    if ((s == F(preferred_facet, 0) && d == F(preferred_facet, 1)) ||
+        (s == F(preferred_facet, 1) && d == F(preferred_facet, 0))) {
+      corner_idx = 2;
+    } else if ((s == F(preferred_facet, 0) && d == F(preferred_facet, 2)) ||
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 0))) {
+      corner_idx = 1;
+    } else if ((s == F(preferred_facet, 1) && d == F(preferred_facet, 2)) ||
+        (s == F(preferred_facet, 2) && d == F(preferred_facet, 1))) {
+      corner_idx = 0;
+    } else {
+      std::cerr << "s: " << s << "\t d:" << d << std::endl;
+      std::cerr << F.row(preferred_facet) << std::endl;
+      throw std::runtime_error(
+          "Invalid connectivity, edge does not belong to facet");
     }
 
-    const size_t num_faces = I.rows();
-    std::vector<Triangle> triangles;
-    for (size_t i=0; i<num_faces; i++) {
-        const Eigen::Vector3i f = F.row(I(i, 0));
-        triangles.emplace_back(
-                Point_3(V(f[0], 0), V(f[0], 1), V(f[0], 2)),
-                Point_3(V(f[1], 0), V(f[1], 1), V(f[1], 2)),
-                Point_3(V(f[2], 0), V(f[2], 1), V(f[2], 2)));
-        if (triangles.back().is_degenerate()) {
-            throw std::runtime_error(
-                    "Input facet components contains degenerated triangles");
-        }
+    auto ueid = EMAP(preferred_facet + corner_idx * F.rows());
+    auto eids = uE2E[ueid];
+    std::vector<size_t> intersected_face_indices;
+    for (auto eid : eids) {
+      const size_t fid = eid % F.rows();
+      if (in_I[fid]) {
+        intersected_face_indices.push_back(fid);
+      }
     }
-    Tree tree(triangles.begin(), triangles.end());
-    tree.accelerate_distance_queries();
-
-    auto on_the_positive_side = [&](size_t fid, const Point_3& p) {
-        const auto& f = F.row(fid).eval();
-        Point_3 v0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
-        Point_3 v1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
-        Point_3 v2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
-        auto ori = CGAL::orientation(v0, v1, v2, p);
-        switch (ori) {
-            case CGAL::POSITIVE:
-                return true;
-            case CGAL::NEGATIVE:
-                return false;
-            case CGAL::COPLANAR:
-                throw std::runtime_error(
-                        "It seems input mesh contains self intersection");
-            default:
-                throw std::runtime_error("Unknown CGAL state.");
-        }
-        return false;
-    };
 
-    auto get_orientation = [&](size_t fid, size_t s, size_t d) -> bool {
-        const auto& f = F.row(fid);
-        if      ((size_t)f[0] == s && (size_t)f[1] == d) return false;
-        else if ((size_t)f[1] == s && (size_t)f[2] == d) return false;
-        else if ((size_t)f[2] == s && (size_t)f[0] == d) return false;
-        else if ((size_t)f[0] == d && (size_t)f[1] == s) return true;
-        else if ((size_t)f[1] == d && (size_t)f[2] == s) return true;
-        else if ((size_t)f[2] == d && (size_t)f[0] == s) return true;
-        else {
-            throw std::runtime_error(
-                    "Cannot compute orientation due to incorrect connectivity");
-            return false;
-        }
-    };
-    auto index_to_signed_index = [&](size_t index, bool ori) -> int{
-        return (index+1) * (ori? 1:-1);
-    };
-    //auto signed_index_to_index = [&](int signed_index) -> size_t {
-    //    return abs(signed_index) - 1;
-    //};
-
-    enum ElementType { VERTEX, EDGE, FACE };
-    auto determine_element_type = [&](const Point_3& p, const size_t fid,
-            size_t& element_index) {
-        const auto& tri = triangles[fid];
-        const Point_3 p0 = tri[0];
-        const Point_3 p1 = tri[1];
-        const Point_3 p2 = tri[2];
-
-        if (p == p0) { element_index = 0; return VERTEX; }
-        if (p == p1) { element_index = 1; return VERTEX; }
-        if (p == p2) { element_index = 2; return VERTEX; }
-        if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
-        if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
-        if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
-
-        element_index = 0;
-        return FACE;
-    };
+    const size_t num_intersected_faces = intersected_face_indices.size();
+    std::vector<int> intersected_face_signed_indices(num_intersected_faces);
+    std::transform(
+        intersected_face_indices.begin(),
+        intersected_face_indices.end(),
+        intersected_face_signed_indices.begin(),
+        [&](size_t index) {
+        return index_to_signed_index(
+          index, get_orientation(index, s,d));
+        });
+
+    assert(num_intersected_faces >= 1);
+    if (num_intersected_faces == 1) {
+      // The edge must be a boundary edge.  Thus, the orientation can be
+      // simply determined by checking if the query point is on the
+      // positive side of the facet.
+      const size_t fid = intersected_face_indices[0];
+      orientation = on_the_positive_side(fid, query_point);
+      return fid;
+    }
 
-    auto process_edge_case = [&](
-            size_t query_idx,
-            const size_t s, const size_t d,
-            size_t preferred_facet,
-            bool& orientation) {
-
-        Point_3 mid_edge_point(
-                (V(s,0) + V(d,0)) * 0.5,
-                (V(s,1) + V(d,1)) * 0.5,
-                (V(s,2) + V(d,2)) * 0.5);
-        Point_3 query_point(
-                P(query_idx, 0),
-                P(query_idx, 1),
-                P(query_idx, 2));
-
-        std::vector<Tree::Primitive_id> intersected_faces;
-        tree.all_intersected_primitives(Segment_3(mid_edge_point, query_point),
-                std::back_inserter(intersected_faces));
-
-        const size_t num_intersected_faces = intersected_faces.size();
-        std::vector<size_t> intersected_face_indices(num_intersected_faces);
-        std::vector<int> intersected_face_signed_indices(num_intersected_faces);
-        std::transform(intersected_faces.begin(),
-                intersected_faces.end(),
-                intersected_face_indices.begin(),
-                [&](const Tree::Primitive_id& itr) -> size_t
-                { return I(itr-triangles.begin(), 0); });
-        std::transform(
-                intersected_face_indices.begin(),
-                intersected_face_indices.end(),
-                intersected_face_signed_indices.begin(),
-                [&](size_t index) {
-                    return index_to_signed_index(
-                        index, get_orientation(index, s,d));
-                });
-
-        assert(num_intersected_faces >= 1);
-        if (num_intersected_faces == 1) {
-            // The edge must be a boundary edge.  Thus, the orientation can be
-            // simply determined by checking if the query point is on the
-            // positive side of the facet.
-            const size_t fid = intersected_face_indices[0];
-            orientation = on_the_positive_side(fid, query_point);
-            return fid;
+    Eigen::VectorXi order;
+    DerivedP pivot = P.row(query_idx).eval();
+    igl::copyleft::cgal::order_facets_around_edge(V, F, s, d,
+        intersected_face_signed_indices,
+        pivot, order);
+
+    // Although first and last are equivalent, make the choice based on
+    // preferred_facet.
+    const size_t first = order[0];
+    const size_t last = order[num_intersected_faces-1];
+    if (intersected_face_indices[first] == preferred_facet) {
+      orientation = intersected_face_signed_indices[first] < 0;
+      return intersected_face_indices[first];
+    } else if (intersected_face_indices[last] == preferred_facet) {
+      orientation = intersected_face_signed_indices[last] > 0;
+      return intersected_face_indices[last];
+    } else {
+      orientation = intersected_face_signed_indices[order[0]] < 0;
+      return intersected_face_indices[order[0]];
+    }
+  };
+
+  auto process_face_case = [&](
+      const size_t query_idx, const Point_3& closest_point,
+      const size_t fid, bool& orientation) {
+    const auto& f = F.row(I(fid, 0));
+    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) {
+    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.
+      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++) {
+        const size_t fid = all_adj_faces[i];
+        if (in_I[fid]) {
+          adj_faces.push_back(fid);
+          adj_face_corners.push_back(all_adj_face_corners[i]);
         }
+      }
+    }
+    const size_t num_adj_faces = adj_faces.size();
+    assert(num_adj_faces > 0);
+
+    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++) {
+      const size_t fid = adj_faces[i];
+      const size_t cid = adj_face_corners[i];
+      const auto& f = F.row(adj_faces[i]);
+      const size_t next = f[(cid+1)%3];
+      const size_t prev = f[(cid+2)%3];
+      adj_vertices_set.insert(next);
+      adj_vertices_set.insert(prev);
+      v2f.insert({{next, fid}, {prev, fid}});
+    }
+    const size_t num_adj_vertices = adj_vertices_set.size();
+    std::vector<size_t> adj_vertices(num_adj_vertices);
+    std::copy(adj_vertices_set.begin(), adj_vertices_set.end(),
+        adj_vertices.begin());
+
+    std::vector<Point_3> adj_points;
+    for (size_t vid : adj_vertices) {
+      adj_points.emplace_back(V(vid,0), V(vid,1), V(vid,2));
+    }
 
-        Eigen::VectorXi order;
-        DerivedP pivot = P.row(query_idx).eval();
-        igl::copyleft::cgal::order_facets_around_edge(V, F, s, d,
-                intersected_face_signed_indices,
-                pivot, order);
-
-        // Although first and last are equivalent, make the choice based on
-        // preferred_facet.
-        const size_t first = order[0];
-        const size_t last = order[num_intersected_faces-1];
-        if (intersected_face_indices[first] == preferred_facet) {
-            orientation = intersected_face_signed_indices[first] < 0;
-            return intersected_face_indices[first];
-        } else if (intersected_face_indices[last] == preferred_facet) {
-            orientation = intersected_face_signed_indices[last] > 0;
-            return intersected_face_indices[last];
-        } else {
-            orientation = intersected_face_signed_indices[order[0]] < 0;
-            return intersected_face_indices[order[0]];
+    // A plane is on the exterior if all adj_points lies on or to
+    // one side of the plane.
+    auto is_on_exterior = [&](const Plane_3& separator) {
+      size_t positive=0;
+      size_t negative=0;
+      size_t coplanar=0;
+      for (const auto& point : adj_points) {
+        switch(separator.oriented_side(point)) {
+          case CGAL::ON_POSITIVE_SIDE:
+            positive++;
+            break;
+          case CGAL::ON_NEGATIVE_SIDE:
+            negative++;
+            break;
+          case CGAL::ON_ORIENTED_BOUNDARY:
+            coplanar++;
+            break;
+          default:
+            throw "Unknown plane-point orientation";
         }
+      }
+      auto query_orientation = separator.oriented_side(query_point);
+      if (query_orientation == CGAL::ON_ORIENTED_BOUNDARY &&
+          (positive == 0 && negative == 0)) {
+        // All adj vertices and query point are coplanar.
+        // In this case, all separators are equally valid.
+        return true;
+      } else {
+        bool r = (positive == 0 && query_orientation == CGAL::POSITIVE)
+          || (negative == 0 && query_orientation == CGAL::NEGATIVE);
+        return r;
+      }
     };
 
-    auto process_face_case = [&](
-            const size_t query_idx, const size_t fid, bool& orientation) {
-        const auto& f = F.row(I(fid, 0));
-        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) {
-        Point_3 closest_point(V(s, 0), V(s, 1), V(s, 2));
-        Point_3 query_point(P(query_idx, 0), P(query_idx, 1), P(query_idx, 2));
-
-        std::vector<Tree::Primitive_id> intersected_faces;
-        tree.all_intersected_primitives(Segment_3(closest_point, query_point),
-                std::back_inserter(intersected_faces));
-
-        const size_t num_intersected_faces = intersected_faces.size();
-        std::vector<size_t> intersected_face_indices(num_intersected_faces);
-        std::transform(intersected_faces.begin(),
-                intersected_faces.end(),
-                intersected_face_indices.begin(),
-                [&](const Tree::Primitive_id& itr) -> size_t
-                { return I(itr-triangles.begin(), 0); });
-
-        std::set<size_t> adj_vertices_set;
-        for (auto fid : intersected_face_indices) {
-            const auto& f = F.row(fid);
-            if ((size_t)f[0] != s) adj_vertices_set.insert(f[0]);
-            if ((size_t)f[1] != s) adj_vertices_set.insert(f[1]);
-            if ((size_t)f[2] != s) adj_vertices_set.insert(f[2]);
+    size_t d = std::numeric_limits<size_t>::max();
+    for (size_t i=0; i<num_adj_vertices; i++) {
+      const size_t vi = adj_vertices[i];
+      for (size_t j=i+1; j<num_adj_vertices; j++) {
+        Plane_3 separator(closest_point, adj_points[i], adj_points[j]);
+        if (separator.is_degenerate()) {
+          continue;
         }
-        const size_t num_adj_vertices = adj_vertices_set.size();
-        std::vector<size_t> adj_vertices(num_adj_vertices);
-        std::copy(adj_vertices_set.begin(), adj_vertices_set.end(),
-                adj_vertices.begin());
-
-        std::vector<Point_3> adj_points;
-        for (size_t vid : adj_vertices) {
-            adj_points.emplace_back(V(vid,0), V(vid,1), V(vid,2));
+        if (is_on_exterior(separator)) {
+          if (!CGAL::collinear(
+                query_point, adj_points[i], closest_point)) {
+            d = vi;
+            break;
+          } else {
+            d = adj_vertices[j];
+            assert(!CGAL::collinear(
+                  query_point, adj_points[j], closest_point));
+            break;
+          }
         }
-
-        // A plane is on the exterior if all adj_points lies on or to
-        // one side of the plane.
-        auto is_on_exterior = [&](const Plane_3& separator) {
-            size_t positive=0;
-            size_t negative=0;
-            size_t coplanar=0;
-            for (const auto& point : adj_points) {
-                switch(separator.oriented_side(point)) {
-                    case CGAL::ON_POSITIVE_SIDE:
-                        positive++;
-                        break;
-                    case CGAL::ON_NEGATIVE_SIDE:
-                        negative++;
-                        break;
-                    case CGAL::ON_ORIENTED_BOUNDARY:
-                        coplanar++;
-                        break;
-                    default:
-                        throw "Unknown plane-point orientation";
-                }
-            }
-            auto query_orientation = separator.oriented_side(query_point);
-            bool r = (positive == 0 && query_orientation == CGAL::POSITIVE)
-                || (negative == 0 && query_orientation == CGAL::NEGATIVE);
-            return r;
-        };
-
-        size_t d = std::numeric_limits<size_t>::max();
-        for (size_t i=0; i<num_adj_vertices; i++) {
-            const size_t vi = adj_vertices[i];
-            for (size_t j=i+1; j<num_adj_vertices; j++) {
-                Plane_3 separator(closest_point, adj_points[i], adj_points[j]);
-                if (separator.is_degenerate()) {
-                    throw std::runtime_error(
-                            "Input mesh contains degenerated faces");
-                }
-                if (is_on_exterior(separator)) {
-                    d = vi;
-                    assert(!CGAL::collinear(
-                                query_point, adj_points[i], closest_point));
-                    break;
-                }
-            }
+      }
+    }
+    if (d == std::numeric_limits<size_t>::max()) {
+      Eigen::MatrixXd tmp_vertices(V.rows(), V.cols());
+      for (size_t i=0; i<V.rows(); i++) {
+        for (size_t j=0; j<V.cols(); j++) {
+          tmp_vertices(i,j) = CGAL::to_double(V(i,j));
         }
-        assert(d != std::numeric_limits<size_t>::max());
-
-        return process_edge_case(query_idx, s, d, preferred_facet, orientation);
-    };
-
-    const size_t num_queries = P.rows();
-    R.resize(num_queries, 1);
-    S.resize(num_queries, 1);
-    for (size_t i=0; i<num_queries; i++) {
-        const Point_3 query(P(i,0), P(i,1), P(i,2));
-        auto projection = tree.closest_point_and_primitive(query);
-        const Point_3 closest_point = projection.first;
-        size_t fid = projection.second - triangles.begin();
-        bool fid_ori = false;
-
-        // Gether all facets sharing the closest point.
-        std::vector<Tree::Primitive_id> intersected_faces;
-        tree.all_intersected_primitives(Segment_3(closest_point, query),
-                std::back_inserter(intersected_faces));
-        const size_t num_intersected_faces = intersected_faces.size();
-        std::vector<size_t> intersected_face_indices(num_intersected_faces);
-        std::transform(intersected_faces.begin(),
-                intersected_faces.end(),
-                intersected_face_indices.begin(),
-                [&](const Tree::Primitive_id& itr) -> size_t
-                { return I(itr-triangles.begin(), 0); });
-
-        size_t element_index;
-        auto element_type = determine_element_type(closest_point, fid,
-                element_index);
-        switch(element_type) {
-            case VERTEX:
-                {
-                    const auto& f = F.row(I(fid, 0));
-                    const size_t s = f[element_index];
-                    fid = process_vertex_case(i, s, I(fid, 0), fid_ori);
-                }
-                break;
-            case EDGE:
-                {
-                    const auto& f = F.row(I(fid, 0));
-                    const size_t s = f[(element_index+1)%3];
-                    const size_t d = f[(element_index+2)%3];
-                    fid = process_edge_case(i, s, d, I(fid, 0), fid_ori);
-                }
-                break;
-            case FACE:
-                {
-                    fid = process_face_case(i, fid, fid_ori);
-                }
-                break;
-            default:
-                throw std::runtime_error("Unknown element type.");
+      }
+      Eigen::MatrixXi tmp_faces(adj_faces.size(), 3);
+      for (size_t i=0; i<adj_faces.size(); i++) {
+        tmp_faces.row(i) = F.row(adj_faces[i]);
+      }
+      igl::writePLY("debug.ply", tmp_vertices, tmp_faces, false);
+      throw std::runtime_error("Invalid vertex neighborhood");
+    }
+    const auto itr = v2f.equal_range(d);
+    assert(itr.first != itr.second);
+
+    return process_edge_case(query_idx, s, d, itr.first->second, orientation);
+  };
+
+  const size_t num_queries = P.rows();
+  R.resize(num_queries, 1);
+  S.resize(num_queries, 1);
+  for (size_t i=0; i<num_queries; i++) {
+    const Point_3 query(P(i,0), P(i,1), P(i,2));
+    auto projection = tree.closest_point_and_primitive(query);
+    const Point_3 closest_point = projection.first;
+    size_t fid = projection.second - triangles.begin();
+    bool fid_ori = false;
+
+    // Gether all facets sharing the closest point.
+    std::vector<Tree::Primitive_id> intersected_faces;
+    tree.all_intersected_primitives(Segment_3(closest_point, query),
+        std::back_inserter(intersected_faces));
+    const size_t num_intersected_faces = intersected_faces.size();
+    std::vector<size_t> intersected_face_indices(num_intersected_faces);
+    std::transform(intersected_faces.begin(),
+        intersected_faces.end(),
+        intersected_face_indices.begin(),
+        [&](const Tree::Primitive_id& itr) -> size_t
+        { return I(itr-triangles.begin(), 0); });
+
+    size_t element_index;
+    auto element_type = determine_element_type(closest_point, fid,
+        element_index);
+    switch(element_type) {
+      case VERTEX:
+        {
+          const auto& f = F.row(I(fid, 0));
+          const size_t s = f[element_index];
+          fid = process_vertex_case(i, s, I(fid, 0), fid_ori);
+        }
+        break;
+      case EDGE:
+        {
+          const auto& f = F.row(I(fid, 0));
+          const size_t s = f[(element_index+1)%3];
+          const size_t d = f[(element_index+2)%3];
+          fid = process_edge_case(i, s, d, I(fid, 0), fid_ori);
         }
+        break;
+      case FACE:
+        {
+          fid = process_face_case(i, closest_point, fid, fid_ori);
+        }
+        break;
+      default:
+        throw std::runtime_error("Unknown element type.");
+    }
 
 
-        R(i,0) = fid;
-        S(i,0) = fid_ori;
-    }
+    R(i,0) = fid;
+    S(i,0) = fid_ori;
+  }
 }
 
 template<
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedP,
-    typename DerivedR,
-    typename DerivedS >
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename uE2EType,
+  typename DerivedEMAP,
+  typename DerivedR,
+  typename DerivedS >
 IGL_INLINE void igl::copyleft::cgal::closest_facet(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedP>& P,
-        Eigen::PlainObjectBase<DerivedR>& R,
-        Eigen::PlainObjectBase<DerivedS>& S) {
-    const size_t num_faces = F.rows();
-    Eigen::VectorXi I(num_faces);
-    I.setLinSpaced(num_faces, 0, num_faces-1);
-    igl::copyleft::cgal::closest_facet(V, F, I, P, R, S);
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    const std::vector<std::vector<uE2EType> >& uE2E,
+    const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+    Eigen::PlainObjectBase<DerivedR>& R,
+    Eigen::PlainObjectBase<DerivedS>& S) {
+  const size_t num_faces = F.rows();
+  Eigen::VectorXi I(num_faces);
+  I.setLinSpaced(num_faces, 0, num_faces-1);
+  igl::copyleft::cgal::closest_facet(V, F, I, P, uE2E, EMAP, R, S);
 }
 
 #ifdef IGL_STATIC_LIBRARY
-template void igl::copyleft::cgal::closest_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::closest_facet<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/copyleft/cgal/closest_facet.h

@@ -11,6 +11,7 @@
 
 #include "../../igl_inline.h"
 #include <Eigen/Core>
+#include <vector>
 
 namespace igl
 {
@@ -35,6 +36,8 @@ namespace igl
           typename DerivedF,
           typename DerivedI,
           typename DerivedP,
+          typename uE2EType,
+          typename DerivedEMAP,
           typename DerivedR,
           typename DerivedS >
       IGL_INLINE void closest_facet(
@@ -42,6 +45,8 @@ namespace igl
         const Eigen::PlainObjectBase<DerivedF>& F,
         const Eigen::PlainObjectBase<DerivedI>& I,
         const Eigen::PlainObjectBase<DerivedP>& P,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
               Eigen::PlainObjectBase<DerivedR>& R,
               Eigen::PlainObjectBase<DerivedS>& S);
 
@@ -49,12 +54,16 @@ namespace igl
           typename DerivedV,
           typename DerivedF,
           typename DerivedP,
+          typename uE2EType,
+          typename DerivedEMAP,
           typename DerivedR,
           typename DerivedS >
       IGL_INLINE void closest_facet(
               const Eigen::PlainObjectBase<DerivedV>& V,
               const Eigen::PlainObjectBase<DerivedF>& F,
               const Eigen::PlainObjectBase<DerivedP>& P,
+              const std::vector<std::vector<uE2EType> >& uE2E,
+              const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
               Eigen::PlainObjectBase<DerivedR>& R,
               Eigen::PlainObjectBase<DerivedS>& S);
     }

include/igl/copyleft/cgal/extract_cells.cpp

@@ -229,7 +229,7 @@ IGL_INLINE size_t igl::copyleft::cgal::extract_cells(
             const auto& I = Is[i];
             Eigen::VectorXi closest_facets, closest_facet_orientations;
             igl::copyleft::cgal::closest_facet(V, F, I, queries,
-                    closest_facets, closest_facet_orientations);
+                    uE2E, EMAP, closest_facets, closest_facet_orientations);
 
             for (size_t j=0; j<num_components; j++) {
                 if (i == j) continue;

include/igl/copyleft/cgal/points_inside_component.cpp

@@ -36,30 +36,6 @@ namespace igl {
             typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
             typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
 
-            enum ElementType { VERTEX, EDGE, FACE };
-            template<typename DerivedV, typename DerivedF, typename DerivedI>
-            ElementType determine_element_type(
-                    const Eigen::PlainObjectBase<DerivedV>& V,
-                    const Eigen::PlainObjectBase<DerivedF>& F,
-                    const Eigen::PlainObjectBase<DerivedI>& I,
-                    const size_t fid, const Point_3& p,
-                    size_t& element_index) {
-                const Eigen::Vector3i f = F.row(I(fid, 0));
-                const Point_3 p0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
-                const Point_3 p1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
-                const Point_3 p2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
-
-                if (p == p0) { element_index = 0; return VERTEX; }
-                if (p == p1) { element_index = 1; return VERTEX; }
-                if (p == p2) { element_index = 2; return VERTEX; }
-                if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
-                if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
-                if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
-
-                element_index = 0;
-                return FACE;
-            }
-
             template<typename DerivedF, typename DerivedI>
             void extract_adj_faces(
                     const Eigen::PlainObjectBase<DerivedF>& F,
@@ -296,6 +272,25 @@ IGL_INLINE void igl::copyleft::cgal::points_inside_component(
     Tree tree(triangles.begin(), triangles.end());
     tree.accelerate_distance_queries();
 
+    enum ElementType { VERTEX, EDGE, FACE };
+    auto determine_element_type = [&](
+            size_t fid, const Point_3& p, size_t& element_index) -> ElementType{
+        const Eigen::Vector3i f = F.row(I(fid, 0));
+        const Point_3 p0(V(f[0], 0), V(f[0], 1), V(f[0], 2));
+        const Point_3 p1(V(f[1], 0), V(f[1], 1), V(f[1], 2));
+        const Point_3 p2(V(f[2], 0), V(f[2], 1), V(f[2], 2));
+
+        if (p == p0) { element_index = 0; return VERTEX; }
+        if (p == p1) { element_index = 1; return VERTEX; }
+        if (p == p2) { element_index = 2; return VERTEX; }
+        if (CGAL::collinear(p0, p1, p)) { element_index = 2; return EDGE; }
+        if (CGAL::collinear(p1, p2, p)) { element_index = 0; return EDGE; }
+        if (CGAL::collinear(p2, p0, p)) { element_index = 1; return EDGE; }
+
+        element_index = 0;
+        return FACE;
+    };
+
     const size_t num_queries = P.rows();
     inside.resize(num_queries, 1);
     for (size_t i=0; i<num_queries; i++) {
@@ -305,8 +300,7 @@ IGL_INLINE void igl::copyleft::cgal::points_inside_component(
         size_t fid = projection.second - triangles.begin();
 
         size_t element_index;
-        switch (determine_element_type(
-                    V, F, I, fid, closest_point, element_index)) {
+        switch (determine_element_type(fid, closest_point, element_index)) {
             case VERTEX:
                 {
                     const size_t s = F(I(fid, 0), element_index);

include/igl/copyleft/cgal/points_inside_component.h

@@ -17,41 +17,49 @@ namespace igl {
   {
     namespace cgal {
 
-        // Determine if queries points P are inside of connected facet component
-        // (V, F, I), where I indicates a subset of facets that forms the
-        // component.
-        //
-        // Precondition:
-        // The input mesh must be a closed, self-intersection free,
-        // non-degenerated surface.  Queries points must be either inside or
-        // outside of the mesh (i.e. not on the surface of the mesh).
-        //
-        // Inputs:
-        //   V  #V by 3 array of vertex positions.
-        //   F  #F by 3 array of triangles.
-        //   I  #I list of triangle indices to consider.
-        //   P  #P by 3 array of query points.
-        //
-        // Outputs:
-        //   inside  #P list of booleans that is true iff the corresponding
-        //           query point is inside of the mesh.
-        template<typename DerivedV, typename DerivedF, typename DerivedI,
-            typename DerivedP, typename DerivedB>
-            IGL_INLINE void points_inside_component(
-                    const Eigen::PlainObjectBase<DerivedV>& V,
-                    const Eigen::PlainObjectBase<DerivedF>& F,
-                    const Eigen::PlainObjectBase<DerivedI>& I,
-                    const Eigen::PlainObjectBase<DerivedP>& P,
-                    Eigen::PlainObjectBase<DerivedB>& inside);
+      // Determine if queries points P are inside of connected facet component
+      // (V, F, I), where I indicates a subset of facets that forms the
+      // component.
+      //
+      // Precondition:
+      // The input mesh must be a closed, self-intersection free,
+      // non-degenerated surface.  Queries points must be either inside or
+      // outside of the mesh (i.e. not on the surface of the mesh).
+      //
+      // Inputs:
+      //   V  #V by 3 array of vertex positions.
+      //   F  #F by 3 array of triangles.
+      //   I  #I list of triangle indices to consider.
+      //   P  #P by 3 array of query points.
+      //
+      // Outputs:
+      //   inside  #P list of booleans that is true iff the corresponding
+      //           query point is inside of the mesh.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedI,
+        typename DerivedP,
+        typename DerivedB>
+      IGL_INLINE void points_inside_component(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedI>& I,
+        const Eigen::PlainObjectBase<DerivedP>& P,
+        Eigen::PlainObjectBase<DerivedB>& inside);
 
-        // Determine if query points P is inside of the mesh (V, F).
-        // See above for precondition and I/O specs.
-        template<typename DerivedV, typename DerivedF, typename DerivedP, typename DerivedB>
-            IGL_INLINE void points_inside_component(
-                    const Eigen::PlainObjectBase<DerivedV>& V,
-                    const Eigen::PlainObjectBase<DerivedF>& F,
-                    const Eigen::PlainObjectBase<DerivedP>& P,
-                    Eigen::PlainObjectBase<DerivedB>& inside);
+      // Determine if query points P is inside of the mesh (V, F).
+      // See above for precondition and I/O specs.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedP,
+        typename DerivedB>
+      IGL_INLINE void points_inside_component(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          const Eigen::PlainObjectBase<DerivedP>& P,
+          Eigen::PlainObjectBase<DerivedB>& inside);
     }
   }
 }

include/igl/copyleft/cgal/propagate_winding_numbers.cpp

@@ -26,235 +26,234 @@
 #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;
+  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;
+      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 DerivedV,
-typename DerivedF,
-typename DerivedL,
-typename DerivedW>
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedL,
+  typename DerivedW>
 IGL_INLINE void igl::copyleft::cgal::propagate_winding_numbers(
-        const Eigen::PlainObjectBase<DerivedV>& V,
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedL>& labels,
-        Eigen::PlainObjectBase<DerivedW>& W) {
-    const size_t num_faces = F.rows();
-    //typedef typename DerivedF::Scalar Index;
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const Eigen::PlainObjectBase<DerivedL>& labels,
+    Eigen::PlainObjectBase<DerivedW>& W) {
+  const size_t num_faces = F.rows();
+  //typedef typename DerivedF::Scalar Index;
+
+  Eigen::MatrixXi E, uE;
+  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;
+  }
 
-    Eigen::MatrixXi E, uE;
-    Eigen::VectorXi EMAP;
-    std::vector<std::vector<size_t> > uE2E;
-    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
-    assert(propagate_winding_numbers_helper::is_orientable(F, uE, uE2E));
+  Eigen::VectorXi P;
+  const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
 
-    Eigen::VectorXi P;
-    const size_t num_patches = igl::extract_manifold_patches(F, EMAP, uE2E, P);
+  DerivedW per_patch_cells;
+  const size_t num_cells =
+    igl::copyleft::cgal::extract_cells(
+        V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
 
-    DerivedW per_patch_cells;
-    const size_t num_cells =
-        igl::copyleft::cgal::extract_cells(V, F, P, E, uE, uE2E, EMAP, per_patch_cells);
+  typedef std::tuple<size_t, bool, size_t> CellConnection;
+  std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
+  for (size_t i=0; i<num_patches; i++) {
+    const int positive_cell = per_patch_cells(i,0);
+    const int negative_cell = per_patch_cells(i,1);
+    cell_adjacency[positive_cell].emplace(negative_cell, false, i);
+    cell_adjacency[negative_cell].emplace(positive_cell, true, i);
+  }
 
-    typedef std::tuple<size_t, bool, size_t> CellConnection;
-    std::vector<std::set<CellConnection> > cell_adjacency(num_cells);
+  auto save_cell = [&](const std::string& filename, size_t cell_id) {
+    std::vector<size_t> faces;
     for (size_t i=0; i<num_patches; i++) {
-        const int positive_cell = per_patch_cells(i,0);
-        const int negative_cell = per_patch_cells(i,1);
-        cell_adjacency[positive_cell].emplace(negative_cell, false, i);
-        cell_adjacency[negative_cell].emplace(positive_cell, true, i);
+      if ((per_patch_cells.row(i).array() == cell_id).any()) {
+        for (size_t j=0; j<num_faces; j++) {
+          if ((size_t)P[j] == i) {
+            faces.push_back(j);
+          }
+        }
+      }
+    }
+    Eigen::MatrixXi cell_faces(faces.size(), 3);
+    for (size_t i=0; i<faces.size(); i++) {
+      cell_faces.row(i) = F.row(faces[i]);
     }
+    Eigen::MatrixXd vertices(V.rows(), 3);
+    for (size_t i=0; i<(size_t)V.rows(); i++) {
+      assign_scalar(V(i,0), vertices(i,0));
+      assign_scalar(V(i,1), vertices(i,1));
+      assign_scalar(V(i,2), vertices(i,2));
+    }
+    writePLY(filename, vertices, cell_faces);
+  };
 
-    auto save_cell = [&](const std::string& filename, size_t cell_id) {
-        std::vector<size_t> faces;
-        for (size_t i=0; i<num_patches; i++) {
-            if ((per_patch_cells.row(i).array() == cell_id).any()) {
-                for (size_t j=0; j<num_faces; j++) {
-                    if ((size_t)P[j] == i) {
-                        faces.push_back(j);
-                    }
-                }
-            }
-        }
-        Eigen::MatrixXi cell_faces(faces.size(), 3);
-        for (size_t i=0; i<faces.size(); i++) {
-            cell_faces.row(i) = F.row(faces[i]);
-        }
-        Eigen::MatrixXd vertices(V.rows(), 3);
-        for (size_t i=0; i<(size_t)V.rows(); i++) {
-            assign_scalar(V(i,0), vertices(i,0));
-            assign_scalar(V(i,1), vertices(i,1));
-            assign_scalar(V(i,2), vertices(i,2));
-        }
-        writePLY(filename, vertices, cell_faces);
+#ifndef NDEBUG
+  {
+    // Check for odd cycle.
+    Eigen::VectorXi cell_labels(num_cells);
+    cell_labels.setZero();
+    Eigen::VectorXi parents(num_cells);
+    parents.setConstant(-1);
+    auto trace_parents = [&](size_t idx) {
+      std::list<size_t> path;
+      path.push_back(idx);
+      while ((size_t)parents[path.back()] != path.back()) {
+        path.push_back(parents[path.back()]);
+      }
+      return path;
     };
-
-    {
-        // Check for odd cycle.
-        Eigen::VectorXi cell_labels(num_cells);
-        cell_labels.setZero();
-        Eigen::VectorXi parents(num_cells);
-        parents.setConstant(-1);
-        auto trace_parents = [&](size_t idx) {
-            std::list<size_t> path;
-            path.push_back(idx);
-            while ((size_t)parents[path.back()] != path.back()) {
-                path.push_back(parents[path.back()]);
-            }
-            return path;
-        };
-        for (size_t i=0; i<num_cells; i++) {
-            if (cell_labels[i] == 0) {
-                cell_labels[i] = 1;
-                std::queue<size_t> Q;
-                Q.push(i);
-                parents[i] = i;
-                while (!Q.empty()) {
-                    size_t curr_idx = Q.front();
-                    Q.pop();
-                    int curr_label = cell_labels[curr_idx];
-                    for (const auto& neighbor : cell_adjacency[curr_idx]) {
-                        if (cell_labels[std::get<0>(neighbor)] == 0) {
-                            cell_labels[std::get<0>(neighbor)] = curr_label * -1;
-                            Q.push(std::get<0>(neighbor));
-                            parents[std::get<0>(neighbor)] = curr_idx;
-                        } else {
-                            if (cell_labels[std::get<0>(neighbor)] !=
-                                    curr_label * -1) {
-                                std::cerr << "Odd cell cycle detected!" << std::endl;
-                                auto path = trace_parents(curr_idx);
-                                path.reverse();
-                                auto path2 = trace_parents(std::get<0>(neighbor));
-                                path.insert(path.end(),
-                                        path2.begin(), path2.end());
-                                for (auto cell_id : path) {
-                                    std::cout << cell_id << " ";
-                                    std::stringstream filename;
-                                    filename << "cell_" << cell_id << ".ply";
-                                    save_cell(filename.str(), cell_id);
-                                }
-                                std::cout << std::endl;
-                            }
-                            assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
-                        }
-                    }
+    for (size_t i=0; i<num_cells; i++) {
+      if (cell_labels[i] == 0) {
+        cell_labels[i] = 1;
+        std::queue<size_t> Q;
+        Q.push(i);
+        parents[i] = i;
+        while (!Q.empty()) {
+          size_t curr_idx = Q.front();
+          Q.pop();
+          int curr_label = cell_labels[curr_idx];
+          for (const auto& neighbor : cell_adjacency[curr_idx]) {
+            if (cell_labels[std::get<0>(neighbor)] == 0) {
+              cell_labels[std::get<0>(neighbor)] = curr_label * -1;
+              Q.push(std::get<0>(neighbor));
+              parents[std::get<0>(neighbor)] = curr_idx;
+            } else {
+              if (cell_labels[std::get<0>(neighbor)] !=
+                  curr_label * -1) {
+                std::cerr << "Odd cell cycle detected!" << std::endl;
+                auto path = trace_parents(curr_idx);
+                path.reverse();
+                auto path2 = trace_parents(std::get<0>(neighbor));
+                path.insert(path.end(),
+                    path2.begin(), path2.end());
+                for (auto cell_id : path) {
+                  std::cout << cell_id << " ";
+                  std::stringstream filename;
+                  filename << "cell_" << cell_id << ".ply";
+                  save_cell(filename.str(), cell_id);
                 }
+                std::cout << std::endl;
+              }
+              assert(cell_labels[std::get<0>(neighbor)] == curr_label * -1);
             }
+          }
         }
+      }
     }
+  }
+#endif
 
-    size_t outer_facet;
-    bool flipped;
-    Eigen::VectorXi I;
-    I.setLinSpaced(num_faces, 0, num_faces-1);
-    igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
+  size_t outer_facet;
+  bool flipped;
+  Eigen::VectorXi I;
+  I.setLinSpaced(num_faces, 0, num_faces-1);
+  igl::copyleft::cgal::outer_facet(V, F, I, outer_facet, flipped);
 
-    const size_t outer_patch = P[outer_facet];
-    const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
+  const size_t outer_patch = P[outer_facet];
+  const size_t infinity_cell = per_patch_cells(outer_patch, flipped?1:0);
 
-    Eigen::VectorXi patch_labels(num_patches);
-    const int INVALID = std::numeric_limits<int>::max();
-    patch_labels.setConstant(INVALID);
-    for (size_t i=0; i<num_faces; i++) {
-        if (patch_labels[P[i]] == INVALID) {
-            patch_labels[P[i]] = labels[i];
-        } else {
-            assert(patch_labels[P[i]] == labels[i]);
-        }
+  Eigen::VectorXi patch_labels(num_patches);
+  const int INVALID = std::numeric_limits<int>::max();
+  patch_labels.setConstant(INVALID);
+  for (size_t i=0; i<num_faces; i++) {
+    if (patch_labels[P[i]] == INVALID) {
+      patch_labels[P[i]] = labels[i];
+    } else {
+      assert(patch_labels[P[i]] == labels[i]);
     }
-    assert((patch_labels.array() != INVALID).all());
-    const size_t num_labels = patch_labels.maxCoeff()+1;
+  }
+  assert((patch_labels.array() != INVALID).all());
+  const size_t num_labels = patch_labels.maxCoeff()+1;
 
-    Eigen::MatrixXi per_cell_W(num_cells, num_labels);
-    per_cell_W.setConstant(INVALID);
-    per_cell_W.row(infinity_cell).setZero();
-    std::queue<size_t> Q;
-    Q.push(infinity_cell);
-    while (!Q.empty()) {
-        size_t curr_cell = Q.front();
-        Q.pop();
-        for (const auto& neighbor : cell_adjacency[curr_cell]) {
-            size_t neighbor_cell, patch_idx;
-            bool direction;
-            std::tie(neighbor_cell, direction, patch_idx) = neighbor;
-            if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
-                per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
-                for (size_t i=0; i<num_labels; i++) {
-                    int inc = (patch_labels[patch_idx] == (int)i) ?
-                        (direction ? -1:1) :0;
-                    per_cell_W(neighbor_cell, i) =
-                        per_cell_W(curr_cell, i) + inc;
-                }
-                Q.push(neighbor_cell);
-            } else {
-                for (size_t i=0; i<num_labels; i++) {
-                    if ((int)i == patch_labels[patch_idx]) {
+  Eigen::MatrixXi per_cell_W(num_cells, num_labels);
+  per_cell_W.setConstant(INVALID);
+  per_cell_W.row(infinity_cell).setZero();
+  std::queue<size_t> Q;
+  Q.push(infinity_cell);
+  while (!Q.empty()) {
+    size_t curr_cell = Q.front();
+    Q.pop();
+    for (const auto& neighbor : cell_adjacency[curr_cell]) {
+      size_t neighbor_cell, patch_idx;
+      bool direction;
+      std::tie(neighbor_cell, direction, patch_idx) = neighbor;
+      if ((per_cell_W.row(neighbor_cell).array() == INVALID).any()) {
+        per_cell_W.row(neighbor_cell) = per_cell_W.row(curr_cell);
+        for (size_t i=0; i<num_labels; i++) {
+          int inc = (patch_labels[patch_idx] == (int)i) ?
+            (direction ? -1:1) :0;
+          per_cell_W(neighbor_cell, i) =
+            per_cell_W(curr_cell, i) + inc;
+        }
+        Q.push(neighbor_cell);
+      } else {
 #ifndef NDEBUG
-                        int inc = direction ? -1:1;
-                        assert(per_cell_W(neighbor_cell, i) ==
-                                per_cell_W(curr_cell, i) + inc);
-#endif
-                    } else {
-                        assert(per_cell_W(neighbor_cell, i) ==
-                                per_cell_W(curr_cell, i));
-                    }
-                }
-            }
+        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);
+          } else {
+            assert(per_cell_W(neighbor_cell, i) ==
+                per_cell_W(curr_cell, i));
+          }
         }
+#endif
+      }
     }
+  }
 
-    W.resize(num_faces, num_labels*2);
-    for (size_t i=0; i<num_faces; i++) {
-        const size_t patch = P[i];
-        const size_t positive_cell = per_patch_cells(patch, 0);
-        const size_t negative_cell = per_patch_cells(patch, 1);
-        for (size_t j=0; j<num_labels; j++) {
-            W(i,j*2  ) = per_cell_W(positive_cell, j);
-            W(i,j*2+1) = per_cell_W(negative_cell, j);
-        }
+  W.resize(num_faces, num_labels*2);
+  for (size_t i=0; i<num_faces; i++) {
+    const size_t patch = P[i];
+    const size_t positive_cell = per_patch_cells(patch, 0);
+    const size_t negative_cell = per_patch_cells(patch, 1);
+    for (size_t j=0; j<num_labels; j++) {
+      W(i,j*2  ) = per_cell_W(positive_cell, j);
+      W(i,j*2+1) = per_cell_W(negative_cell, j);
     }
-
-    //for (size_t i=0; i<num_cells; i++) {
-    //    std::stringstream filename;
-    //    filename << "cell_" << i << "_w_" << per_cell_W(i, 1) << ".ply";
-    //    save_cell(filename.str(), i);
-    //}
+  }
 }
 
 

include/igl/copyleft/cgal/propagate_winding_numbers.h

@@ -24,32 +24,32 @@ namespace igl {
   {
     namespace cgal {
 
-        // Compute winding number on each side of the face.  The input mesh
-        // could contain multiple connected components.  The input mesh must
-        // represent the boundary of a valid 3D volume, which means it is
-        // closed, consistently oriented and induces integer winding numbers.
-        //
-        // Inputs:
-        //   V  #V by 3 list of vertex positions.
-        //   F  #F by 3 list of triangle indices into V.
-        //   labels  #F list of facet labels ranging from 0 to k-1.
-        //
-        // Output:
-        //   W  #F by k*2 list of winding numbers.  ``W(i,j*2)`` is the winding
-        //      number on the positive side of facet ``i`` with respect to the
-        //      facets labeled ``j``.  Similarly, ``W(i,j*2+1)`` is the winding
-        //      number on the negative side of facet ``i`` with respect to the
-        //      facets labeled ``j``.
-        template<
-            typename DerivedV,
-            typename DerivedF,
-            typename DerivedL,
-            typename DerivedW>
-        IGL_INLINE void propagate_winding_numbers(
-                const Eigen::PlainObjectBase<DerivedV>& V,
-                const Eigen::PlainObjectBase<DerivedF>& F,
-                const Eigen::PlainObjectBase<DerivedL>& labels,
-                Eigen::PlainObjectBase<DerivedW>& W);
+      // Compute winding number on each side of the face.  The input mesh
+      // could contain multiple connected components.  The input mesh must
+      // represent the boundary of a valid 3D volume, which means it is
+      // closed, consistently oriented and induces integer winding numbers.
+      //
+      // Inputs:
+      //   V  #V by 3 list of vertex positions.
+      //   F  #F by 3 list of triangle indices into V.
+      //   labels  #F list of facet labels ranging from 0 to k-1.
+      //
+      // Output:
+      //   W  #F by k*2 list of winding numbers.  ``W(i,j*2)`` is the winding
+      //      number on the positive side of facet ``i`` with respect to the
+      //      facets labeled ``j``.  Similarly, ``W(i,j*2+1)`` is the winding
+      //      number on the negative side of facet ``i`` with respect to the
+      //      facets labeled ``j``.
+      template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedL,
+        typename DerivedW>
+      IGL_INLINE void propagate_winding_numbers(
+          const Eigen::PlainObjectBase<DerivedV>& V,
+          const Eigen::PlainObjectBase<DerivedF>& F,
+          const Eigen::PlainObjectBase<DerivedL>& labels,
+          Eigen::PlainObjectBase<DerivedW>& W);
     }
   }
 }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -272,7 +272,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
 
     // Process un-touched faces.
     for (size_t i=0; i<num_faces; i++) {
-        if (!is_offending[i]) {
+        if (!is_offending[i] && !T[i].is_degenerate()) {
             resolved_faces.push_back(
                     { F(i,0), F(i,1), F(i,2) } );
             source_faces.push_back(i);

include/igl/doublearea.cpp

@@ -204,6 +204,8 @@ Eigen::PlainObjectBase<DeriveddblA> & dblA)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

include/igl/edge_lengths.cpp

@@ -72,6 +72,8 @@ IGL_INLINE void igl::edge_lengths(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&);
 // generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);

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>

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

@@ -7,8 +7,9 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "ambient_occlusion.h"
 #include "EmbreeIntersector.h"
-#include <igl/random_dir.h>
-#include <igl/EPS.h>
+#include "../Hit.h"
+#include "../random_dir.h"
+#include "../EPS.h"
 
 template <
   typename DerivedP,
@@ -44,7 +45,7 @@ IGL_INLINE void igl::embree::ambient_occlusion(
         // reverse ray
         d *= -1;
       }
-      igl::embree::Hit hit;
+      igl::Hit hit;
       const float tnear = 1e-4f;
       if(ei.intersectRay(origin,d,hit,tnear))
       {

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

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

include/igl/per_face_normals.cpp

@@ -105,6 +105,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, 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> >&);
 template void igl::per_face_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::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> >&);

include/igl/per_vertex_normals.cpp

@@ -107,4 +107,5 @@ template void igl::per_vertex_normals<Eigen::Matrix<double, -1, -1, 0, -1, -1>,
 template void igl::per_vertex_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&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_vertex_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&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::per_vertex_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&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::per_vertex_normals<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -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> >&);
 #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/resolve_duplicated_faces.cpp

@@ -21,14 +21,14 @@ IGL_INLINE void igl::resolve_duplicated_faces(
     Eigen::PlainObjectBase<DerivedF2>& F2,
     Eigen::PlainObjectBase<DerivedJ>& J) {
 
-  typedef typename DerivedF1::Scalar Index;
+  //typedef typename DerivedF1::Scalar Index;
   Eigen::VectorXi IA,IC;
   DerivedF1 uF;
   igl::unique_simplices(F1,uF,IA,IC);
 
   const size_t num_faces = F1.rows();
   const size_t num_unique_faces = uF.rows();
-  assert(IA.rows() == num_unique_faces);
+  assert((size_t) IA.rows() == num_unique_faces);
   // faces on top of each unique face
   std::vector<std::vector<int> > uF2F(num_unique_faces);
   // signed counts

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,120 @@
+// 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"
+
+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)
+  {
+    // 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))
+      {
+        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 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::__1::function<void (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, std::__1::vector<igl::Hit, std::__1::allocator<igl::Hit> >&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, std::__1::vector<igl::Hit, std::__1::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::__1::function<void (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, std::__1::vector<igl::Hit, std::__1::allocator<igl::Hit> >&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&, std::__1::vector<igl::Hit, std::__1::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::__1::function<void (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, std::__1::vector<igl::Hit, std::__1::allocator<igl::Hit> >&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::__1::vector<igl::Hit, std::__1::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

include/igl/unproject_ray.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_UNPROJECT_RAY_H
+#define IGL_UNPROJECT_RAY_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Construct a ray (source point + direction vector) given a screen space
+  // positions (e.g. mouse) and a model-view projection constellation.
+  //
+  // Inputs:
+  //   pos  2d screen-space position (x,y) 
+  //   model  4x4 model-view matrix
+  //   proj  4x4 projection matrix
+  //   viewport  4-long viewport vector
+  // Outputs:
+  //   s  source of ray (pos unprojected with z=0)
+  ///  dir  direction of ray (d - s) where d is pos unprojected with z=1
+  // 
+  template <
+    typename Derivedpos,
+    typename Derivedmodel,
+    typename Derivedproj,
+    typename Derivedviewport,
+    typename Deriveds,
+    typename Deriveddir>
+  IGL_INLINE void 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);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "unproject_ray.cpp"
+#endif
+#endif

include/igl/vertex_triangle_adjacency.cpp

@@ -50,4 +50,5 @@ template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, 3, 1, -1,
 template void igl::vertex_triangle_adjacency<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&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
 template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, long, long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&, std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >&);
+template void igl::vertex_triangle_adjacency<Eigen::Matrix<int, -1, -1, 0, -1, -1>, unsigned long, unsigned long>(Eigen::Matrix<int, -1, -1, 0, -1, -1>::Scalar, 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> > > >&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > >&);
 #endif

libigl-teaser.pdf.REMOVED.git-id

@@ -1 +1 @@
-6d633e0a6c52bd9130db774013a00a5dbe7adeb6
+d95997c253523e4aaf71fc6c570639f4bfafd399

scripts/update_gh-pages.sh

@@ -1,7 +1,6 @@
 #!/bin/bash
 # Usage: cd $LIBIGL; scripts/update_gh-pages.sh
 set -o xtrace
-git pull && git checkout gh-pages && git rebase master && git pull
 
 HEADER="title: libigl
 author: Alec Jacobson and Daniele Panozzo and others
@@ -14,8 +13,10 @@ html header:   <script type='text/javascript' src='http://cdn.mathjax.org/mathja
 "
 
 echo "$HEADER" \
-  | cat - README.md | multimarkdown -o index.html && \
-  git commit -m "update index.html to match README.md" README.md index.html
+  | cat - README.md | multimarkdown -o index.html
+
+echo "$HEADER" \
+  | cat - style-guidelines.md | multimarkdown -o style-guidelines.html 
 
 HEADER="title: libigl
 author: Alec Jacobson and Daniele Panozzo and others
@@ -28,12 +29,6 @@ html header:   <script type='text/javascript' src='http://cdn.mathjax.org/mathja
 "
 
 echo "$HEADER" \
-  | cat - google-soc/google-soc.md | multimarkdown -o google-soc/index.html && \
-  git commit -m "update google-soc/index.html to match google-soc/google-soc.md" google-soc/google-soc.md google-soc/index.html 
-
-echo "$HEADER" \
-  | cat - optional/README.md | multimarkdown -o optional/index.html && \
-  git commit -m "update index.html to match README.md" optional/README.md \
-  optional/index.html
+  | cat - optional/README.md | multimarkdown -o optional/index.html
 
-git push origin gh-pages && git checkout master && git merge gh-pages && git push
+multimarkdown tutorial/tutorial.md -o tutorial/tutorial.html

scripts/update_gh-pages_nogit.sh

@@ -1,31 +0,0 @@
-#!/bin/bash
-# Usage: cd $LIBIGL; scripts/update_gh-pages.sh
-set -o xtrace
-
-HEADER="title: libigl
-author: Alec Jacobson and Daniele Panozzo and others
-css: tutorial/style.css
-html header:   <script type='text/javascript' src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
-<link rel='stylesheet' href='http://yandex.st/highlightjs/7.3/styles/default.min.css'>
-<script src='http://yandex.st/highlightjs/7.3/highlight.min.js'></script>
-<script>hljs.initHighlightingOnLoad();</script>
-
-"
-
-echo "$HEADER" \
-  | cat - README.md | multimarkdown -o index.html
-
-HEADER="title: libigl
-author: Alec Jacobson and Daniele Panozzo and others
-css: ../tutorial/style.css
-html header:   <script type='text/javascript' src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
-<link rel='stylesheet' href='http://yandex.st/highlightjs/7.3/styles/default.min.css'>
-<script src='http://yandex.st/highlightjs/7.3/highlight.min.js'></script>
-<script>hljs.initHighlightingOnLoad();</script>
-
-"
-
-echo "$HEADER" \
-  | cat - optional/README.md | multimarkdown -o optional/index.html
-
-multimarkdown tutorial/tutorial.md -o tutorial/tutorial.html

style-guidelines.html

@@ -0,0 +1,316 @@
+<!DOCTYPE html>
+<html>
+<head>
+	<meta charset="utf-8"/>
+	<title>libigl</title>
+	<meta name="author" content="Alec Jacobson and Daniele Panozzo and others"/>
+	<link type="text/css" rel="stylesheet" href="tutorial/style.css"/>
+<script type='text/javascript' src='http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script>
+<link rel='stylesheet' href='http://yandex.st/highlightjs/7.3/styles/default.min.css'>
+<script src='http://yandex.st/highlightjs/7.3/highlight.min.js'></script>
+<script>hljs.initHighlightingOnLoad();</script>
+</head>
+<body>
+
+<h1 id="libiglstyleguidelines">Libigl Style Guidelines</h1>
+
+<p>This library is shared by many people. This document highlights some style
+guidelines for <em>developers</em> of the library, but also act as best-practices for
+users.</p>
+
+<blockquote>
+<p>This is still a work in progress and will likely grow in the near future.</p>
+</blockquote>
+
+<h2 id="filefunction">One function, one .h/.cpp pair</h2>
+
+<p>The structure of libigl is very flat and function-based. For every
+function/sub-routine create a single .h and .cpp file. For example, if you have
+a function that determines connected components from a face list <code>F</code> you would
+create the header <code>connected_components.h</code> and <code>connected_components.cpp</code> and the only
+function defined should be <code>void connected_components(const ... F, ... C)</code>. If the
+implementation of <code>connected_components</code> requires a subroutine to compute an
+adjacency matrix then <em>create another pair</em> <code>adjacency_matrix.h</code> and
+<code>adjacency_matrix.cpp</code> with a single function <code>void adjacency_matrix(const ... F, ...
+A)</code>.</p>
+
+<h3 id="example">Example</h3>
+
+<p>Here is an example function that would be defined in
+<code>include/igl/example_fun.h</code> and implemented in <code>include/igl/example_fun.cpp</code>.</p>
+
+<h4 id="example_fun.h"><code>example_fun.h</code></h4>
+
+<pre><code class="cpp">// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 [Your Name] [your email address]
+// 
+// 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_EXAMPLE_FUN_H
+#define IGL_EXAMPLE_FUN_H
+
+#include &quot;igl_inline.h&quot;
+
+namespace igl
+{
+  // This is an example of a function, it takes a templated parameter and
+  // shovels it into cout
+  //
+  // Templates:
+  //   T  type that supports
+  // Input:
+  //   input  some input of a Printable type
+  // Returns true for the sake of returning something
+  template &lt;typename Printable&gt;
+  IGL_INLINE bool example_fun(const Printable &amp; input);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include &quot;example_fun.cpp&quot;
+#endif
+
+#endif
+</code></pre>
+
+<h4 id="example_fun.cpp"><code>example_fun.cpp</code></h4>
+
+<pre><code>// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 [Your Name] [your email address]
+// 
+// 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 &quot;igl/example_fun.h&quot;
+#include &lt;iostream&gt;
+
+template &lt;typename Printable&gt;
+IGL_INLINE bool igl::example_fun(const Printable &amp; input)
+{
+  using namespace std;
+  cout&lt;&lt;&quot;example_fun: &quot;&lt;&lt;input&lt;&lt;endl;
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::example_fun&lt;double&gt;(const double&amp; input);
+template bool igl::example_fun&lt;int&gt;(const int&amp; input);
+#endif
+</code></pre>
+
+<h3 id="avoidstatichelperfunctions">Avoid static &#8220;helper&#8221; functions</h3>
+
+<p>Strive to encapsulate sub-functions that could possibly be useful outside of
+the implementation of your current function. This might mean abstracting the
+interface a bit. If it doesn&#8217;t dramatically effect performance then create a
+new pair of .h/.cpp files with this sub-function.</p>
+
+<h4 id="lambdafunctions">Lambda functions</h4>
+
+<p>If encapsulation in a separate file is not possible or does not make sense,
+then avoid crowding the namespace by creating lambda functions within the
+function implmentation.</p>
+
+<h3 id="avoidhelperclasses">Avoid &#8220;helper&#8221; classes</h3>
+
+<p>Libigl is built around the high-performance paradigm of &#8220;struct of arrays&#8221;
+rather than &#8220;array of structs&#8221;. The way we achieve this is to avoid classes and
+pass &#8220;basic types&#8221; directly. The price we pay is long function interfaces, but
+this increases code reuse dramatically. A &#8220;basic type&#8221; in our context is a
+Eigen type, stl type, or basic C type.</p>
+
+<h2 id="headerdocumentation">Header Documentation</h2>
+
+<p>Each function prototype should be well documented in its corresponding .h
+header file. A typical documentation consists of four parts:</p>
+
+<pre><code class="cpp">// [A human readable description of what the function does.]
+//
+// Inputs:
+//   [variable name of first (const) input]   [dimensions and
+//     description of this input variable]
+//   [variable name of second (const) input]   [dimensions and
+//     description of this input variable]
+//   ...
+// Outputs:
+//   [variable name of first output ]   [dimensions and
+//     description of this output variable]
+//   [variable name of second output ]   [dimensions and
+//     description of this output variable]
+//   ...
+// Returns [description of return value]
+</code></pre>
+
+<h3 id="example">Example</h3>
+
+<p>For example the header <code>barycenter.h</code></p>
+
+<pre><code>// Computes the barycenter of every simplex
+//
+// Inputs:
+//   V  #V x dim matrix of vertex coordinates
+//   F  #F x simplex_size  matrix of indices of simplex corners into V
+// Output:
+//   BC  #F x dim matrix of 3d vertices
+//
+</code></pre>
+
+<h2 id="constinputs">Const inputs</h2>
+
+<p>All input parameters should be demarcated <code>const</code>. If an input is also an
+output than consider exposing two parameters (one <code>const</code>) or be sure to list
+the variable under both <code>// Inputs:</code> and <code>// Outputs:</code> in the header comments.</p>
+
+<h2 id="referenceparameters">Reference parameters</h2>
+
+<p>All but simple types should be passed by reference (e.g. <code>Matrix &amp; mat</code>) rather
+than pointers (e.g. <code>Matrix * mat</code>) or value (e.g. <code>Matrix mat</code>).</p>
+
+<h2 id="returnsvsoutputparameters">Returns vs output parameters</h2>
+
+<p>All functions should be implemented with at least one overload that has a
+<code>void</code> or simple return type (e.g. <code>bool</code> on success/failure). With this
+implementation its then possible to write an overload that returns a single
+output.</p>
+
+<p>For example:</p>
+
+<pre><code class="cpp">template &lt;typename Atype&gt;
+void adjacency_matrix(const ... &amp; F, Eigen::SparseMatrix&lt;AType&gt; &amp; A);
+
+template &lt;typename Atype&gt;
+Eigen::SparseMatrix&lt;Atype&gt; adjacency_matrix(const ... &amp; F);
+</code></pre>
+
+<h2 id="functionnamingconventions">Function naming conventions</h2>
+
+<p>Functions (and <a href="#filefunction">thus also files</a>) should have simple,
+descriptive names using lowercase letters and underscores between words. Avoid
+unnecessary prefaces. For example, instead of <code>compute_adjacency_matrix</code>,
+<code>construct_adjacency_matrix</code>, <code>extract_adjacency_matrix</code>,
+<code>get_adjacency_matrix</code>, or <code>set_adjacency_matrix</code> just call the function
+<code>adjacency_matrix</code>.</p>
+
+<h2 id="variablenamingconventions">Variable naming conventions</h2>
+
+<p>Libigl prefers short (even single character) variable names <em>with heavy
+documentation</em> in the comments in the header file or above the declaration of
+the function. When possible use <code>V</code> to mean a list of vertex positions and <code>F</code>
+to mean a list of faces/triangles.</p>
+
+<h2 id="classnamingconventions">Class naming conventions</h2>
+
+<p>Classes should be avoided. When naming a class use CamelCase (e.g.
+SortableRow.h).</p>
+
+<h2 id="enumnamingconvertion">Enum naming convertion</h2>
+
+<p>Enums types should be placed in the appropriate <code>igl::</code> namespace and should be
+named in CamelCase (e.g. <code>igl::SolverStatus</code>) and instances should be named in
+ALL_CAPS with underscores between words and prefaced with the name of the enum.
+For example:</p>
+
+<pre><code class="cpp">namespace igl
+{
+  enum SolverStatus
+  {
+    // Good
+    SOLVER_STATUS_CONVERGED = 0,
+    // OK
+    SOLVER_STATUS_MAX_ITER = 1,
+    // Bad
+    SOLVER_STATUS_ERROR = 2,
+    NUM_SOLVER_STATUSES = 3,
+  };
+};
+</code></pre>
+
+<h3 id="exceptionforfileio">Exception for file IO</h3>
+
+<p>For legacy reasons, file reading and writing functions use a different naming
+convention. A functions reading a <code>.xyz</code> file should be named <code>readXYZ</code> and a
+function writing <code>.xyz</code> files should be names <code>writeXYZ</code>.</p>
+
+<h2 id="usingnamespace...inglobalscope"><code>using namespace ...</code> in global scope</h2>
+
+<p>Writing <code>using namespace std;</code>, <code>using namespace Eigen;</code> etc. outside of a
+global scope is strictly forbidden. Place these lines at the top of each
+function instead.</p>
+
+<h2 id="namespacesandexternaldependencies">Namespaces and external dependencies</h2>
+
+<p>Functions in the main library (directly in <code>include/igl</code>) should only depend on
+Eigen and stl. These functions should have the <code>igl::</code> namespace.</p>
+
+<p>Functions with other dependencies should be placed into
+appropriate sub-directories (e.g. if <code>myfunction</code> depends on tetgen then create
+<code>igl/copyleft/tetgen/myfunction.h</code> and <code>igl/copyleft/tetgen/myfunction.cpp</code> and give the function
+the namespace <code>igl::copyleft::tetgen::myfunction</code>.</p>
+
+<h3 id="copyleftsubdirectorynamespace">copyleft subdirectory/namespace</h3>
+
+<p>Dependencies that require users of libigl to release their projects open source
+(e.g. GPL) are considered aggressively &#8220;copyleft&#8221; and should be placed in the
+<code>include/igl/copyleft/</code> sub-directory and <code>igl::copyleft::</code> namespace.</p>
+
+<h2 id="assertions">Assertions</h2>
+
+<p>Be generous with assertions and always identify the assertion with strings:</p>
+
+<pre><code class="cpp">assert(m &lt; n &amp;&amp; &quot;m must be less than n&quot;);
+</code></pre>
+
+<h2 id="ifndefincludeguard">ifndef include guard</h2>
+
+<p>Every header file should be wrapped in an <code>#ifndef</code> compiler directive. The
+name of the guard should be in direct correspondence with the path of the .h
+file. For example, <code>include/igl/copyleft/tetgen/tetrahedralize.h</code> should be</p>
+
+<pre><code class="cpp">#ifndef IGL_COPYLEFT_TETGEN_TETRAHEDRALIZE_H
+#define IGL_COPYLEFT_TETGEN_TETRAHEDRALIZE_H
+...
+#endif
+</code></pre>
+
+<h2 id="spacesvs.tabsindentation">Spaces vs. tabs indentation</h2>
+
+<p>Do not use tabs. Use 2 spaces for each indentation level.</p>
+
+<h2 id="maxlinelength">Max line length</h2>
+
+<p>Limit lines to 80 characters. Break up long lines into many operations (this
+also helps performance).</p>
+
+<h2 id="includeorder">Include order</h2>
+
+<p><code>#include</code> directives at the top of a .h or .cpp file should be sorted
+according to a simple principle: place headers of files most likely to be
+edited by you first. This means for
+<code>include/igl/copyleft/tetgen/tetrahedralize.cpp</code> you might see</p>
+
+<pre><code class="cpp">// [Includes of headers in this directory]
+#include &quot;tetrahedralize.h&quot;
+#include &quot;mesh_to_tetgenio.h&quot;
+#include &quot;tetgenio_to_tetmesh.h&quot;
+// [Includes of headers in this project]
+#include &quot;../../matrix_to_list.h&quot;
+#include &quot;../../list_to_matrix.h&quot;
+#include &quot;../../boundary_facets.h&quot;
+// [Includes of headers of related projects]
+#include &lt;Eigen/Core&gt;
+// [Includes of headers of standard libraries]
+#include &lt;cassert&gt;
+#include &lt;iostream&gt;
+</code></pre>
+
+<h2 id="placementofincludes">Placement of includes</h2>
+
+<p>Whenever possible <code>#include</code> directives should be placed in the <code>.cpp</code>
+implementation file rather than the <code>.h</code> header file.</p>
+
+<h2 id="eigentemplates">Eigen templates</h2>
+
+</body>
+</html>

style-guidelines.md

@@ -1,22 +1,87 @@
 # Libigl Style Guidelines
 
-This library is shared by many people. This document highlights some style
-guidelines for _developers_ of the library, but also act as best-practices for
-users.
-
-> This is still a work in progress and will likely grow in the near future.
+Libigl is used and developed by many people. This document highlights some
+style guidelines for _developers_ of the library, but also acts as
+best-practices for users.
 
 ## One function, one .h/.cpp pair [filefunction]
 
 The structure of libigl is very flat and function-based. For every
-function/sub-routine create a single .h and .cpp file. For example, if you have
+function/sub-routine, create a single .h and .cpp file. For example, if you have
 a function that determines connected components from a face list `F` you would
 create the header `connected_components.h` and `connected_components.cpp` and the only
 function defined should be `void connected_components(const ... F, ... C)`. If the
 implementation of `connected_components` requires a subroutine to compute an
 adjacency matrix then _create another pair_ `adjacency_matrix.h` and
-`adjacency_matrix.cpp` with a single function `void adjacency_matrix(const ... F, ...
-A)`.
+`adjacency_matrix.cpp` with a single function `void adjacency_matrix(const ... F, ... A)`.
+
+### Example
+Here is an example function that would be defined in
+`include/igl/example_fun.h` and implemented in `include/igl/example_fun.cpp`.
+
+#### `example_fun.h`
+
+```cpp
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 [Your Name] [your email address]
+// 
+// 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_EXAMPLE_FUN_H
+#define IGL_EXAMPLE_FUN_H
+
+#include "igl_inline.h"
+
+namespace igl
+{
+  // This is an example of a function, it takes a templated parameter and
+  // shovels it into cout
+  //
+  // Templates:
+  //   T  type that supports
+  // Input:
+  //   input  some input of a Printable type
+  // Returns true for the sake of returning something
+  template <typename Printable>
+  IGL_INLINE bool example_fun(const Printable & input);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "example_fun.cpp"
+#endif
+
+#endif
+```
+
+#### `example_fun.cpp`
+
+```
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 [Your Name] [your email address]
+// 
+// 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 "igl/example_fun.h"
+#include <iostream>
+
+template <typename Printable>
+IGL_INLINE bool igl::example_fun(const Printable & input)
+{
+  using namespace std;
+  cout<<"example_fun: "<<input<<endl;
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::example_fun<double>(const double& input);
+template bool igl::example_fun<int>(const int& input);
+#endif
+```
+
 
 ### Avoid static "helper" functions
 
@@ -49,16 +114,16 @@ header file. A typical documentation consists of four parts:
 // [A human readable description of what the function does.]
 //
 // Inputs:
-//   [variable name of first (const) input]   [dimensions and
-//     description of this input variable]
-//   [variable name of second (const) input]   [dimensions and
-//     description of this input variable]
+//   [variable name of first (const) input]   [dimensions and description of
+//     this input variable]
+//   [variable name of second (const) input]   [dimensions and description of
+//     this input variable]
 //   ...
 // Outputs:
-//   [variable name of first output ]   [dimensions and
-//     description of this output variable]
-//   [variable name of second output ]   [dimensions and
-//     description of this output variable]
+//   [variable name of first output ]   [dimensions and description of this
+//     output variable]
+//   [variable name of second output ]   [dimensions and description of this
+//     output variable]
 //   ...
 // Returns [description of return value]
 ```
@@ -71,10 +136,10 @@ For example the header `barycenter.h`
 // Computes the barycenter of every simplex
 //
 // Inputs:
-//   V  #V x dim matrix of vertex coordinates
-//   F  #F x simplex_size  matrix of indices of simplex corners into V
+//   V  #V by dim matrix of vertex coordinates
+//   F  #F by simplex_size  matrix of indices of simplex corners into V
 // Output:
-//   BC  #F x dim matrix of 3d vertices
+//   BC  #F by dim matrix of 3d vertices
 //
 ```
 
@@ -233,7 +298,7 @@ edited by you first. This means for
 #include <iostream>
 ```
 
-## Eigen templates
-
-
+## Placement of includes
 
+Whenever possible `#include` directives should be placed in the `.cpp`
+implementation file rather than the `.h` header file.

style_guidelines.html

@@ -1,212 +0,0 @@
-<!DOCTYPE HTML>
-<html>
-  <head>
-    <title>libigl - Style Guidelines</title>
-    <link href="./style.css" rel="stylesheet" type="text/css">
-  <body>
-    <div id=container>
-    <div class=article_inner>
-    <a href=.><img src=libigl-logo.jpg alt="igl logo" class=center></a>
-    <h1>libigl Style Guidelines</h1>
-    <p>
-This library is shared by many people. This document highlights some style
-guidelines for <i>developers</i> of the <a href=.>libigl</a> library.
-   </p>
-   <p>
-Each function prototype should be well documented.  Write a summary of what
-the function does and a description of each template, input and output in
-each prototype. 
-    </p>
-
-    <h2>Example</h2>
-    <p>
-Here is an example function defined in <code>include/igl/example_fun.h</code> and
-implemented in <code>include/igl/example_fun.cpp</code>.
-    </p>
-    <h3>example_fun.h</h3>
-    <pre><code>
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2013 Alec Jacobson &lt;alecjacobson@gmail.com&gt;
-// 
-// 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_EXAMPLE_FUN_H
-#define IGL_EXAMPLE_FUN_H
-
-#include "igl_inline.h"
-
-namespace igl
-{
-  // This is an example of a function, it takes a templated parameter and
-  // shovels it into cout
-  //
-  // Templates:
-  //   T  type that supports
-  // Input:
-  //   input  some input of a Printable type
-  // Returns true for the sake of returning something
-  template &lt;typename Printable&gt;
-  IGL_INLINE bool example_fun(const Printable &amp; input);
-}
-
-#ifndef IGL_STATIC_LIBRARY
-#  include "example_fun.cpp"
-#endif
-
-#endif
-    </code></pre>
-    <h3>example_fun.cpp</h3>
-    <pre><code>
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2013 Alec Jacobson &lt;alecjacobson@gmail.com&gt;
-// 
-// 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 "igl/example_fun.h"
-#include &lt;iostream&gt;
-
-template &lt;typename Printable&gt;
-IGL_INLINE bool igl::example_fun(const Printable &amp; input)
-{
-  using namespace std;
-  cout&lt;&lt;"example_fun: "&lt;&lt;input&lt;&lt;endl;
-  return true;
-}
-
-#ifdef IGL_STATIC_LIBRARY
-template bool igl::example_fun&lt;double&gt;(const double&amp; input);
-template bool igl::example_fun&lt;int&gt;(const int&amp; input);
-#endif
-    </code></pre>
-
-
-
-    <h2 id=general_rules>General rules</h2>
-    <ul>
-      <li> Use a single .h/.cpp pair with the same name as the function </li>
-      <li>
-At least one version of the function should use references for all outputs
-      </li>
-      <li>
-Use wrappers and additional prototypes for returning single-output functions'
-outputs, but the default is to only use outputs
-      </li>
-      <li> All inputs should be <code>const</code>.</li>
-      <li>
-Use C++ references (e.g. <code>Matrix &amp; mat</code>) for inputs and outputs
-rather than pointers (e.g. <code>Matrix * mat</code>) or pass-by-copy (e.g.
-<code>Matrix mat</code>).
-      </li>
-      <li>
-Write multiple prototypes if you'd like to have optional parameters with
-default values.
-      </li>
-      <li>
-External dependencies (besides Eigen and OpenGL) are only allowed within extras.
-      </li>
-      <li>
-New extras and their dependencies must be discussed first.
-      </li>
-      <li>
-Hard-to-compile or obscure external dependencies can go in the <code>external/</code> directory.
-      </li>
-      <li>
-Do not use the <code>using namespace</code> directive anywhere outside a local scope. This
-means never write: <code>using namespace std;</code> or <code>using namespace
-igl;</code> etc. at the top of a file.
-      </li>
-      <li>
-Function names should either be the same as the corresponding MATLAB function
-or should use all lowercase separated by underscores: e.g.
-<code>my_function_name</code>
-      </li>
-      <li> Classes should be in <code>CamelCase</code></li>
-      <li> No tabs, only two spaces per indentation level </li>
-      <li> Be generous with assertions and always identify them with strings:
-      <br> e.g. <code>assert(m&lt;n &amp;&amp; "m must be less than n");</code></li>
-  </ul>
-
-    <h2>Specific rules</h2>
-    <p>
-Each file should be wrapped in an ifndef compiler directive. If the
-file's (and function's) name is example.h, then the ifndef should
-always begin with IGL_, then the function/file name in all caps then
-_H. As in:
-    </p>
-    <pre><code>
-#ifndef IGL_EXAMPLE_H
-#define IGL_EXAMPLE_H
-...
-#endif</code></pre>
-  
-    <p>
-Each file should begin with prototypes *without implementations* of
-each version of the function. All defined in the igl namespace. Each
-prototype should have its own comments describing what it is doing,
-and its templates, inputs, outputs.
-    </p>
-
-    <p>
-Implementation should be separated from prototypes in a .cpp file.
-    </p>
-
-    <p>
-Any includes, such as std libraries etc. should be in the .cpp file not the
-header .h file (whenever feasibly possible).
-    </p>
-
-    <p>
-Be generous by templating your inputs and outputs. If you do use
-templates, you must document the template just as you document inputs
-and outputs.
-    </p>
-
-    <h2>Useful scripts to check style</h2>
-  
-    <table>
-      <tr>
-        <th>script</th>
-        <th>Description</th>
-      <tr class=d0>
-        <td><code>grep -L IGL_INLINE *</code></td>
-        <td>Find files that aren't using "IGL_INLINE"</td>
-      </tr>
-      <tr class=d1>
-        <td><code>grep -L "namespace igl" *</code></td>
-        <td>Find files that aren't using igl namespace</td>
-      </tr>
-      <tr class=d0>
-        <td><code>grep -lP '\t' *</code></td>
-        <td>Find files using [TAB] character</td>
-      </tr>
-      <tr class=d1>
-        <td><code>grep -l "^using namespace" *.cpp</code></td>
-        <td>Find .cpp files that contain ^using namespace</td>
-      </tr>
-      <tr class=d0>
-        <td><code>grep -l 'assert([^"]*);' *</code></td>
-        <td>Find files using asserts without identifier strings</td>
-      </tr>
-      <tr class=d1>
-        <td><code>grep -l "ifndef IGL_.*[^H] *$" *.h</code></td>
-        <td>Find .h files that contain ifndef IGL_*[^H]</td>
-      </tr>
-      <tr class=d0>
-        <td><code>grep -L "^#ifndef IGL_" *</code></td>
-        <td>Find files that don't contain #ifndef IGL_</td>
-      </tr>
-      <tr class=d1>
-        <td><code>grep -L "^// This file is part of libigl" *</code></td>
-      </tr>
-    </table>
-
-    <p>See also: <a href=tutorial.html>tutorial</a>, <a
-    href=doc.html>auto-documentation</a>, <a href=file-formats/index.html>file
-    formats</a></p>
-  </div>
-  </div>
-  </body>
-</html>

tutorial/CMakeLists.txt

@@ -61,6 +61,12 @@ ENDIF(MSVC)
 
 option(ENABLE_TUTORIALS " " OFF)
 
+IF(EMBREE_FOUND)
+  add_subdirectory("../external/embree/" "embree")		
+  include_directories("../external/embree/include")		
+  list(APPEND SHARED_LIBRARIES "embree")		
+ENDIF(EMBREE_FOUND)
+
 # endif(True)
 
 # if(True)

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-5d18dcc0e2feb40028d7020f495e49d98bc1c2ad
+658461f27515ecf34dba2bb605a927096b74da75