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

Conflicts:
	.gitignore
	.gitmodules

Former-commit-id: 0bd05e387683e98f4efde2a85aa61c55dab17ba4

.gitignore

@@ -85,3 +85,5 @@ python/py_igl/todo
 python/__pycache__
 iglhelpers.pyc
 python/build2
+tests/build
+tests/bin

.gitmodules

@@ -9,3 +9,6 @@ url=https://github.com/libigl/nanogui.git
 [submodule "external/pybind11"]
 	path = external/pybind11
 	url = https://github.com/wjakob/pybind11.git
+[submodule "tests/googletest"]
+	path = tests/googletest
+	url = https://github.com/google/googletest.git

include/igl/boolean/mesh_boolean.cpp

@@ -7,6 +7,8 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "mesh_boolean.h"
 #include <igl/per_face_normals.h>
+#include <igl/boundary_facets.h>
+#include <igl/exterior_edges.h>
 #include <igl/cgal/peel_outer_hull_layers.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/remove_unreferenced.h>
@@ -187,13 +189,6 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     J = CJ;
     return;
   }
-  MatrixX3S N,CN;
-  per_face_normals_stable(V,F,N);
-  CN.resize(CF.rows(),3);
-  for(size_t f = 0;f<(size_t)CN.rows();f++)
-  {
-    CN.row(f) = N.row(CJ(f));
-  }
 
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"peel..."<<endl;
@@ -202,7 +197,7 @@ IGL_INLINE void igl::boolean::mesh_boolean(
   // peel layers keeping track of odd and even flips
   VectorXi I;
   Matrix<bool,Dynamic,1> flip;
-  peel_outer_hull_layers(EV,CF,CN,I,flip);
+  peel_outer_hull_layers(EV,CF,I,flip);
   // 0 is "first" iteration, so it's odd
   Array<bool,Dynamic,1> odd = igl::mod(I,2).array()==0;
 
@@ -251,6 +246,14 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     G.row(g) = Gflip[g] ? CF.row(vG[g]).reverse().eval() : CF.row(vG[g]);
     GJ(g) = CJ(vG[g]);
   }
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  {
+    MatrixXd O;
+    boundary_facets(FC,O);
+    cout<<"# boundary: "<<O.rows()<<endl;
+  }
+  cout<<"# exterior: "<<exterior_edges(FC).rows()<<endl;
+#endif
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"clean..."<<endl;
 #endif
@@ -264,6 +267,7 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     vector<vector<Index> > uG2G(uG.rows());
     // signed counts
     VectorXi counts = VectorXi::Zero(uG.rows());
+    VectorXi ucounts = VectorXi::Zero(uG.rows());
     // loop over all faces
     for(Index g = 0;g<gm;g++)
     {
@@ -276,6 +280,7 @@ IGL_INLINE void igl::boolean::mesh_boolean(
         (G(g,0) == uG(ug,1) && G(g,1) == uG(ug,2) && G(g,2) == uG(ug,0)) ||
         (G(g,0) == uG(ug,2) && G(g,1) == uG(ug,0) && G(g,2) == uG(ug,1));
       counts(ug) += consistent ? 1 : -1;
+      ucounts(ug)++;
     }
     MatrixX3I oldG = G;
     // Faces of output vG[i] = j means ith face of output should be jth face in
@@ -285,17 +290,33 @@ IGL_INLINE void igl::boolean::mesh_boolean(
     {
       // if signed occurrences is zero or ±two then keep none
       // else if signed occurrences is ±one then keep just one facet
-      if(abs(counts(ug)) == 1)
+      switch(abs(counts(ug)))
       {
-        assert(uG2G.size() > 0);
-        vG.push_back(uG2G[ug][0]);
-      }
+        case 1:
+          assert(uG2G[ug].size() > 0);
+          vG.push_back(uG2G[ug][0]);
 #ifdef IGL_MESH_BOOLEAN_DEBUG
-      else
-      {
-        cout<<"Skipping "<<uG2G.size()<<" facets..."<<endl;
-      }
+          if(abs(ucounts(ug)) != 1)
+          {
+            cout<<"count,ucount of "<<counts(ug)<<","<<ucounts(ug)<<endl;
+          }
 #endif
+          break;
+        case 0:
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+          cout<<"Skipping "<<uG2G[ug].size()<<" facets..."<<endl;
+          if(abs(ucounts(ug)) != 0)
+          {
+            cout<<"count,ucount of "<<counts(ug)<<","<<ucounts(ug)<<endl;
+          }
+#endif
+          break;
+        default:
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+          cout<<"Didn't expect to be here."<<endl;
+#endif
+          assert(false && "Shouldn't count be -1/0/1 ?");
+      }
     }
     G.resize(vG.size(),3);
     J.resize(vG.size());
@@ -311,6 +332,14 @@ IGL_INLINE void igl::boolean::mesh_boolean(
   //cerr<<"warning not removing unref"<<endl;
   //VC = CV;
   //FC = G;
+#ifdef IGL_MESH_BOOLEAN_DEBUG
+  {
+    MatrixXd O;
+    boundary_facets(FC,O);
+    cout<<"# boundary: "<<O.rows()<<endl;
+  }
+  cout<<"# exterior: "<<exterior_edges(FC).rows()<<endl;
+#endif
 }
 
 #ifdef IGL_STATIC_LIBRARY
@@ -319,9 +348,10 @@ IGL_INLINE void igl::boolean::mesh_boolean(
 #include <igl/remove_unreferenced.cpp>
 template void igl::remove_unreferenced<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #include <igl/cgal/peel_outer_hull_layers.cpp>
-template unsigned long igl::cgal::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<bool, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template unsigned long
+igl::cgal::peel_outer_hull_layers<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 #include <igl/cgal/outer_hull.cpp>
-template void igl::cgal::outer_hull<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
+template void igl::cgal::outer_hull<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Lazy_exact_nt<CGAL::Gmpq>, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 // Explicit template specialization
 template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::boolean::mesh_boolean<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::boolean::MeshBooleanType const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/cgal/order_facets_around_edge.cpp

@@ -0,0 +1,211 @@
+#include "order_facets_around_edge.h"
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+namespace igl {
+    namespace cgal {
+        namespace order_facets_around_edges_helper {
+            template<typename T>
+            std::vector<size_t> index_sort(const std::vector<T>& data) {
+                const size_t len = data.size();
+                std::vector<size_t> order(len);
+                for (size_t i=0; i<len; i++) order[i] = i;
+
+                auto comp = [&](size_t i, size_t j) {
+                    return data[i] < data[j];
+                };
+                std::sort(order.begin(), order.end(), comp);
+                return order;
+            }
+        }
+    }
+}
+
+// adj_faces contains signed index starting from +- 1.
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedI >
+void igl::cgal::order_facets_around_edge(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    size_t s, size_t d, 
+    const std::vector<int>& adj_faces,
+    Eigen::PlainObjectBase<DerivedI>& order) {
+
+    using namespace igl::cgal::order_facets_around_edges_helper;
+
+    // Although we only need exact predicates in the algorithm,
+    // exact constructions are needed to avoid degeneracies due to
+    // casting to double.
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+    typedef K::Point_3 Point_3;
+    typedef K::Plane_3 Plane_3;
+
+    auto get_face_index = [&](int adj_f)->size_t{
+        return abs(adj_f) - 1;
+    };
+
+    auto get_opposite_vertex = [&](size_t fid)->size_t {
+        if (F(fid, 0) != s && F(fid, 0) != d) return F(fid, 0);
+        if (F(fid, 1) != s && F(fid, 1) != d) return F(fid, 1);
+        if (F(fid, 2) != s && F(fid, 2) != d) return F(fid, 2);
+        assert(false);
+        return -1;
+    };
+
+    // Handle base cases
+    if (adj_faces.size() == 0) {
+        order.resize(0, 1);
+        return;
+    } else if (adj_faces.size() == 1) {
+        order.resize(1, 1);
+        order(0, 0) = 0;
+        return;
+    } else if (adj_faces.size() == 2) {
+        const size_t o1 =
+            get_opposite_vertex(get_face_index(adj_faces[0]));
+        const size_t o2 =
+            get_opposite_vertex(get_face_index(adj_faces[1]));
+        const Point_3 ps(V(s, 0), V(s, 1), V(s, 2));
+        const Point_3 pd(V(d, 0), V(d, 1), V(d, 2));
+        const Point_3 p1(V(o1, 0), V(o1, 1), V(o1, 2));
+        const Point_3 p2(V(o2, 0), V(o2, 1), V(o2, 2));
+        order.resize(2, 1);
+        switch (CGAL::orientation(ps, pd, p1, p2)) {
+            case CGAL::POSITIVE:
+                order(0, 0) = 1;
+                order(1, 0) = 0;
+                break;
+            case CGAL::NEGATIVE:
+                order(0, 0) = 0;
+                order(1, 0) = 1;
+                break;
+            case CGAL::COPLANAR:
+                order(0, 0) = adj_faces[0] < adj_faces[1] ? 0:1;
+                order(1, 0) = adj_faces[0] < adj_faces[1] ? 1:0;
+                break;
+            default:
+                assert(false);
+        }
+        return;
+    }
+
+    const size_t num_adj_faces = adj_faces.size();
+    const size_t o = get_opposite_vertex(
+            get_face_index(adj_faces[0]));
+    const Point_3 p_s(V(s, 0), V(s, 1), V(s, 2));
+    const Point_3 p_d(V(d, 0), V(d, 1), V(d, 2));
+    const Point_3 p_o(V(o, 0), V(o, 1), V(o, 2));
+    const Plane_3 separator(p_s, p_d, p_o);
+    assert(!separator.is_degenerate());
+
+    std::vector<Point_3> opposite_vertices;
+    for (size_t i=0; i<num_adj_faces; i++) {
+        const size_t o = get_opposite_vertex(
+                get_face_index(adj_faces[i]));
+        opposite_vertices.emplace_back(
+                V(o, 0), V(o, 1), V(o, 2));
+    }
+
+    std::vector<int> positive_side;
+    std::vector<int> negative_side;
+    std::vector<int> tie_positive_oriented;
+    std::vector<int> tie_negative_oriented;
+
+    std::vector<size_t> positive_side_index;
+    std::vector<size_t> negative_side_index;
+    std::vector<size_t> tie_positive_oriented_index;
+    std::vector<size_t> tie_negative_oriented_index;
+
+    for (size_t i=0; i<num_adj_faces; i++) {
+        const int f = adj_faces[i];
+        const Point_3& p_a = opposite_vertices[i];
+        auto orientation = separator.oriented_side(p_a);
+        switch (orientation) {
+            case CGAL::ON_POSITIVE_SIDE:
+                positive_side.push_back(f);
+                positive_side_index.push_back(i);
+                break;
+            case CGAL::ON_NEGATIVE_SIDE:
+                negative_side.push_back(f);
+                negative_side_index.push_back(i);
+                break;
+            case CGAL::ON_ORIENTED_BOUNDARY:
+                {
+                    const Plane_3 other(p_s, p_d, p_a);
+                    const auto target_dir = separator.orthogonal_direction();
+                    const auto query_dir = other.orthogonal_direction();
+                    if (target_dir == query_dir) {
+                        tie_positive_oriented.push_back(f);
+                        tie_positive_oriented_index.push_back(i);
+                    } else if (target_dir == -query_dir) {
+                        tie_negative_oriented.push_back(f);
+                        tie_negative_oriented_index.push_back(i);
+                    } else {
+                        assert(false);
+                    }
+                }
+                break;
+            default:
+                // Should not be here.
+                assert(false);
+        }
+    }
+
+    Eigen::PlainObjectBase<DerivedI> positive_order, negative_order;
+    order_facets_around_edge(V, F, s, d, positive_side, positive_order);
+    order_facets_around_edge(V, F, s, d, negative_side, negative_order);
+    std::vector<size_t> tie_positive_order =
+        index_sort(tie_positive_oriented);
+    std::vector<size_t> tie_negative_order =
+        index_sort(tie_negative_oriented);
+
+    // Copy results into order vector.
+    const size_t tie_positive_size = tie_positive_oriented.size();
+    const size_t tie_negative_size = tie_negative_oriented.size();
+    const size_t positive_size = positive_order.size();
+    const size_t negative_size = negative_order.size();
+
+    order.resize(tie_positive_size + positive_size +
+            tie_negative_size + negative_size, 1);
+
+    size_t count=0;
+    for (size_t i=0; i<tie_positive_size; i++) {
+        order(count+i, 0) =
+            tie_positive_oriented_index[tie_positive_order[i]];
+    }
+    count += tie_positive_size;
+
+    for (size_t i=0; i<negative_size; i++) {
+        order(count+i, 0) = negative_side_index[negative_order(i, 0)];
+    }
+    count += negative_size;
+
+    for (size_t i=0; i<tie_negative_size; i++) {
+        order(count+i, 0) =
+            tie_negative_oriented_index[tie_negative_order[i]];
+    }
+    count += tie_negative_size;
+
+    for (size_t i=0; i<positive_size; i++) {
+        order(count+i, 0) = positive_side_index[positive_order(i, 0)];
+    }
+    count += positive_size;
+    assert(count == num_adj_faces);
+
+    // Find the correct start point.
+    size_t start_idx = 0;
+    for (size_t i=0; i<num_adj_faces; i++) {
+        const Point_3& p_a = opposite_vertices[order(i, 0)];
+        const Point_3& p_b =
+            opposite_vertices[order((i+1)%num_adj_faces, 0)];
+        if (CGAL::orientation(p_s, p_d, p_a, p_b) == CGAL::POSITIVE) {
+            start_idx = (i+1)%num_adj_faces;
+            break;
+        }
+    }
+    DerivedI circular_order = order;
+    for (size_t i=0; i<num_adj_faces; i++) {
+        order(i, 0) = circular_order((start_idx + i)%num_adj_faces, 0);
+    }
+}

include/igl/cgal/order_facets_around_edge.h

@@ -0,0 +1,53 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+//
+#ifndef ORDER_FACETS_AROUND_EDGE_H
+#define ORDER_FACETS_AROUND_EDGE_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl {
+    namespace cgal {
+        // Given a directed edge, sort its adjacent faces.  Assuming the
+        // directed edge is (s, d).  Sort the adjacent faces clockwise around the
+        // axis (d - s), i.e. left-hand rule.  An adjacent face is consistently
+        // oriented if it contains (d, s) as a directed edge.
+        //
+        // For overlapping faces, break the tie using signed face index, smaller
+        // signed index comes before the larger signed index.  Signed index is
+        // computed as (consistent? 1:-1) * (face_index + 1).
+        //
+        // Inputs:
+        //   V          #V by 3 list of vertices.
+        //   F          #F by 3 list of faces
+        //   s          Index of source vertex.
+        //   d          Index of desination vertex.
+        //   adj_faces  List of adjacent face signed indices.
+        //
+        // Output:
+        //   order      List of face indices that orders adjacent faces around
+        //              edge (s, d) clockwise.
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedI >
+        IGL_INLINE
+        void order_facets_around_edge(
+                const Eigen::PlainObjectBase<DerivedV>& V,
+                const Eigen::PlainObjectBase<DerivedF>& F,
+                size_t s, size_t d, 
+                const std::vector<int>& adj_faces,
+                Eigen::PlainObjectBase<DerivedI>& order);
+    }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "order_facets_around_edge.cpp"
+#endif
+#endif

include/igl/cgal/order_facets_around_edges.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 "order_facets_around_edges.h"
+#include "order_facets_around_edge.h"
 #include "../sort_angles.h"
 #include <type_traits>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
@@ -252,3 +253,80 @@ igl::cgal::order_facets_around_edges(
         }
     }
 }
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename uE2oEType,
+    typename uE2CType >
+IGL_INLINE void igl::cgal::order_facets_around_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedE>& E,
+        const Eigen::PlainObjectBase<DeriveduE>& uE,
+        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+        const std::vector<std::vector<uE2EType> >& uE2E,
+        std::vector<std::vector<uE2oEType> >& uE2oE,
+        std::vector<std::vector<uE2CType > >& uE2C ) {
+
+    typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
+    const size_t num_faces = F.rows();
+    const size_t num_undirected_edges = uE.rows();
+
+    auto edge_index_to_face_index = [&](size_t ei) { return ei % num_faces; };
+    auto edge_index_to_corner_index = [&](size_t ei) { return ei / num_faces; };
+
+    uE2oE.resize(num_undirected_edges);
+    uE2C.resize(num_undirected_edges);
+
+    for(size_t ui = 0;ui<num_undirected_edges;ui++)
+    {
+        const auto& adj_edges = uE2E[ui];
+        const size_t edge_valance = adj_edges.size();
+        assert(edge_valance > 0);
+
+        const auto ref_edge = adj_edges[0];
+        const auto ref_face = edge_index_to_face_index(ref_edge);
+
+        const auto ref_corner_o = edge_index_to_corner_index(ref_edge);
+        const auto ref_corner_s = (ref_corner_o+1)%3;
+        const auto ref_corner_d = (ref_corner_o+2)%3;
+
+        const typename DerivedF::Scalar o = F(ref_face, ref_corner_o);
+        const typename DerivedF::Scalar s = F(ref_face, ref_corner_s);
+        const typename DerivedF::Scalar d = F(ref_face, ref_corner_d);
+
+        std::vector<bool> cons(edge_valance);
+        std::vector<int> adj_faces(edge_valance);
+        for (size_t fei=0; fei<edge_valance; fei++) {
+            const auto fe = adj_edges[fei];
+            const auto f = edge_index_to_face_index(fe);
+            const auto c = edge_index_to_corner_index(fe);
+            cons[fei] = (d == F(f, (c+1)%3));
+            adj_faces[fei] = (f+1) * (cons[fei] ? 1:-1);
+
+            assert( cons[fei] ||  (d == F(f,(c+2)%3)));
+            assert(!cons[fei] || (s == F(f,(c+2)%3)));
+            assert(!cons[fei] || (d == F(f,(c+1)%3)));
+        }
+
+        Eigen::VectorXi order;
+        order_facets_around_edge(V, F, s, d, adj_faces, order);
+        assert(order.size() == edge_valance);
+
+        auto& ordered_edges = uE2oE[ui];
+        auto& consistency = uE2C[ui];
+
+        ordered_edges.resize(edge_valance);
+        consistency.resize(edge_valance);
+        for (size_t fei=0; fei<edge_valance; fei++) {
+            ordered_edges[fei] = adj_edges[order[fei]];
+            consistency[fei] = cons[order[fei]];
+        }
+    }
+}
+

include/igl/cgal/order_facets_around_edges.h

@@ -10,12 +10,20 @@
 #include "../igl_inline.h"
 #include <Eigen/Core>
 #include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <vector>
 
 namespace igl
 {
   namespace cgal
   {
-    // For each undirected edge, sort its adjacent faces.
+    // For each undirected edge, sort its adjacent faces.  Assuming the
+    // undirected edge is (s, d).  Sort the adjacent faces clockwise around the
+    // axis (d - s), i.e. left-hand rule.  An adjacent face is consistently
+    // oriented if it contains (d, s) as a directed edge.
+    //
+    // For overlapping faces, break the tie using signed face index, smaller
+    // signed index comes before the larger signed index.  Signed index is
+    // computed as (consistent? 1:-1) * index.
     //
     // Inputs:
     //   V    #V by 3 list of vertices.
@@ -79,6 +87,27 @@ namespace igl
             const std::vector<std::vector<uE2EType> >& uE2E,
             std::vector<std::vector<uE2oEType> >& uE2oE,
             std::vector<std::vector<uE2CType > >& uE2C );
+
+    // Order faces around each edge. Only exact predicate is used in the algorithm.
+    // Normal is not needed.
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedE,
+        typename DeriveduE,
+        typename DerivedEMAP,
+        typename uE2EType,
+        typename uE2oEType,
+        typename uE2CType >
+    IGL_INLINE void order_facets_around_edges(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DerivedE>& E,
+            const Eigen::PlainObjectBase<DeriveduE>& uE,
+            const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+            const std::vector<std::vector<uE2EType> >& uE2E,
+            std::vector<std::vector<uE2oEType> >& uE2oE,
+            std::vector<std::vector<uE2CType > >& uE2C );
   }
 }
 

include/igl/cgal/outer_facet.cpp

@@ -0,0 +1,74 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "outer_facet.h"
+#include "../outer_element.h"
+#include "order_facets_around_edge.h"
+#include <algorithm>
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedI,
+    typename IndexType
+    >
+IGL_INLINE void igl::cgal::outer_facet(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+        const Eigen::PlainObjectBase<DerivedI> & I,
+        IndexType & f,
+        bool & flipped) {
+
+    // Algorithm:
+    //
+    //    1. Find an outer edge (s, d).
+    //
+    //    2. Order adjacent facets around this edge. Because the edge is an
+    //    outer edge, there exists a plane passing through it such that all its
+    //    adjacent facets lie on the same side. The implementation of
+    //    order_facets_around_edge() will find a natural start facet such that
+    //    The first and last facets according to this order are on the outside.
+    //
+    //    3. Because the vertex s is an outer vertex by construction (see
+    //    implemnetation of outer_edge()). The first adjacent facet is facing
+    //    outside (i.e. flipped=false) if it contains directed edge (s, d).  
+    //
+    typedef typename DerivedV::Scalar Scalar; typedef typename DerivedV::Index
+        Index; const size_t INVALID = std::numeric_limits<size_t>::max();
+
+    Index s,d;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> incident_faces;
+    outer_edge(V, F, I, s, d, incident_faces);
+    assert(incident_faces.size() > 0);
+
+    auto convert_to_signed_index = [&](size_t fid) -> int{
+        if ((F(fid, 0) == s && F(fid, 1) == d) ||
+            (F(fid, 1) == s && F(fid, 2) == d) ||
+            (F(fid, 2) == s && F(fid, 0) == d) ) {
+            return int(fid+1) * -1;
+        } else {
+            return int(fid+1);
+        }
+    };
+
+    auto signed_index_to_index = [&](int signed_id) -> size_t {
+        return size_t(abs(signed_id) - 1);
+    };
+
+    std::vector<int> adj_faces(incident_faces.size());
+    std::transform(incident_faces.data(),
+            incident_faces.data() + incident_faces.size(),
+            adj_faces.begin(),
+            convert_to_signed_index);
+
+    Eigen::VectorXi order;
+    order_facets_around_edge(V, F, s, d, adj_faces, order);
+
+    f = signed_index_to_index(adj_faces[order[0]]);
+    flipped = adj_faces[order[0]] > 0;
+}

include/igl/cgal/outer_facet.h

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_CGAL_OUTER_FACET_H
+#define IGL_CGAL_OUTER_FACET_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+    namespace cgal
+    {
+        // Find a facet that is reachable from infinity without crossing any faces.
+        // Such facet is called "outer facet."
+        //
+        // Precondition: The input mesh must have all self-intersection resolved.  I.e
+        // there is no duplicated vertices, no overlapping edge and no intersecting
+        // faces (the only exception is there could be topologically duplicated faces).
+        // See cgal::remesh_self_intersections.h for how to obtain such input.
+        //
+        // This function differ from igl::outer_facet() in the fact this
+        // funciton is more robust because it does not rely on facet normals.
+        //
+        // Inputs:
+        //   V  #V by 3 list of vertex positions
+        //   F  #F by 3 list of triangle indices into V
+        //   I  #I list of facets to consider
+        // Outputs:
+        //   f  Index of the outer facet.
+        //   flipped  true iff the normal of f points inwards.
+        template<
+            typename DerivedV,
+            typename DerivedF,
+            typename DerivedI,
+            typename IndexType
+            >
+        IGL_INLINE void outer_facet(
+                const Eigen::PlainObjectBase<DerivedV> & V,
+                const Eigen::PlainObjectBase<DerivedF> & F,
+                const Eigen::PlainObjectBase<DerivedI> & I,
+                IndexType & f,
+                bool & flipped);
+    }
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "outer_facet.cpp"
+#endif
+#endif

include/igl/cgal/outer_hull.cpp

@@ -7,7 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "outer_hull.h"
 #include "order_facets_around_edges.h"
-#include "../outer_element.h"
+#include "outer_facet.h"
 #include "../sortrows.h"
 #include "../facet_components.h"
 #include "../winding_number.h"
@@ -30,14 +30,12 @@
 template <
   typename DerivedV,
   typename DerivedF,
-  typename DerivedN,
   typename DerivedG,
   typename DerivedJ,
   typename Derivedflip>
 IGL_INLINE void igl::cgal::outer_hull(
   const Eigen::PlainObjectBase<DerivedV> & V,
   const Eigen::PlainObjectBase<DerivedF> & F,
-  const Eigen::PlainObjectBase<DerivedN> & N,
   Eigen::PlainObjectBase<DerivedG> & G,
   Eigen::PlainObjectBase<DerivedJ> & J,
   Eigen::PlainObjectBase<Derivedflip> & flip)
@@ -53,7 +51,6 @@ IGL_INLINE void igl::cgal::outer_hull(
   typedef Matrix<typename DerivedF::Scalar,Dynamic,DerivedF::ColsAtCompileTime> MatrixXF;
   typedef Matrix<typename DerivedG::Scalar,Dynamic,DerivedG::ColsAtCompileTime> MatrixXG;
   typedef Matrix<typename DerivedJ::Scalar,Dynamic,DerivedJ::ColsAtCompileTime> MatrixXJ;
-  typedef Matrix<typename DerivedN::Scalar,1,3> RowVector3N;
   const Index m = F.rows();
 
   // UNUSED:
@@ -94,7 +91,7 @@ IGL_INLINE void igl::cgal::outer_hull(
 
   std::vector<std::vector<typename DerivedF::Index> > uE2oE;
   std::vector<std::vector<bool> > uE2C;
-  order_facets_around_edges(V, F, N, E, uE, EMAP, uE2E, uE2oE, uE2C);
+  order_facets_around_edges(V, F, E, uE, EMAP, uE2E, uE2oE, uE2C);
   uE2E = uE2oE;
   VectorXI diIM(3*m);
   for (auto ue : uE2E) {
@@ -160,7 +157,7 @@ IGL_INLINE void igl::cgal::outer_hull(
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"outer facet..."<<endl;
 #endif
-    outer_facet(V,F,N,IM,f,f_flip);
+  igl::cgal::outer_facet(V,F,IM,f,f_flip);
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"outer facet: "<<f<<endl;
   //cout << V.row(F(f, 0)) << std::endl;
@@ -223,19 +220,6 @@ IGL_INLINE void igl::cgal::outer_hull(
               << uE2E[EMAP(e)][i] % m * (uE2C[EMAP(e)][i] ? 1:-1) << ")" << std::endl;
       }
 #endif
-      //// find overlapping face-edges
-      //const auto & neighbors = uE2E[EMAP(e)];
-      //// normal after possible flipping
-      //const auto & fN = (flip(f)?-1.:1.)*N.row(f);
-      //// Edge vector according to f's (flipped) orientation.
-      ////const auto & eV = (V.row(fd)-V.row(fs)).normalized();
-
-//#warning "EXPERIMENTAL, DO NOT USE"
-      //// THIS IS WRONG! The first face is---after sorting---no longer the face
-      //// used for orienting the sort.
-      //const auto ui = EMAP(e);
-      //const auto fe0 = uE2E[ui][0];
-      //const auto es = F(fe0%m,((fe0/m)+1)%3);
 
       // is edge consistent with edge of face used for sorting
       const int e_cons = (uE2C[EMAP(e)][diIM(e)] ? 1: -1);
@@ -245,11 +229,6 @@ IGL_INLINE void igl::cgal::outer_hull(
       {
         const int nfei_new = (diIM(e) + 2*val + e_cons*step*(flip(f)?-1:1))%val;
         const int nf = uE2E[EMAP(e)][nfei_new] % m;
-        // Don't consider faces with identical dihedral angles
-        //if ((di[EMAP(e)][diIM(e)].array() != di[EMAP(e)][nfei_new].array()).any())
-        //if((di[EMAP(e)][diIM(e)] != di[EMAP(e)][nfei_new]))
-//#warning "THIS IS HACK, FIX ME"
-//        if( abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new]) < 1e-15 )
         {
 #ifdef IGL_OUTER_HULL_DEBUG
         //cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
@@ -281,64 +260,6 @@ IGL_INLINE void igl::cgal::outer_hull(
       }
 
       int max_ne = -1;
-      //// Loop over and find max dihedral angle
-      //typename DerivedV::Scalar max_di = -1;
-      //for(const auto & ne : neighbors)
-      //{
-      //  const int nf = ne%m;
-      //  if(nf == f)
-      //  {
-      //    continue;
-      //  }
-      //  // Corner of neighbor
-      //  const int nc = ne/m;
-      //  // Is neighbor oriented consistently with (flipped) f?
-      //  //const int ns = F(nf,(nc+1)%3);
-      //  const int nd = F(nf,(nc+2)%3);
-      //  const bool cons = (flip(f)?fd:fs) == nd;
-      //  // Normal after possibly flipping to match flip or orientation of f
-      //  const auto & nN = (cons? (flip(f)?-1:1.) : (flip(f)?1.:-1.) )*N.row(nf);
-      //  // Angle between n and f
-      //  const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN));
-      //  if(ndi>=max_di)
-      //  {
-      //    max_ne = ne;
-      //    max_di = ndi;
-      //  }
-      //}
-
-      ////cout<<(max_ne != max_ne_2)<<" =?= "<<e_cons<<endl;
-      //if(max_ne != max_ne_2)
-      //{
-      //  cout<<(f+1)<<" ---> "<<(max_ne%m)+1<<" != "<<(max_ne_2%m)+1<<" ... "<<e_cons<<" "<<flip(f)<<endl;
-      //  typename DerivedV::Scalar max_di = -1;
-      //  for(size_t nei = 0;nei<neighbors.size();nei++)
-      //  {
-      //    const auto & ne = neighbors[nei];
-      //    const int nf = ne%m;
-      //    if(nf == f)
-      //    {
-      //      cout<<"  "<<(ne%m)+1<<":\t"<<0<<"\t"<<di[EMAP[e]][nei]<<" "<<diIM(ne)<<endl;
-      //      continue;
-      //    }
-      //    // Corner of neighbor
-      //    const int nc = ne/m;
-      //    // Is neighbor oriented consistently with (flipped) f?
-      //    //const int ns = F(nf,(nc+1)%3);
-      //    const int nd = F(nf,(nc+2)%3);
-      //    const bool cons = (flip(f)?fd:fs) == nd;
-      //    // Normal after possibly flipping to match flip or orientation of f
-      //    const auto & nN = (cons? (flip(f)?-1:1.) : (flip(f)?1.:-1.) )*N.row(nf);
-      //    // Angle between n and f
-      //    const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN));
-      //    cout<<"  "<<(ne%m)+1<<":\t"<<ndi<<"\t"<<di[EMAP[e]][nei]<<" "<<diIM(ne)<<endl;
-      //    if(ndi>=max_di)
-      //    {
-      //      max_ne = ne;
-      //      max_di = ndi;
-      //    }
-      //  }
-      //}
       if(nfei >= 0)
       {
         max_ne = uE2E[EMAP(e)][nfei];
@@ -513,33 +434,9 @@ IGL_INLINE void igl::cgal::outer_hull(
   }
 }
 
-template <
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedG,
-  typename DerivedJ,
-  typename Derivedflip>
-IGL_INLINE void igl::cgal::outer_hull(
-  const Eigen::PlainObjectBase<DerivedV> & V,
-  const Eigen::PlainObjectBase<DerivedF> & F,
-  Eigen::PlainObjectBase<DerivedG> & G,
-  Eigen::PlainObjectBase<DerivedJ> & J,
-  Eigen::PlainObjectBase<Derivedflip> & flip)
-{
-  Eigen::Matrix<typename DerivedV::Scalar,DerivedF::RowsAtCompileTime,3> N;
-  per_face_normals_stable(V,F,N);
-  return outer_hull(V,F,N,G,J,flip);
-}
-
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-#undef IGL_STATIC_LIBRARY
-#include <igl/barycenter.cpp>
-#include <igl/outer_element.cpp>
-#include <igl/cgal/order_facets_around_edges.cpp>
-#define IGL_STATIC_LIBRARY
 template void igl::cgal::outer_hull<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 template void igl::cgal::outer_hull<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template void igl::cgal::outer_hull<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<bool, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<bool, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/cgal/outer_hull.h

@@ -25,26 +25,11 @@ namespace igl
     // Inputs:
     //   V  #V by 3 list of vertex positions
     //   F  #F by 3 list of triangle indices into V
-    //   N  #F by 3 list of per-face normals
     // Outputs:
     //   G  #G by 3 list of output triangle indices into V
     //   J  #G list of indices into F
     //   flip  #F list of whether facet was added to G **and** flipped orientation
     //     (false for faces not added to G)
-    template <
-      typename DerivedV,
-      typename DerivedF,
-      typename DerivedN,
-      typename DerivedG,
-      typename DerivedJ,
-      typename Derivedflip>
-    IGL_INLINE void outer_hull(
-      const Eigen::PlainObjectBase<DerivedV> & V,
-      const Eigen::PlainObjectBase<DerivedF> & F,
-      const Eigen::PlainObjectBase<DerivedN> & N,
-      Eigen::PlainObjectBase<DerivedG> & G,
-      Eigen::PlainObjectBase<DerivedJ> & J,
-      Eigen::PlainObjectBase<Derivedflip> & flip);
     template <
       typename DerivedV,
       typename DerivedF,

include/igl/cgal/peel_outer_hull_layers.cpp

@@ -6,7 +6,6 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "peel_outer_hull_layers.h"
-#include "../per_face_normals.h"
 #include "outer_hull.h"
 #include <vector>
 #include <iostream>
@@ -20,13 +19,11 @@
 template <
   typename DerivedV,
   typename DerivedF,
-  typename DerivedN,
   typename DerivedI,
   typename Derivedflip>
 IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
   const Eigen::PlainObjectBase<DerivedV > & V,
   const Eigen::PlainObjectBase<DerivedF > & F,
-  const Eigen::PlainObjectBase<DerivedN > & N,
   Eigen::PlainObjectBase<DerivedI> & I,
   Eigen::PlainObjectBase<Derivedflip > & flip)
 {
@@ -34,7 +31,6 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
   using namespace std;
   typedef typename DerivedF::Index Index;
   typedef Matrix<typename DerivedF::Scalar,Dynamic,DerivedF::ColsAtCompileTime> MatrixXF;
-  typedef Matrix<typename DerivedN::Scalar,Dynamic,DerivedN::ColsAtCompileTime> MatrixXN;
   typedef Matrix<Index,Dynamic,1> MatrixXI;
   typedef Matrix<typename Derivedflip::Scalar,Dynamic,Derivedflip::ColsAtCompileTime> MatrixXflip;
   const Index m = F.rows();
@@ -51,7 +47,6 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
 #endif
   // keep track of iteration parity and whether flipped in hull
   MatrixXF Fr = F;
-  MatrixXN Nr = N;
   I.resize(m,1);
   flip.resize(m,1);
   // Keep track of index map
@@ -67,15 +62,24 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
     MatrixXI Jo;
     MatrixXflip flipr;
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
-  cout<<"calling outer hull..."<<endl;
-  writePLY(STR("outer-hull-input-"<<iter<<".ply"),V,Fr);
-  writeDMAT(STR("outer-hull-input-"<<iter<<".dmat"),Nr);
+  {
+      cout<<"calling outer hull..." << iter <<endl;
+      std::stringstream ss;
+      ss << "outer_hull_" << iter << ".ply";
+      Eigen::MatrixXd vertices(V.rows(), V.cols());
+      std::transform(V.data(), V.data() + V.rows()*V.cols(),
+              vertices.data(),
+              [](typename DerivedV::Scalar val)
+              {return CGAL::to_double(val); });
+      writePLY(ss.str(), vertices, Fr);
+  }
 #endif
-    outer_hull(V,Fr,Nr,Fo,Jo,flipr);
+    outer_hull(V,Fr,Fo,Jo,flipr);
 #ifdef IGL_PEEL_OUTER_HULL_LAYERS_DEBUG
   writePLY(STR("outer-hull-output-"<<iter<<".ply"),V,Fo);
   cout<<"reindex, flip..."<<endl;
 #endif
+    assert(Fo.rows() != 0);
     assert(Fo.rows() == Jo.rows());
     // all faces in Fo of Fr
     vector<bool> in_outer(Fr.rows(),false);
@@ -89,10 +93,8 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
     // Fr = Fr - Fo
     // update IM
     MatrixXF prev_Fr = Fr;
-    MatrixXN prev_Nr = Nr;
     MatrixXI prev_IM = IM;
     Fr.resize(prev_Fr.rows() - Fo.rows(),F.cols());
-    Nr.resize(Fr.rows(),3);
     IM.resize(Fr.rows());
     {
       Index g = 0;
@@ -101,7 +103,6 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
         if(!in_outer[f])
         {
           Fr.row(g) = prev_Fr.row(f);
-          Nr.row(g) = prev_Nr.row(f);
           IM(g) = prev_IM(f);
           g++;
         }
@@ -112,26 +113,7 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
   return iter;
 }
 
-template <
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedI,
-  typename Derivedflip>
-IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
-  const Eigen::PlainObjectBase<DerivedV > & V,
-  const Eigen::PlainObjectBase<DerivedF > & F,
-  Eigen::PlainObjectBase<DerivedI > & I,
-  Eigen::PlainObjectBase<Derivedflip > & flip)
-{
-  using namespace std;
-  Eigen::Matrix<typename DerivedV::Scalar,DerivedF::RowsAtCompileTime,3> N;
-  per_face_normals(V,F,N);
-  return peel_outer_hull_layers(V,F,N,I,flip);
-}
-
-
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template size_t igl::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
-template size_t igl::cgal::peel_outer_hull_layers<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/cgal/peel_outer_hull_layers.h

@@ -20,24 +20,11 @@ namespace igl
     // Inputs:
     //   V  #V by 3 list of vertex positions
     //   F  #F by 3 list of triangle indices into V
-    //   N  #F by 3 list of per-face normals
     // Outputs:
     //   I  #F list of which peel Iation a facet belongs 
     //   flip  #F list of whether a facet's orientation was flipped when facet
     //     "peeled" into its associated outer hull layer.
     // Returns number of peels
-    template <
-      typename DerivedV,
-      typename DerivedF,
-      typename DerivedN,
-      typename DerivedI,
-      typename Derivedflip>
-    IGL_INLINE size_t peel_outer_hull_layers(
-      const Eigen::PlainObjectBase<DerivedV > & V,
-      const Eigen::PlainObjectBase<DerivedF > & F,
-      const Eigen::PlainObjectBase<DerivedN > & N,
-      Eigen::PlainObjectBase<DerivedI > & I,
-      Eigen::PlainObjectBase<Derivedflip > & flip);
     template <
       typename DerivedV,
       typename DerivedF,

include/igl/file_dialog_open.cpp

@@ -68,9 +68,8 @@ IGL_INLINE std::string igl::file_dialog_open()
       buffer[pos] = (char)ofn.lpstrFile[pos];
       pos++;
     }
-    buffer[pos] = 0;
   } 
-  
+  buffer[pos] = 0;
 #else
   
   // For linux use zenity

include/igl/normal_derivative.cpp

@@ -64,7 +64,7 @@ IGL_INLINE void igl::normal_derivative(
           2);
       DDV *= S;
 
-      IJV.resize(DDV.size());
+      IJV.reserve(DDV.size());
       for(size_t f = 0;f<6*4;f++)
       {
         for(size_t e = 0;e<m;e++)
@@ -95,7 +95,7 @@ IGL_INLINE void igl::normal_derivative(
         slice(C,(VectorXi(12)<<1,1,2,2,2,2,0,0,0,0,1,1).finished(),2);
       DDV *= S;
 
-      IJV.resize(DDV.size());
+      IJV.reserve(DDV.size());
       for(size_t f = 0;f<12;f++)
       {
         for(size_t e = 0;e<m;e++)

tests/CMakeLists.txt

@@ -0,0 +1,28 @@
+# This file is part of libigl, a simple c++ geometry processing library.
+# 
+# Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
+# 
+# This Source Code Form is subject to the terms of the Mozilla Public License 
+# v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+# obtain one at http://mozilla.org/MPL/2.0/.
+#
+# This file is based on PyMesh's unit test setup.
+CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
+
+# Set compiler
+IF(APPLE)
+    SET(CMAKE_C_COMPILER "clang")
+    SET(CMAKE_CXX_COMPILER "clang++")
+ELSEIF(UNIX)
+    SET(CMAKE_C_COMPILER "gcc")
+    SET(CMAKE_CXX_COMPILER "g++")
+ELSEIF(WIN32)
+    # TODO: not how the code would work on Windows.
+ENDIF()
+
+PROJECT(libigl_unit_tests)
+
+INCLUDE(Settings.cmake)
+
+# Process code in each subdirectories
+ADD_SUBDIRECTORY(${PROJECT_SOURCE_DIR}/include/igl)

tests/Settings.cmake

@@ -0,0 +1,54 @@
+# This file is part of libigl, a simple c++ geometry processing library.
+# 
+# Copyright (C) 2015 Qingnan Zhou <qnzhou@gmail.com>
+# 
+# This Source Code Form is subject to the terms of the Mozilla Public License 
+# v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+# obtain one at http://mozilla.org/MPL/2.0/.
+#
+# This file is based on PyMesh's unit test setup.
+
+# Include directories to search for source.
+INCLUDE_DIRECTORIES(${PROJECT_SOURCE_DIR}/src)
+GET_FILENAME_COMPONENT(LIBIGL_PATH ${PROJECT_SOURCE_DIR} DIRECTORY)
+INCLUDE_DIRECTORIES(${LIBIGL_PATH}/include/)
+
+# Set build type.
+SET(CMAKE_BUILD_TYPE Debug)
+#SET(CMAKE_BUILD_TYPE Release)
+
+# Create 64 bits binary.  32 bits support is dropped.
+SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
+SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0")
+SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Os")
+SET(CMAKE_LIBRARY_PATH /opt/local/lib ${CMAKE_LIBRARY_PATH})
+
+# Set output directories
+SET(EXECUTABLE_OUTPUT_PATH ${PROJECT_SOURCE_DIR}/bin)
+MAKE_DIRECTORY(${EXECUTABLE_OUTPUT_PATH})
+
+LINK_DIRECTORIES(/opt/local/lib)
+
+SET(CMAKE_MACOSX_RPATH ON)
+
+# Set TEST_DIR definition
+ADD_DEFINITIONS(-DTEST_DIR="${PROJECT_SOURCE_DIR}")
+
+# Include current directory
+INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR})
+
+# Include Eigen
+#SET(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
+SET(CMAKE_MODULE_PATH ${LIBIGL_PATH}/tutorial/cmake)
+FIND_PACKAGE(Eigen REQUIRED)
+INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS})
+INCLUDE_DIRECTORIES(${EIGEN_INCLUDE_DIRS}/unsupported)
+ADD_DEFINITIONS(-DEIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET)
+
+# Add googletest googlemock support
+ADD_SUBDIRECTORY(${PROJECT_SOURCE_DIR}/googletest/googlemock)
+SET(GTEST_BOTH_LIBRARIES gtest gmock)
+INCLUDE_DIRECTORIES(${gmock_SOURCE_DIR})
+INCLUDE_DIRECTORIES(${gmock_SOURCE_DIR}/include)
+INCLUDE_DIRECTORIES(${gtest_SOURCE_DIR})
+INCLUDE_DIRECTORIES(${gtest_SOURCE_DIR}/include)

tests/data/cube.obj

@@ -0,0 +1,21 @@
+# Generated with PyMesh
+v -0.5 -0.5 -0.5 
+v 0.5 -0.5 -0.5 
+v 0.5 0.5 -0.5 
+v -0.5 0.5 -0.5 
+v -0.5 -0.5 0.5 
+v 0.5 -0.5 0.5 
+v 0.5 0.5 0.5 
+v -0.5 0.5 0.5 
+f 1 3 2 
+f 1 4 3 
+f 6 1 2 
+f 6 5 1 
+f 3 7 2 
+f 4 7 3 
+f 6 2 7 
+f 6 7 5 
+f 1 8 4 
+f 1 5 8 
+f 4 8 7 
+f 7 8 5 

tests/googletest

@@ -0,0 +1 @@
+Subproject commit de411c3e80120f8dcc2a3f4f62f3ca692c0431d7

tests/include/igl/CMakeLists.txt

@@ -0,0 +1,9 @@
+FILE(GLOB TEST_SRC_FILES *.cpp main.cpp)
+FILE(GLOB TEST_INC_FILES *.h *.inl)
+
+ADD_EXECUTABLE(igl_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+TARGET_LINK_LIBRARIES(igl_tests ${GTEST_BOTH_LIBRARIES})
+ADD_CUSTOM_COMMAND(TARGET igl_tests POST_BUILD COMMAND igl_tests)
+
+ADD_SUBDIRECTORY(cgal)
+ADD_SUBDIRECTORY(boolean)

tests/include/igl/boolean/CMakeLists.txt

@@ -0,0 +1,19 @@
+# Check if CGAL is available
+IF (NOT CGAL_FOUND)
+    IF (DEFINED ENV{CGAL_PATH} AND NOT DEFINED ENV{CGAL_DIR})
+        SET(CGAL_DIR $ENV{CGAL_PATH})
+    ENDIF (DEFINED ENV{CGAL_PATH} AND NOT DEFINED ENV{CGAL_DIR})
+    SET(CGAL_DONT_OVERRIDE_CMAKE_FLAGS TRUE CACHE BOOL
+        "Disable CGAL from overwriting my cmake flags")
+    FIND_PACKAGE(CGAL QUIET)
+    INCLUDE(${CGAL_USE_FILE})
+ENDIF (NOT CGAL_FOUND)
+
+IF (CGAL_FOUND)
+    FILE(GLOB TEST_SRC_FILES *.cpp main.cpp)
+    FILE(GLOB TEST_INC_FILES *.h *.inl)
+
+    ADD_EXECUTABLE(igl_boolean_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+    TARGET_LINK_LIBRARIES(igl_boolean_tests ${GTEST_BOTH_LIBRARIES} ${CGAL_LIBRARIES})
+    ADD_CUSTOM_COMMAND(TARGET igl_boolean_tests POST_BUILD COMMAND igl_boolean_tests)
+ENDIF (CGAL_FOUND)

tests/include/igl/boolean/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/boolean/mesh_boolean.cpp

@@ -0,0 +1,26 @@
+#include <test_common.h>
+
+#include <igl/boolean/mesh_boolean.h>
+#include <igl/boolean/MeshBooleanType.h>
+#include <igl/exterior_edges.h>
+
+TEST(MeshBoolean, TwoCubes) {
+    Eigen::MatrixXd V1;
+    Eigen::MatrixXi F1;
+    test_common::load_mesh("two-boxes-bad-self-union.ply", V1, F1);
+
+    Eigen::MatrixXd V2(0, 3);
+    Eigen::MatrixXi F2(0, 3);
+
+    Eigen::MatrixXd Vo;
+    Eigen::MatrixXi Fo;
+
+    igl::boolean::mesh_boolean(V1, F1, V2, F2,
+            igl::boolean::MESH_BOOLEAN_TYPE_UNION,
+            Vo, Fo);
+
+    Eigen::MatrixXi Eb;
+    igl::exterior_edges(Fo, Eb);
+
+    ASSERT_EQ(0, Eb.rows());
+}

tests/include/igl/cgal/CMakeLists.txt

@@ -0,0 +1,17 @@
+# Check if CGAL is available
+IF (DEFINED ENV{CGAL_PATH} AND NOT DEFINED ENV{CGAL_DIR})
+    SET(CGAL_DIR $ENV{CGAL_PATH})
+ENDIF (DEFINED ENV{CGAL_PATH} AND NOT DEFINED ENV{CGAL_DIR})
+SET(CGAL_DONT_OVERRIDE_CMAKE_FLAGS TRUE CACHE BOOL
+    "Disable CGAL from overwriting my cmake flags")
+FIND_PACKAGE(CGAL QUIET)
+INCLUDE(${CGAL_USE_FILE})
+
+IF (CGAL_FOUND)
+    FILE(GLOB TEST_SRC_FILES *.cpp main.cpp)
+    FILE(GLOB TEST_INC_FILES *.h *.inl)
+
+    ADD_EXECUTABLE(igl_cgal_tests ${TEST_SRC_FILES} ${TEST_INC_FILES})
+    TARGET_LINK_LIBRARIES(igl_cgal_tests ${GTEST_BOTH_LIBRARIES} ${CGAL_LIBRARIES})
+    ADD_CUSTOM_COMMAND(TARGET igl_cgal_tests POST_BUILD COMMAND igl_cgal_tests)
+ENDIF (CGAL_FOUND)

tests/include/igl/cgal/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/cgal/order_facets_around_edges.cpp

@@ -0,0 +1,187 @@
+#include <test_common.h>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <igl/cgal/order_facets_around_edges.h>
+#include <igl/unique_edge_map.h>
+#include <igl/readDMAT.h>
+#include <igl/per_face_normals.h>
+
+namespace order_facets_around_edges_test {
+
+typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+
+template<typename T>
+size_t index_of(const std::vector<T>& array, T val) {
+    auto loc = std::find(array.begin(), array.end(), val);
+    assert(loc != array.end());
+    return loc - array.begin();
+}
+
+void assert_consistently_oriented(size_t num_faces,
+        const std::vector<int>& expected_face_order,
+        const std::vector<int>& e_order) {
+    const size_t num_items = expected_face_order.size();
+    ASSERT_EQ(num_items, e_order.size());
+
+    std::vector<int> order(num_items);
+    std::transform(e_order.begin(), e_order.end(), order.begin(),
+            [=](int val) { return val % num_faces; });
+
+    size_t ref_start = index_of(order, expected_face_order[0]);
+    for (size_t i=0; i<num_items; i++) {
+        ASSERT_EQ(expected_face_order[i], order[(ref_start + i) % num_items]);
+    }
+}
+
+template<typename DerivedV, typename DerivedF>
+void assert_order(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        size_t v0, size_t v1,
+        std::vector<int> expected_order, const std::string& normal="") {
+    Eigen::MatrixXi E, uE, EMAP;
+    std::vector<std::vector<int> > uE2E;
+    igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+
+    std::vector<std::vector<int> > uE2oE;
+    std::vector<std::vector<bool> > uE2C;
+
+    if (normal != "") {
+        Eigen::MatrixXd N;
+        //igl::per_face_normals_stable(V, F, N);
+        //igl::per_face_normals(V, F, N);
+        test_common::load_matrix(normal, N);
+        igl::cgal::order_facets_around_edges(V, F, N, E, uE, EMAP, uE2E, uE2oE, uE2C);
+    } else {
+        igl::cgal::order_facets_around_edges(V, F, E, uE, EMAP, uE2E, uE2oE, uE2C);
+    }
+
+    const size_t num_uE = uE.rows();
+    for (size_t i=0; i<num_uE; i++) {
+        const auto& order = uE2oE[i];
+        const auto& cons  = uE2C[i];
+        Eigen::VectorXi e = uE.row(i);
+        if (order.size() <= 1) continue;
+        if (e[0] != v0 && e[0] != v1) continue;
+        if (e[1] != v0 && e[1] != v1) continue;
+        if (e[0] == v1 && e[1] == v0) {
+            std::reverse(expected_order.begin(), expected_order.end());
+        }
+        assert_consistently_oriented(F.rows(), expected_order, order);
+    }
+}
+
+TEST(OrderFacetsAroundEdges, Simple) {
+    Eigen::MatrixXd V(4, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0;
+    Eigen::MatrixXi F(2, 3);
+    F << 0, 1, 2,
+         2, 1, 3;
+
+    assert_order(V, F, 1, 2, {0, 1});
+}
+
+TEST(OrderFacetsAroundEdges, TripletFaces) {
+    Eigen::MatrixXd V(5, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::MatrixXi F(3, 3);
+    F << 0, 1, 2,
+         2, 1, 3,
+         1, 2, 4;
+
+    assert_order(V, F, 1, 2, {0, 1, 2});
+}
+
+TEST(OrderFacetsAroundEdges, DuplicatedFaces) {
+    Eigen::MatrixXd V(5, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::MatrixXi F(4, 3);
+    F << 0, 1, 2,
+         2, 1, 3,
+         1, 2, 4,
+         4, 1, 2;
+
+    assert_order(V, F, 1, 2, {0, 1, 3, 2});
+}
+
+TEST(OrderFacetsAroundEdges, MultipleDuplicatedFaces) {
+    Eigen::MatrixXd V(5, 3);
+    V << 0.0, 0.0, 0.0,
+         1.0, 0.0, 0.0,
+         0.0, 1.0, 0.0,
+         1.0, 1.0, 0.0,
+         0.0, 0.0, 1.0;
+    Eigen::MatrixXi F(6, 3);
+    F << 0, 1, 2,
+         1, 2, 0,
+         2, 1, 3,
+         1, 3, 2,
+         1, 2, 4,
+         4, 1, 2;
+
+    assert_order(V, F, 1, 2, {1, 0, 2, 3, 5, 4});
+}
+
+TEST(OrderFacetsAroundEdges, Debug) {
+    Eigen::MatrixXd V(5, 3);
+    V <<
+        -44.3205080756887781, 4.22994972382184579e-15, 75,
+        -27.933756729740665, -48.382685902179837, 75,
+        -55.8675134594812945, -2.81996648254789745e-15, 75,
+        -27.933756729740665, -48.382685902179837, 70,
+        -31.4903810567666049, -42.2224318643354408, 85;
+
+    Eigen::MatrixXi F(3, 3);
+    F << 1, 0, 2,
+         2, 3, 1,
+         4, 1, 2;
+
+    assert_order(V, F, 1, 2, {0, 2, 1});
+}
+
+TEST(OrderFacetsAroundEdges, Debug2) {
+    Eigen::MatrixXd V(5, 3);
+    V <<
+        -22.160254037844382, 38.3826859021798441, 75,
+        -27.9337567297406331, 48.3826859021798654, 75,
+        27.9337567297406544, 48.3826859021798512, 75,
+        27.9337567297406544, 48.3826859021798512, 70,
+        20.8205080756887924, 48.3826859021798512, 85;
+    Eigen::MatrixXi F(3, 3);
+    F << 1, 0, 2,
+         3, 1, 2,
+         2, 4, 1;
+
+    assert_order(V, F, 1, 2, {1, 0, 2});
+}
+
+TEST(OrderFacetsAroundEdges, NormalSensitivity) {
+    // This example shows that epsilon difference in normal vectors could
+    // results in very different ordering of facets.
+
+    Eigen::MatrixXd V;
+    test_common::load_matrix("duplicated_faces_V.dmat", V);
+    Eigen::MatrixXi F;
+    test_common::load_matrix("duplicated_faces_F.dmat", F);
+
+    assert_order(V, F, 223, 224, {2, 0, 3, 1}, "duplicated_faces_N1.dmat");
+    assert_order(V, F, 223, 224, {0, 3, 2, 1}, "duplicated_faces_N2.dmat");
+}
+
+
+}

tests/include/igl/cgal/peel_outer_hull_layers.cpp

@@ -0,0 +1,51 @@
+#include <test_common.h>
+#include <iostream>
+#include <Eigen/Dense>
+
+#include <igl/cgal/peel_outer_hull_layers.h>
+#include <igl/cgal/remesh_self_intersections.h>
+#include <igl/cgal/RemeshSelfIntersectionsParam.h>
+#include <igl/per_face_normals.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/writeOBJ.h>
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+TEST(PeelOuterHullLayers, TwoCubes) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("two-boxes-bad-self-union.ply", V, F);
+    ASSERT_EQ(486, V.rows());
+    ASSERT_EQ(708, F.rows());
+
+    typedef CGAL::Exact_predicates_exact_constructions_kernel K;
+    typedef K::FT Scalar;
+    typedef Eigen::Matrix<Scalar,
+            Eigen::Dynamic,
+            Eigen::Dynamic> MatrixXe;
+
+    MatrixXe Vs;
+    Eigen::MatrixXi Fs, IF;
+    Eigen::VectorXi J, IM;
+    igl::cgal::RemeshSelfIntersectionsParam param;
+    igl::cgal::remesh_self_intersections(V, F, param, Vs, Fs, IF, J, IM);
+
+    std::for_each(Fs.data(),Fs.data()+Fs.size(),
+            [&IM](int & a){ a=IM(a); });
+    MatrixXe Vt;
+    Eigen::MatrixXi Ft;
+    igl::remove_unreferenced(Vs,Fs,Vt,Ft,IM);
+    const size_t num_faces = Ft.rows();
+
+    Eigen::VectorXi I, flipped;
+    size_t num_peels = igl::cgal::peel_outer_hull_layers(Vt, Ft, I, flipped);
+
+    Eigen::MatrixXd vertices(Vt.rows(), Vt.cols());
+    std::transform(Vt.data(), Vt.data() + Vt.rows() * Vt.cols(),
+            vertices.data(), [](Scalar v) { return CGAL::to_double(v); });
+    igl::writeOBJ("debug.obj", vertices, Ft);
+
+    ASSERT_EQ(num_faces, I.rows());
+    ASSERT_EQ(0, I.minCoeff());
+    ASSERT_EQ(1, I.maxCoeff());
+}

tests/include/igl/main.cpp

@@ -0,0 +1,6 @@
+#include <gtest/gtest.h>
+
+int main(int argc, char **argv) {
+      ::testing::InitGoogleTest(&argc, argv);
+        return RUN_ALL_TESTS();
+}

tests/include/igl/readDMAT.cpp

@@ -0,0 +1,19 @@
+#include <test_common.h>
+
+TEST(readDMAT, Comp) {
+    Eigen::MatrixXd N1, N2;
+    test_common::load_matrix("duplicated_faces_N1.dmat", N1);
+    test_common::load_matrix("duplicated_faces_N2.dmat", N2);
+
+    ASSERT_EQ(N1.rows(), N2.rows());
+    ASSERT_EQ(N1.cols(), N2.cols());
+    ASSERT_FALSE(((N1-N2).array() != 0.0).all());
+
+    const size_t rows = N1.rows();
+    const size_t cols = N1.cols();
+    for (size_t i=0; i<rows; i++) {
+        for (size_t j=0; j<cols; j++) {
+            ASSERT_FLOAT_EQ(N1(i,j), N2(i,j));
+        }
+    }
+}

tests/include/igl/readOBJ.cpp

@@ -0,0 +1,9 @@
+#include <test_common.h>
+
+TEST(readOBJ, simple) {
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    test_common::load_mesh("cube.obj", V, F);
+    ASSERT_EQ(8, V.rows());
+    ASSERT_EQ(12, F.rows());
+}

tests/test_common.h

@@ -0,0 +1,27 @@
+#pragma once
+#include <string>
+
+#include <Eigen/Core>
+#include <gtest/gtest.h>
+
+#include <igl/read_triangle_mesh.h>
+#include <igl/readDMAT.h>
+
+namespace test_common {
+    template<typename DerivedV, typename DerivedF>
+    void load_mesh(
+            const std::string& filename, 
+            Eigen::PlainObjectBase<DerivedV>& V,
+            Eigen::PlainObjectBase<DerivedF>& F) {
+        auto find_file = [&](const std::string& val) {
+            return std::string(TEST_DIR) + "/data/" + val;
+        };
+        igl::read_triangle_mesh(find_file(filename), V, F);
+    }
+
+    template<typename Derived>
+    void load_matrix(const std::string& filename,
+            Eigen::PlainObjectBase<Derived>& M) {
+        igl::readDMAT(std::string(TEST_DIR) + "/data/" + filename, M);
+    }
+}

tutorial/cmake/FindGLEW.cmake

@@ -32,7 +32,7 @@ if(GLEW_FOUND)
   set(GLEW_SOURCES ${GLEW_INCLUDE_DIR}/../src/glew.c)
   message(STATUS "Found GLEW: ${GLEW_INCLUDE_DIR}")
 else(GLEW_FOUND)
-  message(FATAL_ERROR "could NOT find glew")
+  message(WARNING "could NOT find glew")
 endif(GLEW_FOUND)
 
 MARK_AS_ADVANCED(GLEW_INCLUDE_DIR)

tutorial/cmake/FindGLFWH.cmake

@@ -37,7 +37,7 @@ if(GLFW_FOUND)
   message(STATUS "Found GLFW: ${GLFW_INCLUDE_DIR} -- HEADERS ONLY")
 else(GLFW_FOUND)
   if (NOT GLFW_FIND_QUIETLY)
-    message(FATAL_ERROR "could NOT find GLFW")
+    message(WARNING "could NOT find GLFW")
   endif (NOT GLFW_FIND_QUIETLY)
 endif(GLFW_FOUND)