Remove normal dependencies from outer hull and mesh boolean.

Former-commit-id: 8594dde52a16cae5b5576dd1da6764d5783377b8

include/igl/boolean/mesh_boolean.cpp

@@ -187,13 +187,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 +195,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;
 

include/igl/cgal/order_facets_around_edge.cpp

@@ -0,0 +1,210 @@
+#include "order_facets_around_edge.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);
+        return;
+    } else if (adj_faces.size() == 1) {
+        order.resize(1);
+        order[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);
+        switch (CGAL::orientation(ps, pd, p1, p2)) {
+            case CGAL::POSITIVE:
+                order[0] = 1;
+                order[1] = 0;
+                break;
+            case CGAL::NEGATIVE:
+                order[0] = 0;
+                order[1] = 1;
+                break;
+            case CGAL::COPLANAR:
+                order[0] = adj_faces[0] < adj_faces[1] ? 0:1;
+                order[1] = 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);
+
+    size_t count=0;
+    for (size_t i=0; i<tie_positive_size; i++) {
+        order[count+i] =
+            tie_positive_oriented_index[tie_positive_order[i]];
+    }
+    count += tie_positive_size;
+
+    for (size_t i=0; i<negative_size; i++) {
+        order[count+i] = negative_side_index[negative_order[i]];
+    }
+    count += negative_size;
+
+    for (size_t i=0; i<tie_negative_size; i++) {
+        order[count+i] =
+            tie_negative_oriented_index[tie_negative_order[i]];
+    }
+    count += tie_negative_size;
+
+    for (size_t i=0; i<positive_size; i++) {
+        order[count+i] = positive_side_index[positive_order[i]];
+    }
+    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]];
+        const Point_3& p_b =
+            opposite_vertices[order[(i+1)%num_adj_faces]];
+        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] = circular_order[(start_idx + i)%num_adj_faces];
+    }
+}

include/igl/cgal/order_facets_around_edge.h

@@ -0,0 +1,51 @@
+// 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"
+
+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>
@@ -256,215 +257,6 @@ igl::cgal::order_facets_around_edges(
     }
 }
 
-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 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) {
-
-                // 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);
-                    return;
-                } else if (adj_faces.size() == 1) {
-                    order.resize(1);
-                    order[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);
-                    switch (CGAL::orientation(ps, pd, p1, p2)) {
-                        case CGAL::POSITIVE:
-                            order[0] = 1;
-                            order[1] = 0;
-                            break;
-                        case CGAL::NEGATIVE:
-                            order[0] = 0;
-                            order[1] = 1;
-                            break;
-                        case CGAL::COPLANAR:
-                            order[0] = adj_faces[0] < adj_faces[1] ? 0:1;
-                            order[1] = 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;
-                igl::cgal::order_facets_around_edges_helper::order_facets_around_edge(
-                        V, F, s, d, positive_side, positive_order);
-                igl::cgal::order_facets_around_edges_helper::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);
-
-                size_t count=0;
-                for (size_t i=0; i<tie_positive_size; i++) {
-                    order[count+i] =
-                        tie_positive_oriented_index[tie_positive_order[i]];
-                }
-                count += tie_positive_size;
-
-                for (size_t i=0; i<negative_size; i++) {
-                    order[count+i] = negative_side_index[negative_order[i]];
-                }
-                count += negative_size;
-
-                for (size_t i=0; i<tie_negative_size; i++) {
-                    order[count+i] =
-                        tie_negative_oriented_index[tie_negative_order[i]];
-                }
-                count += tie_negative_size;
-
-                for (size_t i=0; i<positive_size; i++) {
-                    order[count+i] = positive_side_index[positive_order[i]];
-                }
-                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]];
-                    const Point_3& p_b =
-                        opposite_vertices[order[(i+1)%num_adj_faces]];
-                    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] = circular_order[(start_idx + i)%num_adj_faces];
-                }
-            }
-        }
-    }
-}
-
 template<
     typename DerivedV,
     typename DerivedF,
@@ -526,8 +318,7 @@ IGL_INLINE void igl::cgal::order_facets_around_edges(
         }
 
         Eigen::VectorXi order;
-        igl::cgal::order_facets_around_edges_helper::order_facets_around_edge(
-                V, F, s, d, adj_faces, order);
+        order_facets_around_edge(V, F, s, d, adj_faces, order);
         assert(order.size() == edge_valance);
 
         auto& ordered_edges = uE2oE[ui];

include/igl/cgal/outer_facet.cpp

@@ -0,0 +1,68 @@
+// 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:
+    //    Find an outer edge.
+    //    Order adjacent facets around the edge.
+    //    For consecutive facets, check if the tetrahedron formed by them is
+    //    positively oriented.
+    //    If yes, check the next pair.
+    //    If not, return either face.
+    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"
@@ -31,14 +31,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)
@@ -54,7 +52,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:
@@ -95,7 +92,6 @@ 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);
@@ -162,7 +158,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;
@@ -225,19 +221,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);
@@ -247,11 +230,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)<<", |"<<
@@ -283,64 +261,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];
@@ -515,24 +435,6 @@ 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
@@ -543,5 +445,4 @@ IGL_INLINE void igl::cgal::outer_hull(
 #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
@@ -69,9 +64,8 @@ IGL_INLINE size_t igl::cgal::peel_outer_hull_layers(
 #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);
 #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;
@@ -89,10 +83,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 +93,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 +103,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,