Merge branch 'qnzhou-master'

Former-commit-id: 1a538759cfebac47c4fd27b89b6852cbcb749d5a

include/igl/boolean/mesh_boolean.cpp

@@ -6,6 +6,8 @@
 #include <igl/unique_simplices.h>
 #include <iostream>
 
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
 //#define IGL_MESH_BOOLEAN_DEBUG
 
 template <
@@ -93,9 +95,12 @@ IGL_INLINE void igl::mesh_boolean(
   using namespace igl;
   MeshBooleanType eff_type = type;
   // Concatenate A and B into a single mesh
+  typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+  typedef Kernel::FT ExactScalar;
   typedef typename DerivedVC::Scalar Scalar;
   typedef typename DerivedFC::Scalar Index;
   typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
+  typedef Matrix<ExactScalar,Dynamic,3> MatrixX3ES;
   typedef Matrix<Index,Dynamic,3> MatrixX3I;
   typedef Matrix<Index,Dynamic,2> MatrixX2I;
   typedef Matrix<Index,Dynamic,1> VectorXI;
@@ -131,11 +136,11 @@ IGL_INLINE void igl::mesh_boolean(
   const auto & libigl_resolve = [](
     const MatrixX3S & V,
     const MatrixX3I & F,
-    MatrixX3S & CV,
+    MatrixX3ES & CV,
     MatrixX3I & CF,
     VectorXJ & J)
   {
-    MatrixX3S SV;
+    MatrixX3ES SV;
     MatrixX3I SF;
     MatrixX2I SIF;
     VectorXI SIM,UIM;
@@ -151,6 +156,7 @@ IGL_INLINE void igl::mesh_boolean(
   cout<<"resolve..."<<endl;
 #endif
   MatrixX3S CV;
+  MatrixX3ES EV;
   MatrixX3I CF;
   VectorXJ CJ;
   if(resolve_fun)
@@ -158,7 +164,12 @@ IGL_INLINE void igl::mesh_boolean(
     resolve_fun(V,F,CV,CF,CJ);
   }else
   {
-    libigl_resolve(V,F,CV,CF,CJ);
+    libigl_resolve(V,F,EV,CF,CJ);
+    CV.resize(EV.rows(), EV.cols());
+    std::transform(EV.data(), EV.data() + EV.rows()*EV.cols(),
+            CV.data(), [&](ExactScalar val) {
+            return CGAL::to_double(val);
+            });
   }
 
   if(type == MESH_BOOLEAN_TYPE_RESOLVE)
@@ -183,7 +194,7 @@ IGL_INLINE void igl::mesh_boolean(
   // peel layers keeping track of odd and even flips
   Matrix<bool,Dynamic,1> odd;
   Matrix<bool,Dynamic,1> flip;
-  peel_outer_hull_layers(CV,CF,CN,odd,flip);
+  peel_outer_hull_layers(EV,CF,CN,odd,flip);
 
 #ifdef IGL_MESH_BOOLEAN_DEBUG
   cout<<"categorize..."<<endl;
@@ -234,6 +245,7 @@ IGL_INLINE void igl::mesh_boolean(
   cout<<"clean..."<<endl;
 #endif
   // Deal with duplicate faces
+  if (false)
   {
     VectorXi IA,IC;
     MatrixX3I uG;

include/igl/order_facets_around_edges.cpp

@@ -0,0 +1,246 @@
+#include "order_facets_around_edges.h"
+#include <type_traits>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename uE2oEType,
+    typename uE2CType >
+IGL_INLINE
+typename std::enable_if<!std::is_same<typename DerivedV::Scalar,
+typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+igl::order_facets_around_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedN>& N,
+        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 DerivedN::Scalar, 3, 1> Vector3F;
+    const typename DerivedV::Scalar EPS = 1e-12;
+    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);
+        Vector3F ref_normal = N.row(ref_face);
+
+        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);
+
+        Vector3F edge = V.row(d) - V.row(s);
+        auto edge_len = edge.norm();
+        bool degenerated = edge_len < EPS;
+        if (degenerated) {
+            if (edge_valance <= 2) {
+                // There is only one way to order 2 or less faces.
+                edge.setZero();
+            } else {
+                edge.setZero();
+                Eigen::Matrix<typename DerivedN::Scalar, Eigen::Dynamic, 3>
+                    normals(edge_valance, 3);
+                for (size_t fei=0; fei<edge_valance; fei++) {
+                    const auto fe = adj_edges[fei];
+                    const auto f = edge_index_to_face_index(fe);
+                    normals.row(fei) = N.row(f);
+                }
+                for (size_t i=0; i<edge_valance; i++) {
+                    size_t j = (i+1) % edge_valance;
+                    Vector3F ni = normals.row(i);
+                    Vector3F nj = normals.row(j);
+                    edge = ni.cross(nj);
+                    edge_len = edge.norm();
+                    if (edge_len >= EPS) {
+                        edge.normalize();
+                        break;
+                    }
+                }
+
+                // Ensure edge direction are consistent with reference face.
+                Vector3F in_face_vec = V.row(o) - V.row(s);
+                if (edge.cross(in_face_vec).dot(ref_normal) < 0) {
+                    edge *= -1;
+                }
+
+                if (edge.norm() < EPS) {
+                    std::cerr << "=====================================" << std::endl;
+                    std::cerr << "  ui: " << ui << std::endl;
+                    std::cerr << "edge: " << ref_edge << std::endl;
+                    std::cerr << "face: " << ref_face << std::endl;
+                    std::cerr << "  vs: " << V.row(s) << std::endl;
+                    std::cerr << "  vd: " << V.row(d) << std::endl;
+                    std::cerr << "adj face normals: " << std::endl;
+                    std::cerr << normals << std::endl;
+                    std::cerr << "Very degenerated case detected:" << std::endl;
+                    std::cerr << "Near zero edge surrounded by "
+                        << edge_valance << " neearly colinear faces" <<
+                        std::endl;
+                    std::cerr << "=====================================" << std::endl;
+                }
+            }
+        } else {
+            edge.normalize();
+        }
+
+        Eigen::MatrixXd angle_data(edge_valance, 3);
+        std::vector<bool> cons(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));
+            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)));
+            Vector3F n = N.row(f);
+            angle_data(fei, 0) = ref_normal.cross(n).dot(edge);
+            angle_data(fei, 1) = ref_normal.dot(n);
+            if (cons[fei]) {
+                angle_data(fei, 0) *= -1;
+                angle_data(fei, 1) *= -1;
+            }
+            angle_data(fei, 0) *= -1; // Sort clockwise.
+            angle_data(fei, 2) = (cons[fei]?1.:-1.)*(f+1);
+        }
+
+        Eigen::VectorXi order;
+        igl::sort_angles(angle_data, order);
+
+        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]];
+        }
+    }
+}
+
+template<
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedN,
+    typename DerivedE,
+    typename DeriveduE,
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename uE2oEType,
+    typename uE2CType >
+IGL_INLINE 
+typename std::enable_if<std::is_same<typename DerivedV::Scalar,
+typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+igl::order_facets_around_edges(
+        const Eigen::PlainObjectBase<DerivedV>& V,
+        const Eigen::PlainObjectBase<DerivedF>& F,
+        const Eigen::PlainObjectBase<DerivedN>& N,
+        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 DerivedN::Scalar, 3, 1> Vector3F;
+    typedef Eigen::Matrix<typename DerivedV::Scalar, 3, 1> Vector3E;
+    const typename DerivedV::Scalar EPS = 1e-12;
+    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);
+        Vector3F ref_normal = N.row(ref_face);
+
+        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);
+
+        Vector3E exact_edge = V.row(d) - V.row(s);
+        exact_edge.array() /= exact_edge.squaredNorm();
+        Vector3F edge(
+                CGAL::to_double(exact_edge[0]),
+                CGAL::to_double(exact_edge[1]),
+                CGAL::to_double(exact_edge[2]));
+        edge.normalize();
+
+        Eigen::MatrixXd angle_data(edge_valance, 3);
+        std::vector<bool> cons(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));
+            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)));
+            Vector3F n = N.row(f);
+            angle_data(fei, 0) = ref_normal.cross(n).dot(edge);
+            angle_data(fei, 1) = ref_normal.dot(n);
+            if (cons[fei]) {
+                angle_data(fei, 0) *= -1;
+                angle_data(fei, 1) *= -1;
+            }
+            angle_data(fei, 0) *= -1; // Sort clockwise.
+            angle_data(fei, 2) = (cons[fei]?1.:-1.)*(f+1);
+        }
+
+        Eigen::VectorXi order;
+        igl::sort_angles(angle_data, order);
+
+        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/order_facets_around_edges.h

@@ -0,0 +1,79 @@
+#ifndef ORDER_FACETS_AROUND_EDGES
+#define ORDER_FACETS_AROUND_EDGES
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+namespace igl {
+    // For each undirected edge, sort its adjacent faces.
+    //
+    // Inputs:
+    //   V    #V by 3 list of vertices.
+    //   F    #F by 3 list of faces
+    //   N    #F by 3 list of face normals.
+    //   E    #F*3 by 2 list vertex indices, represents directed edges.
+    //  uE    #uE by 2 list of vertex_indices, represents undirected edges.
+    //  EMAP  #F*3 list of indices that maps E to uE. (a many-to-one map)
+    //  uE2E  #uE list of lists that maps uE to E. (a one-to-many map)
+    //
+    // Outputs:
+    //   uE2oE #uE list of lists that maps uE to an ordered list of E. (a
+    //         one-to-many map)
+    //   uE2C  #uE list of lists of bools indicates whether each face in
+    //         uE2oE[i] is consistently orientated as the ordering.
+    //
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedN,
+        typename DerivedE,
+        typename DeriveduE,
+        typename DerivedEMAP,
+        typename uE2EType,
+        typename uE2oEType,
+        typename uE2CType >
+    IGL_INLINE
+    typename std::enable_if<!std::is_same<typename DerivedV::Scalar,
+    typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+    order_facets_around_edges(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DerivedN>& N,
+            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 );
+
+    template<
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedN,
+        typename DerivedE,
+        typename DeriveduE,
+        typename DerivedEMAP,
+        typename uE2EType,
+        typename uE2oEType,
+        typename uE2CType >
+    IGL_INLINE 
+    typename std::enable_if<std::is_same<typename DerivedV::Scalar,
+    typename CGAL::Exact_predicates_exact_constructions_kernel::FT>::value, void>::type
+    order_facets_around_edges(
+            const Eigen::PlainObjectBase<DerivedV>& V,
+            const Eigen::PlainObjectBase<DerivedF>& F,
+            const Eigen::PlainObjectBase<DerivedN>& N,
+            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 );
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "order_facets_around_edges.cpp"
+#endif
+
+#endif

include/igl/outer_facet.cpp

@@ -30,6 +30,10 @@ IGL_INLINE void igl::outer_facet(
   assert(V.cols() == 3);
   assert(N.cols() == 3);
   Index max_v = -1;
+  auto generic_fabs = [&](Scalar val) { 
+      if (val >= 0) return val;
+      else return -val;
+  };
   for(size_t d = 0;d<(size_t)V.cols();d++)
   {
     Scalar max_d = -1e26;
@@ -38,19 +42,37 @@ IGL_INLINE void igl::outer_facet(
     {
       const Index f = I(i);
       const Scalar nd = N(f,d);
-      if(fabs(nd)>0)
+      if(generic_fabs(nd)>0)
       {
         for(Index c = 0;c<3;c++)
         {
           const Index v = F(f,c);
           if(v == max_v)
           {
-            if(fabs(nd) > max_nd)
+            if(generic_fabs(nd) > max_nd)
             {
               // Just update max face and normal
               max_f = f;
-              max_nd = fabs(nd);
+              max_nd = generic_fabs(nd);
               flip = nd<0;
+            } else if (generic_fabs(nd) == max_nd) {
+                if (nd == max_nd) {
+                    if (flip) {
+                        max_f = f;
+                        max_nd = nd;
+                        flip = false;
+                    } else if (f > max_f){
+                        max_f = f;
+                        max_nd = nd;
+                        flip = false;
+                    }
+                } else {
+                    if (flip && f < max_f) {
+                        max_f = f;
+                        max_nd = generic_fabs(nd);
+                        flip = true;
+                    }
+                }
             }
           }else
           {
@@ -61,7 +83,7 @@ IGL_INLINE void igl::outer_facet(
               max_v = v;
               max_d = vd;
               max_f = f;
-              max_nd = fabs(nd);
+              max_nd = generic_fabs(nd);
               flip = nd<0;
             }
           }

include/igl/outer_hull.cpp

@@ -7,12 +7,16 @@
 #include "unique_edge_map.h"
 #include "barycenter.h"
 #include "per_face_normals.h"
+#include "writePLY.h"
+#include "sort_angles.h"
+#include "order_facets_around_edges.h"
 
 #include <Eigen/Geometry>
 #include <vector>
 #include <map>
 #include <queue>
 #include <iostream>
+#include <CGAL/number_utils.h>
 //#define IGL_OUTER_HULL_DEBUG
 
 template <
@@ -30,6 +34,9 @@ IGL_INLINE void igl::outer_hull(
   Eigen::PlainObjectBase<DerivedJ> & J,
   Eigen::PlainObjectBase<Derivedflip> & flip)
 {
+#ifdef IGL_OUTER_HULL_DEBUG
+  std::cerr << "Extracting outer hull" << std::endl;
+#endif
   using namespace Eigen;
   using namespace std;
   typedef typename DerivedF::Index Index;
@@ -61,118 +68,31 @@ IGL_INLINE void igl::outer_hull(
 #endif
   typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
   typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
+  typedef Matrix<typename DerivedV::Scalar, 3, 1> Vector3F;
   MatrixX2I E,uE;
   VectorXI EMAP;
   vector<vector<typename DerivedF::Index> > uE2E;
   unique_edge_map(F,E,uE,EMAP,uE2E);
-
-  // TODO:
-  // uE --> face-edge index, sorted CCW around edge according to normal
-  // uE --> sorted order index 
-  // uE --> bool, whether needed to flip face to make "consistent" with unique
-  //   edge
-  // Place order of each half-edge in its corresponding sorted list around edge
-  VectorXI diIM(3*m);
-  // Whether face's edge used for sorting is consistent with unique edge
-  VectorXI dicons(3*m);
-  // dihedral angles of faces around edge with face of edge in dicons
-  vector<vector<typename DerivedV::Scalar> > di(uE2E.size());
-  // For each list of face-edges incide on a unique edge
-  for(size_t ui = 0;ui<(size_t)uE.rows();ui++)
-  {
-    // Base normal vector to orient against
-    const auto fe0 = uE2E[ui][0];
-    const RowVector3N & eVp = N.row(fe0%m);
-    MatrixXd di_I(uE2E[ui].size(),2);
-
-    const typename DerivedF::Scalar d = F(fe0%m,((fe0/m)+2)%3);
-    const typename DerivedF::Scalar s = F(fe0%m,((fe0/m)+1)%3);
-    // Edge vector
-    const auto & eV = (V.row(d)-V.row(s)).normalized();
-
-    vector<bool> cons(uE2E[ui].size());
-    // Loop over incident face edges
-    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
-    {
-      const auto & fe = uE2E[ui][fei];
-      const auto f = fe % m;
-      const auto c = fe / m;
-      // source should match destination to be consistent
-      cons[fei] = (d == F(f,(c+1)%3));
-      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)));
-      // Angle between n and f
-      const RowVector3N & n = N.row(f);
-      di_I(fei,0) = M_PI - atan2( eVp.cross(n).dot(eV), eVp.dot(n));
 #ifdef IGL_OUTER_HULL_DEBUG
-      if(di_I(fei,0) != di_I(fei,0) )
-      {
-        cout<<"NaN from face: "<<(f+1)<<endl;
-        cout<<"  n: "<<n<<endl;
-        cout<<"  eVp: "<<eVp<<endl;
-        cout<<"  eV: "<<eV<<endl;
-        cout<<"  eVp x n . eV: "<<(eVp.cross(n).dot(eV))<<endl;
-        cout<<"  eVp . n: "<<(eVp.dot(n))<<endl;
+  for (size_t ui=0; ui<uE.rows(); ui++) {
+      std::cout << ui << ": " << uE2E[ui].size() << " -- (";
+      for (size_t i=0; i<uE2E[ui].size(); i++) {
+          std::cout << uE2E[ui][i] << ", ";
       }
+      std::cout << ")" << std::endl;
+  }
 #endif
-      assert(di_I(fei,0) == di_I(fei,0) && "NaN Alert!");
-      if(!cons[fei])
-      {
-        di_I(fei,0) = di_I(fei,0) + M_PI;
-        if(di_I(fei,0)>=2.*M_PI)
-        {
-          di_I(fei,0) = di_I(fei,0) - 2.*M_PI;
-        }
-      }
-      // This signing is very important to make sure different edges sort
-      // duplicate faces the same way, regardless of their orientations
-      di_I(fei,1) = (cons[fei]?1.:-1.)*f;
-    }
 
-    // Despite the effort to get stable normals the atan2 up doesn't
-    // compute (exactly) -θ for -n if it computes θ for n. So just
-    // explicitly check if there's a duplicate face
-    // Shitty O(val^2) implementation
-    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
-    {
-      const auto & fe = uE2E[ui][fei];
-      const auto f = fe % m;
-      for(size_t gei = fei+1;gei<uE2E[ui].size();gei++)
-      {
-        const auto & ge = uE2E[ui][gei];
-        const auto g = ge % m;
-        if(duplicate_simplex(f,g))
-        {
-#ifdef IGL_OUTER_HULL_DEBUG
-          cout<<"Forcing duplicate: "<<(f+1)<<","<<(g+1)<<endl;
-#endif
-          di_I(gei,0) = di_I(fei,0);
-        }
+  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);
+  uE2E = uE2oE;
+  VectorXI diIM(3*m);
+  for (auto ue : uE2E) {
+      for (size_t i=0; i<ue.size(); i++) {
+          auto fe = ue[i];
+          diIM[fe] = i;
       }
-    }
-    VectorXi IM;
-    //igl::sort(di[ui],true,di[ui],IM);
-    // Sort, but break ties using "signed index" to ensure that duplicates
-    // always show up in same order.
-    MatrixXd s_di_I;
-    igl::sortrows(di_I,true,s_di_I,IM);
-    di[ui].resize(uE2E[ui].size());
-    for(size_t i = 0;i<di[ui].size();i++)
-    {
-      di[ui][i] = s_di_I(i,0);
-    }
-
-    // copy old list
-    vector<typename DerivedF::Index> temp = uE2E[ui];
-    for(size_t fei = 0;fei<uE2E[ui].size();fei++)
-    {
-      uE2E[ui][fei] = temp[IM(fei)];
-      const auto & fe = uE2E[ui][fei];
-      diIM(fe) = fei;
-      dicons(fe) = cons[IM(fei)];
-    }
-
   }
 
   vector<vector<vector<Index > > > TT,_1;
@@ -197,6 +117,7 @@ IGL_INLINE void igl::outer_hull(
   vector<MatrixXG> vG(ncc);
   vector<MatrixXJ> vJ(ncc);
   vector<MatrixXJ> vIM(ncc);
+  size_t face_count = 0;
   for(size_t id = 0;id<ncc;id++)
   {
     vIM[id].resize(counts[id],1);
@@ -233,6 +154,9 @@ IGL_INLINE void igl::outer_hull(
     outer_facet(V,F,N,IM,f,f_flip);
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"outer facet: "<<f<<endl;
+  cout << V.row(F(f, 0)) << std::endl;
+  cout << V.row(F(f, 1)) << std::endl;
+  cout << V.row(F(f, 2)) << std::endl;
 #endif
     int FHcount = 1;
     FH[f] = true;
@@ -242,6 +166,8 @@ IGL_INLINE void igl::outer_hull(
     Q.push(f+1*m);
     Q.push(f+2*m);
     flip(f) = f_flip;
+    //std::cout << "face " << face_count++ << ": " << f << std::endl;
+    //std::cout << "f " << F.row(f).array()+1 << std::endl;
     //cout<<"flip("<<f<<") = "<<(flip(f)?"true":"false")<<endl;
 #ifdef IGL_OUTER_HULL_DEBUG
   cout<<"BFS..."<<endl;
@@ -255,6 +181,12 @@ IGL_INLINE void igl::outer_hull(
       const int f = e%m;
       // corner
       const int c = e/m;
+#ifdef IGL_OUTER_HULL_DEBUG
+      std::cout << "edge: " << e << ", ue: " << EMAP(e) << std::endl;
+      std::cout << "face: " << f << std::endl;
+      std::cout << "corner: " << c << std::endl;
+      std::cout << "consistent: " << uE2C[EMAP(e)][diIM[e]] << std::endl;
+#endif
       // Should never see edge again...
       if(EH[e] == true)
       {
@@ -267,6 +199,21 @@ IGL_INLINE void igl::outer_hull(
       const int fd = flip(f)?F(f,(c+1)%3):F(f,(c+2)%3);
       // edge valence
       const size_t val = uE2E[EMAP(e)].size();
+#ifdef IGL_OUTER_HULL_DEBUG
+      std::cout << "vd: " << V.row(fd) << std::endl;
+      std::cout << "vs: " << V.row(fs) << std::endl;
+      std::cout << "edge: " << V.row(fd) - V.row(fs) << std::endl;
+      for (size_t i=0; i<val; i++) {
+          if (i == diIM(e)) {
+              std::cout << "* ";
+          } else {
+              std::cout << "  ";
+          }
+          std::cout << i << ": "
+              << " (e: " << uE2E[EMAP(e)][i] << ", f: "
+              << 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 
@@ -282,7 +229,7 @@ IGL_INLINE void igl::outer_hull(
       //const auto es = F(fe0%m,((fe0/m)+1)%3);
 
       // is edge consistent with edge of face used for sorting
-      const int e_cons = (dicons(e) ? 1: -1);
+      const int e_cons = (uE2C[EMAP(e)][diIM(e)] ? 1: -1);
       int nfei = -1;
       // Loop once around trying to find suitable next face
       for(size_t step = 1; step<val+2;step++)
@@ -290,15 +237,16 @@ IGL_INLINE void igl::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)] != di[EMAP(e)][nfei_new]))
+        //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-16 )
+//        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)<<", |"<<
-          di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei_new]<<"| = "<<
-            abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new])
-            <<endl;
+        //cout<<"Next facet: "<<(f+1)<<" --> "<<(nf+1)<<", |"<<
+        //  di[EMAP(e)][diIM(e)]<<" - "<<di[EMAP(e)][nfei_new]<<"| = "<<
+        //    abs(di[EMAP(e)][diIM(e)] - di[EMAP(e)][nfei_new])
+        //    <<endl;
 #endif
         
 
@@ -307,8 +255,16 @@ IGL_INLINE void igl::outer_hull(
           if(!FH[nf])
           {
             nfei = nfei_new;
+          //} else {
+          //    std::cout << "skipping face " << nfei_new << " because it is seen before"
+          //        << std::endl;
           }
           break;
+        //} else {
+        //    std::cout << di[EMAP(e)][diIM(e)].transpose() << std::endl;
+        //    std::cout << di[EMAP(e)][diIM(nfei_new)].transpose() << std::endl;
+        //    std::cout << "skipping face " << nfei_new << " with identical dihedral angle" 
+        //        << std::endl;
         }
 //#ifdef IGL_OUTER_HULL_DEBUG
 //        cout<<"Skipping co-planar facet: "<<(f+1)<<" --> "<<(nf+1)<<endl;
@@ -391,6 +347,8 @@ IGL_INLINE void igl::outer_hull(
         }
 #endif
         FH[nf] = true;
+        //std::cout << "face " << face_count++ << ": " << nf << std::endl;
+        //std::cout << "f " << F.row(nf).array()+1 << std::endl;
         FHcount++;
         // corner of neighbor
         const int nc = max_ne/m;
@@ -438,18 +396,19 @@ IGL_INLINE void igl::outer_hull(
   // Is A inside B? Assuming A and B are consistently oriented but closed and
   // non-intersecting.
   const auto & is_component_inside_other = [](
-    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> & V,
     const MatrixXV & BC,
     const MatrixXG & A,
     const MatrixXJ & AJ,
     const MatrixXG & B)->bool
   {
+      typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Matrix;
     const auto & bounding_box = [](
-      const Eigen::PlainObjectBase<DerivedV> & V,
+      const Matrix & V,
       const MatrixXG & F)->
-      MatrixXV
+      Matrix
     {
-      MatrixXV BB(2,3);
+      Matrix BB(2,3);
       BB<<
          1e26,1e26,1e26,
         -1e26,-1e26,-1e26;
@@ -467,8 +426,8 @@ IGL_INLINE void igl::outer_hull(
     };
     // A lot of the time we're dealing with unrelated, distant components: cull
     // them.
-    MatrixXV ABB = bounding_box(V,A);
-    MatrixXV BBB = bounding_box(V,B);
+    Matrix ABB = bounding_box(V,A);
+    Matrix BBB = bounding_box(V,B);
     if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
         (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
     {
@@ -479,34 +438,33 @@ IGL_INLINE void igl::outer_hull(
     // POTENTIAL ROBUSTNESS WEAK AREA
     ////////////////////////////////////////////////////////////////////////
     //
-    // q could be so close (<~1e-16) to B that the winding number is not a robust way to
+
+    // winding_number_3 expects colmajor
+    // q could be so close (<~1e-15) to B that the winding number is not a robust way to
     // determine inside/outsideness. We could try to find a _better_ q which is
     // farther away, but couldn't they all be bad?
-    MatrixXV q = BC.row(AJ(0));
+    double q[3] = {
+        CGAL::to_double(BC(AJ(0), 0)),
+        CGAL::to_double(BC(AJ(0), 1)),
+        CGAL::to_double(BC(AJ(0), 2)) };
     // In a perfect world, it's enough to test a single point.
     double w;
-
-    // winding_number_3 expects colmajor
-    const typename DerivedV::Scalar * Vdata;
+    const double * Vdata;
     Vdata = V.data();
-    Matrix<
-      typename DerivedV::Scalar,
-      DerivedV::RowsAtCompileTime,
-      DerivedV::ColsAtCompileTime,
-      ColMajor> Vcol;
-    if(DerivedV::IsRowMajor)
-    {
-      // copy to convert to colmajor
-      Vcol = V;
-      Vdata = Vcol.data();
-    }
     winding_number_3(
       Vdata,V.rows(),
       B.data(),B.rows(),
-      q.data(),1,&w);
-    return fabs(w)>0.5;
+      q,1,&w);
+    return w > 0.5 || w < -0.5;
   };
 
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Vcol(V.rows(), V.cols());
+  for (size_t i=0; i<V.rows(); i++) {
+      for (size_t j=0; j<V.cols(); j++) {
+          Vcol(i, j) = CGAL::to_double(V(i, j));
+      }
+  }
+
   // Reject components which are completely inside other components
   vector<bool> keep(ncc,true);
   size_t nG = 0;
@@ -519,7 +477,7 @@ IGL_INLINE void igl::outer_hull(
       {
         continue;
       }
-      const bool inside = is_component_inside_other(V,BC,vG[id],vJ[id],vG[oid]);
+      const bool inside = is_component_inside_other(Vcol,BC,vG[id],vJ[id],vG[oid]);
 #ifdef IGL_OUTER_HULL_DEBUG
       cout<<id<<" is inside "<<oid<<" ? "<<inside<<endl;
 #endif
@@ -567,6 +525,8 @@ IGL_INLINE void igl::outer_hull(
   return outer_hull(V,F,N,G,J,flip);
 }
 
+
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::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> >&);

include/igl/sort_angles.cpp

@@ -0,0 +1,99 @@
+#include "sort_angles.h"
+#include <algorithm>
+
+template <typename DerivedM, typename DerivedR>
+IGL_INLINE void igl::sort_angles(
+        const Eigen::PlainObjectBase<DerivedM>& M,
+        Eigen::PlainObjectBase<DerivedR>& R) {
+    const size_t num_rows = M.rows();
+    const size_t num_cols = M.cols();
+    assert(num_cols >= 2);
+
+    R.resize(num_rows);
+    R.setLinSpaced(num_rows, 0, num_rows-1);
+
+    //              |
+    // (pi/2, pi)   | (0, pi/2)
+    //              |
+    // -------------+--------------
+    //              |
+    // (-pi, -pi/2) | (-pi/2, 0)
+    //              |
+    auto comp = [&](size_t i, size_t j) {
+        auto yi = M(i, 0);
+        auto xi = M(i, 1);
+        auto yj = M(j, 0);
+        auto xj = M(j, 1);
+
+        if (xi == xj && yi == yj) {
+            for (size_t idx=2; idx<num_cols; idx++) {
+                auto i_val = M(i, idx);
+                auto j_val = M(j, idx);
+                if (i_val != j_val) {
+                    return i_val < j_val;
+                }
+            }
+            // If the entire rows are equal, use the row index.
+            return i < j;
+        }
+
+        if (xi >= 0 && yi >= 0) {
+            if (xj >=0 && yj >= 0) {
+                if (xi != xj) {
+                    return xi > xj;
+                } else {
+                    return yi < yj;
+                }
+            } else if (xj < 0 && yj >= 0) {
+                return true;
+            } else if (xj < 0 && yj < 0) {
+                return false;
+            } else {
+                return false;
+            }
+        } else if (xi < 0 && yi >= 0) {
+            if (xj >= 0 && yj >= 0) {
+                return false;
+            } else if (xj < 0 && yj >= 0) {
+                if (xi != xj) {
+                    return xi > xj;
+                } else {
+                    return yi > yj;
+                }
+            } else if (xj < 0 && yj < 0) {
+                return false;
+            } else {
+                return false;
+            }
+        } else if (xi < 0 && yi < 0) {
+            if (xj >= 0 && yj >= 0) {
+                return true;
+            } else if (xj < 0 && yj >= 0) {
+                return true;
+            } else if (xj < 0 && yj < 0) {
+                if (xi != xj) {
+                    return xi < xj;
+                } else {
+                    return yi > yj;
+                }
+            } else {
+                return true;
+            }
+        } else {
+            if (xj >= 0 && yj >= 0) {
+                return true;
+            } else if (xj < 0 && yj >= 0) {
+                return true;
+            } else if (xj < 0 && yj < 0) {
+                return false;
+            } else {
+                if (xi != xj) {
+                    return xi < xj;
+                } else {
+                    return yi < yj;
+                }
+            }
+        }
+    };
+    std::sort(R.data(), R.data() + num_rows, comp);
+}

include/igl/sort_angles.h

@@ -0,0 +1,27 @@
+#ifndef SORT_ANGLES_H
+#define SORT_ANGLES_H
+
+#include "igl_inline.h"
+
+namespace igl {
+    // Sort angles in ascending order in a numerically robust way.
+    //
+    // Instead of computing angles using atan2(y, x), sort directly on (y, x).
+    //
+    // Inputs:
+    //   M: m by n matrix of scalars. (n >= 2).  Assuming the first column of M
+    //      contains values for y, and the second column is x.  Using the rest
+    //      of the columns as tie-breaker.
+    //   R: an array of m indices.  M.row(R[i]) contains the i-th smallest
+    //      angle.
+    template<typename DerivedM, typename DerivedR>
+    IGL_INLINE void sort_angles(
+            const Eigen::PlainObjectBase<DerivedM>& M,
+            Eigen::PlainObjectBase<DerivedR>& R);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "sort_angles.cpp"
+#endif
+
+#endif

include/igl/winding_number.cpp

@@ -43,15 +43,15 @@ IGL_INLINE void igl::winding_number(
   }
 }
 
-template <typename DerivedF>
+template <typename Scalar, typename DerivedF>
 IGL_INLINE void igl::winding_number_3(
-  const double * V,
+  const Scalar * V,
   const int n,
   const DerivedF * F,
   const int m,
-  const double * O,
+  const Scalar * O,
   const int no,
-  double * S)
+  Scalar * S)
 {
   // Initialize output
   // loop over origins
@@ -68,7 +68,7 @@ IGL_INLINE void igl::winding_number_3(
   for(int f = 0;f<m;f++)
   {
     // Gather corners 
-    double C[3][3];
+    Scalar C[3][3];
     // loop around triangle
     for(int t=0;t<3;t++)
     {
@@ -84,8 +84,8 @@ IGL_INLINE void igl::winding_number_3(
     for(int o = 0;o<no;o++)
     {
       // Gather vectors to corners
-      double v[3][3];
-      double vl[3];
+      Scalar v[3][3];
+      Scalar vl[3];
       // loop around triangle
       for(int t=0;t<3;t++)
       {
@@ -106,7 +106,7 @@ IGL_INLINE void igl::winding_number_3(
       }
       //printf("\n");
       // Compute determinant
-      double detf = 
+      Scalar detf = 
         v[0][0]*v[1][1]*v[2][2]+
         v[1][0]*v[2][1]*v[0][2]+
         v[2][0]*v[0][1]*v[1][2]-
@@ -114,7 +114,7 @@ IGL_INLINE void igl::winding_number_3(
         v[1][0]*v[0][1]*v[2][2]-
         v[0][0]*v[2][1]*v[1][2];
       // Compute pairwise dotproducts
-      double dp[3];
+      Scalar dp[3];
       dp[0] = v[1][0]*v[2][0];
       dp[0] += v[1][1]*v[2][1];
       dp[0] += v[1][2]*v[2][2];
@@ -126,7 +126,7 @@ IGL_INLINE void igl::winding_number_3(
       dp[2] += v[0][2]*v[1][2];
       // Compute winding number
       // Only divide by TWO_PI instead of 4*pi because there was a 2 out front
-      double val = atan2(detf,
+      Scalar val = atan2(detf,
         vl[0]*vl[1]*vl[2] + 
         dp[0]*vl[0] +
         dp[1]*vl[1] +

include/igl/winding_number.h

@@ -39,15 +39,15 @@ namespace igl
   //   no  number of origins
   // Outputs:
   //   S  no by 1 list of winding numbers
-  template <typename DerivedF>
+  template <typename Scalar, typename DerivedF>
   IGL_INLINE void winding_number_3(
-    const double * V,
+    const Scalar * V,
     const int n,
     const DerivedF * F,
     const int m,
-    const double * O,
+    const Scalar * O,
     const int no,
-    double * S);
+    Scalar * S);
   //// Only one evaluation origin
   //template <typename DerivedF>
   //IGL_INLINE void winding_number_3(