From ee6ea7e3d026d96f0e7ec965459cbff1d7a68056 Mon Sep 17 00:00:00 2001
From: Alec Jacobson <alecjacobson@gmail.com>
Date: Tue, 11 Nov 2014 18:19:10 -0500
Subject: [PATCH] fixed bug in outer_hull causing incorrect booleans for
 non-overlapping components

---
 examples/skeleton-posing/example.cpp |  52 +++-
 include/igl/WindingNumberTree.h      |   4 +-
 include/igl/barycenter.cpp           |   1 +
 include/igl/exterior_edges.h         |   2 +-
 include/igl/facet_components.cpp     |  19 +-
 include/igl/facet_components.h       |   6 +-
 include/igl/outer_hull.cpp           | 340 ++++++++++++++++-----------
 include/igl/sort.cpp                 |   1 +
 include/igl/unique_edge_map.h        |   4 +-
 include/igl/winding_number.h         |  18 +-
 10 files changed, 282 insertions(+), 165 deletions(-)

diff --git a/examples/skeleton-posing/example.cpp b/examples/skeleton-posing/example.cpp
index 7f6b8d8dd..9f1b57ab0 100644
--- a/examples/skeleton-posing/example.cpp
+++ b/examples/skeleton-posing/example.cpp
@@ -106,7 +106,7 @@ bool is_rotating = false;
 bool centroid_is_visible = true;
 int down_x,down_y;
 igl::Camera down_camera;
-std::string output_prefix;
+std::string output_weights_filename,output_pose_prefix;
 
 struct CameraAnimation
 {
@@ -594,13 +594,13 @@ void redo()
   }
 }
 
-bool save()
+bool save_pose()
 {
   using namespace std;
   using namespace igl;
   using namespace Eigen;
   string output_filename;
-  next_filename(output_prefix,4,".dmat",output_filename);
+  next_filename(output_pose_prefix,4,".dmat",output_filename);
   MatrixXd T;
   forward_kinematics(C,BE,P,s.mouse.rotations(),T);
   if(writeDMAT(output_filename,T))
@@ -614,6 +614,23 @@ bool save()
   }
 }
 
+bool save_weights()
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  if(writeDMAT(output_weights_filename,W))
+  {
+    cout<<GREENGIN("Current weights written to "+
+      output_weights_filename+".")<<endl;
+    return true;
+  }else
+  {
+    cout<<REDRUM("Writing to "+output_weights_filename+" failed.")<<endl;
+    return false;
+  }
+}
+
 void key(unsigned char key, int mouse_x, int mouse_y)
 {
   using namespace std;
@@ -660,7 +677,13 @@ void key(unsigned char key, int mouse_x, int mouse_y)
     case 'S':
     case 's':
     {
-      save();
+      save_pose();
+      break;
+    }
+    case 'W':
+    case 'w':
+    {
+      save_weights();
       break;
     }
     case 'z':
@@ -848,11 +871,11 @@ int main(int argc, char * argv[])
   string filename = "../shared/cheburashka.off";
   string skel_filename = "../shared/cheburashka.tgf";
   string weights_filename = "";
-  output_prefix = "";
+  output_pose_prefix = "";
   switch(argc)
   {
     case 5:
-      output_prefix = argv[4];
+      output_pose_prefix = argv[4];
       //fall through
     case 4:
       weights_filename = argv[3];
@@ -872,9 +895,10 @@ int main(int argc, char * argv[])
   cout<<"⌥ +[Click] and [drag]  Rotate secene."<<endl;
   cout<<"⌫                      Delete selected node(s) and incident bones."<<endl;
   cout<<"D,d                    Deselect all."<<endl;
-  cout<<"S,s                    Save current pose."<<endl;
   cout<<"R                      Reset selected rotation."<<endl;
   cout<<"r                      Reset all rotations."<<endl;
+  cout<<"S,s                    Save current pose."<<endl;
+  cout<<"W,w                    Save current weights."<<endl;
   cout<<"Z,z                    Snap to canonical view."<<endl;
   cout<<"⌘ Z                    Undo."<<endl;
   cout<<"⇧ ⌘ Z                  Redo."<<endl;
@@ -883,15 +907,21 @@ int main(int argc, char * argv[])
   string dir,_1,_2,name;
   read_triangle_mesh(filename,V,F,dir,_1,_2,name);
 
-  if(output_prefix.size() == 0)
+  if(output_weights_filename.size() == 0)
   {
-    output_prefix = dir+"/"+name+"-pose-";
+    output_weights_filename = dir+"/"+name+"-weights.dmat";
+  }
+  if(output_pose_prefix.size() == 0)
+  {
+    output_pose_prefix = dir+"/"+name+"-pose-";
   }
 
   {
     string output_filename;
-    next_filename(output_prefix,4,".dmat",output_filename);
-    cout<<BLUEGIN("Output set to start with "<<output_filename)<<endl;
+    next_filename(output_pose_prefix,4,".dmat",output_filename);
+    cout<<BLUEGIN("Output weights set to start with "<<
+      output_weights_filename)<<endl;
+    cout<<BLUEGIN("Output poses set to start with "<<output_filename)<<endl;
   }
 
   // Read in skeleton and precompute hierarchy
diff --git a/include/igl/WindingNumberTree.h b/include/igl/WindingNumberTree.h
index 6ea478dd4..1bc6cec29 100644
--- a/include/igl/WindingNumberTree.h
+++ b/include/igl/WindingNumberTree.h
@@ -170,7 +170,7 @@ inline void igl::WindingNumberTree<Point>::set_mesh(
   // Q: Would we gain even more by remove almost exactly duplicate vertices?
   Eigen::MatrixXi SF,SVI,SVJ;
   igl::remove_duplicate_vertices(_V,_F,0.0,SV,SVI,SVJ,F);
-  triangle_fan(exterior_edges(F),cap);
+  triangle_fan(igl::exterior_edges(F),cap);
   V = SV;
 }
 
@@ -183,7 +183,7 @@ inline igl::WindingNumberTree<Point>::WindingNumberTree(
   V(parent.V),
   SV(),
   F(_F),
-  cap(triangle_fan(exterior_edges(_F)))
+  cap(triangle_fan(igl::exterior_edges(_F)))
 {
 }
 
diff --git a/include/igl/barycenter.cpp b/include/igl/barycenter.cpp
index 79738c822..2b5805078 100644
--- a/include/igl/barycenter.cpp
+++ b/include/igl/barycenter.cpp
@@ -37,4 +37,5 @@ template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::M
 template void igl::barycenter<Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(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, 3, 0, -1, 3> >&);
 template void igl::barycenter<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::barycenter<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 #endif
diff --git a/include/igl/exterior_edges.h b/include/igl/exterior_edges.h
index b003f86ee..2acea1955 100644
--- a/include/igl/exterior_edges.h
+++ b/include/igl/exterior_edges.h
@@ -17,7 +17,7 @@ namespace igl
     const Eigen::MatrixXi & F,
     Eigen::MatrixXi & E);
   // Inline version
-  Eigen::MatrixXi exterior_edges( const Eigen::MatrixXi & F);
+  IGL_INLINE Eigen::MatrixXi exterior_edges( const Eigen::MatrixXi & F);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "exterior_edges.cpp"
diff --git a/include/igl/facet_components.cpp b/include/igl/facet_components.cpp
index ff2a78ce8..e3a00d6f4 100644
--- a/include/igl/facet_components.cpp
+++ b/include/igl/facet_components.cpp
@@ -18,10 +18,12 @@ IGL_INLINE void igl::facet_components(
 
 template <
   typename TTIndex, 
-  typename DerivedC>
+  typename DerivedC,
+  typename Derivedcounts>
 IGL_INLINE void igl::facet_components(
   const std::vector<std::vector<std::vector<TTIndex > > > & TT,
-  Eigen::PlainObjectBase<DerivedC> & C)
+  Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<Derivedcounts> & counts)
 {
   using namespace std;
   using namespace igl;
@@ -30,12 +32,14 @@ IGL_INLINE void igl::facet_components(
   C.resize(m,1);
   vector<bool> seen(m,false);
   Index id = 0;
+  vector<Index> vcounts;
   for(Index g = 0;g<m;g++)
   {
     if(seen[g])
     {
       continue;
     }
+    vcounts.push_back(0);
     queue<Index> Q;
     Q.push(g);
     while(!Q.empty())
@@ -47,6 +51,7 @@ IGL_INLINE void igl::facet_components(
         continue;
       }
       seen[f] = true;
+      vcounts[id]++;
       C(f,0) = id;
       // Face f's neighbor lists opposite opposite each corner
       for(const auto & c : TT[f])
@@ -63,9 +68,17 @@ IGL_INLINE void igl::facet_components(
     }
     id++;
   }
+  assert(id == vcounts.size());
+  const size_t ncc = vcounts.size();
+  assert(C.maxCoeff()+1 == ncc);
+  counts.resize(ncc,1);
+  for(size_t i = 0;i<ncc;i++)
+  {
+    counts(i) = vcounts[i];
+  }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::facet_components<long, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
+template void igl::facet_components<long, Eigen::Matrix<long, -1, 1, 0, -1, 1>, Eigen::Matrix<long, -1, 1, 0, -1, 1> >(std::vector<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > >, std::allocator<std::vector<std::vector<long, std::allocator<long> >, std::allocator<std::vector<long, std::allocator<long> > > > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<long, -1, 1, 0, -1, 1> >&);
 #endif
diff --git a/include/igl/facet_components.h b/include/igl/facet_components.h
index 254486f97..96ea5db03 100644
--- a/include/igl/facet_components.h
+++ b/include/igl/facet_components.h
@@ -22,10 +22,12 @@ namespace igl
   //   C  #F list of connected component ids
   template <
     typename TTIndex, 
-    typename DerivedC>
+    typename DerivedC,
+    typename Derivedcounts>
   IGL_INLINE void facet_components(
     const std::vector<std::vector<std::vector<TTIndex > > > & TT,
-    Eigen::PlainObjectBase<DerivedC> & C);
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<Derivedcounts> & counts);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "facet_components.cpp"
diff --git a/include/igl/outer_hull.cpp b/include/igl/outer_hull.cpp
index 63856799c..6b7fea2e1 100644
--- a/include/igl/outer_hull.cpp
+++ b/include/igl/outer_hull.cpp
@@ -1,9 +1,10 @@
 #include "outer_hull.h"
 #include "outer_facet.h"
 #include "facet_components.h"
-
+#include "winding_number.h"
 #include "triangle_triangle_adjacency.h"
 #include "unique_edge_map.h"
+#include "barycenter.h"
 #include "per_face_normals.h"
 #include "all_edges.h"
 #include "colon.h"
@@ -33,10 +34,10 @@ IGL_INLINE void igl::outer_hull(
   using namespace igl;
   typedef typename DerivedF::Index Index;
   Matrix<Index,DerivedF::RowsAtCompileTime,1> C;
+  typedef Matrix<typename DerivedV::Scalar,Dynamic,DerivedV::ColsAtCompileTime> MatrixXV;
   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 Derivedflip::Scalar,Dynamic,Derivedflip::ColsAtCompileTime> MatrixXflip;
   const Index m = F.rows();
 
   typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
@@ -48,8 +49,10 @@ IGL_INLINE void igl::outer_hull(
 
   vector<vector<vector<Index > > > TT,_1;
   triangle_triangle_adjacency(E,EMAP,uE2E,false,TT,_1);
-  facet_components(TT,C);
-  const Index ncc = C.maxCoeff()+1;
+  VectorXI counts;
+  facet_components(TT,C,counts);
+  assert(C.maxCoeff()+1 == counts.rows());
+  const size_t ncc = counts.rows();
   G.resize(0,F.cols());
   J.resize(0,1);
   flip.setConstant(m,1,false);
@@ -64,162 +67,229 @@ IGL_INLINE void igl::outer_hull(
   vector<bool> EH(3*m,false);
   vector<MatrixXG> vG(ncc);
   vector<MatrixXJ> vJ(ncc);
-
-  // Total size of G (and J)
-  size_t nG = 0;
-  // This is O( (n+m) * nnc) !
-  for(Index id = 0;id<ncc;id++)
+  vector<MatrixXJ> vIM(ncc);
+  for(size_t id = 0;id<ncc;id++)
   {
-    // Determine number of facets in component id
-    Index num_id = 0;
-    for(Index f = 0;f<m;f++)
+    vIM[id].resize(counts[id],1);
+  }
+  // current index into each IM
+  vector<size_t> g(ncc,0);
+  // place order of each face in its respective component
+  for(Index f = 0;f<m;f++)
+  {
+    vIM[C(f)](g[C(f)]++) = f;
+  }
+
+  // assumes that "resolve" has handled any coplanar cases correctly and nearly
+  // coplanar cases can be sorted based on barycenter.
+  MatrixXV BC;
+  barycenter(V,F,BC);
+
+  for(Index id = 0;id<(Index)ncc;id++)
+  {
+    auto & IM = vIM[id];
+    // starting face that's guaranteed to be on the outer hull and in this
+    // component
+    int f;
+    bool f_flip;
+    outer_facet(V,F,N,IM,f,f_flip);
+    int FHcount = 0;
+    // Q contains list of face edges to continue traversing upong
+    queue<int> Q;
+    Q.push(f+0*m);
+    Q.push(f+1*m);
+    Q.push(f+2*m);
+    flip[f] = f_flip;
+    while(!Q.empty())
     {
-      if(C(f) == id)
+      // face-edge
+      const int e = Q.front();
+      Q.pop();
+      // face
+      const int f = e%m;
+      // corner
+      const int c = e/m;
+      // Should never see edge again...
+      if(EH[e] == true)
       {
-        num_id++;
+        continue;
       }
-    }
-    //MatrixXF Fc(num_id,F.cols());
-    MatrixXJ IM(num_id,1);
-    if(C.maxCoeff() == 0)
-    {
-      assert(num_id == F.rows());
-      //Fc = F;
-      IM = MatrixXJ::LinSpaced(F.rows(),0,F.rows()-1);
-    }else
-    {
-      int g = 0;
-      for(Index f = 0;f<m;f++)
+      EH[e] = true;
+      // first time seeing face
+      if(!FH[f])
       {
-        if(C(f) == id)
-        {
-          //Fc.row(g) = F.row(f);
-          IM(g) = f;
-          g++;
-        }
+        FH[f] = true;
+        FHcount++;
       }
-    }
-    {
-      // starting face that's guaranteed to be on the outer hull and in this
-      // component
-      int f;
-      bool f_flip;
-      outer_facet(V,F,N,IM,f,f_flip);
-      // Q contains list of face edges to continue traversing upong
-      queue<int> Q;
-      Q.push(f+0*m);
-      Q.push(f+1*m);
-      Q.push(f+2*m);
-      flip[f] = f_flip;
-      int FHcount = 0;
-      while(!Q.empty())
+      // find overlapping face-edges
+      const auto & neighbors = uE2E[EMAP(e)];
+      const auto & fN = (flip[f]?-1.:1.)*N.row(f);
+      // source of edge according to f
+      const int fs = flip[f]?F(f,(c+2)%3):F(f,(c+1)%3);
+      // destination of edge according to f
+      const int fd = flip[f]?F(f,(c+1)%3):F(f,(c+2)%3);
+      const auto & eV = (V.row(fd)-V.row(fs)).normalized();
+      // Loop over and find max dihedral angle
+      typename DerivedV::Scalar max_di = -1;
+      int max_ne = -1;
+      typename Eigen::Matrix< typename DerivedV::Scalar, 1, 3> max_nN;
+      for(const auto & ne : neighbors)
       {
-        // face-edge
-        const int e = Q.front();
-        Q.pop();
-        // face
-        const int f = e%m;
-        // corner
-        const int c = e/m;
-        // Should never see edge again...
-        if(EH[e] == true)
+        const int nf = ne%m;
+        if(nf == f)
         {
           continue;
         }
-        EH[e] = true;
-        // first time seeing face
-        if(!FH[f])
+        const int nc = ne/m;
+        // are faces consistently oriented
+        //const int ns = F(nf,(nc+1)%3);
+        const int nd = F(nf,(nc+2)%3);
+        const bool cons = (flip[f]?fd:fs) == nd;
+        const auto & nN = (cons? (flip[f]?-1:1.) : (flip[f]?1.:-1.) )*N.row(nf);
+        const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN));
+        if(ndi>=max_di)
         {
-          FH[f] = true;
-          FHcount++;
-        }
-        // find overlapping face-edges
-        const auto & neighbors = uE2E[EMAP(e)];
-        const auto & fN = (flip[f]?-1.:1.)*N.row(f);
-        // source of edge according to f
-        const int fs = flip[f]?F(f,(c+2)%3):F(f,(c+1)%3);
-        // destination of edge according to f
-        const int fd = flip[f]?F(f,(c+1)%3):F(f,(c+2)%3);
-        const auto & eV = (V.row(fd)-V.row(fs)).normalized();
-        // Loop over and find max dihedral angle
-        typename DerivedV::Scalar max_di = -1;
-        int max_ne = -1;
-        typename Eigen::Matrix< typename DerivedV::Scalar, 1, 3> max_nN;
-        for(const auto & ne : neighbors)
-        {
-          const int nf = ne%m;
-          if(nf == f)
-          {
-            continue;
-          }
-          const int nc = ne/m;
-          // are faces consistently oriented
-          //const int ns = F(nf,(nc+1)%3);
-          const int nd = F(nf,(nc+2)%3);
-          const bool cons = (flip[f]?fd:fs) == nd;
-          const auto & nN = (cons? (flip[f]?-1:1.) : (flip[f]?1.:-1.) )*N.row(nf);
-          const auto & ndi = M_PI - atan2( fN.cross(nN).dot(eV), fN.dot(nN));
-          if(ndi>=max_di)
-          {
-            max_ne = ne;
-            max_di = ndi;
-            max_nN = nN;
-          }
-        }
-        if(max_ne>=0)
-        {
-          const int nf = max_ne%m;
-          const int nc = max_ne/m;
-          const int nd = F(nf,(nc+2)%3);
-          const bool cons = (flip[f]?fd:fs) == nd;
-          flip[nf] = (cons ? flip[f] : !flip[f]);
-          const int ne1 = nf+((nc+1)%3)*m;
-          const int ne2 = nf+((nc+2)%3)*m;
-          if(!EH[ne1])
-          {
-            Q.push(ne1);
-          }
-          if(!EH[ne2])
-          {
-            Q.push(ne2);
-          }
+          max_ne = ne;
+          max_di = ndi;
+          max_nN = nN;
         }
       }
-      
+      if(max_ne>=0)
       {
-        vG[id].resize(FHcount,3);
-        vJ[id].resize(FHcount,1);
-        nG += FHcount;
-        size_t h = 0;
-        for(int i = 0;i<num_id;i++)
+        const int nf = max_ne%m;
+        const int nc = max_ne/m;
+        const int nd = F(nf,(nc+2)%3);
+        const bool cons = (flip[f]?fd:fs) == nd;
+        flip[nf] = (cons ? flip[f] : !flip[f]);
+        const int ne1 = nf+((nc+1)%3)*m;
+        const int ne2 = nf+((nc+2)%3)*m;
+        if(!EH[ne1])
         {
-          const size_t f = IM(i);
-          //if(f_flip)
-          //{
-          //  flip[f] = !flip[f];
-          //}
-          if(FH[f])
-          {
-            vG[id].row(h) = (flip[f]?F.row(f).reverse().eval():F.row(f));
-            vJ[id](h,0) = f;
-            h++;
-          }
+          Q.push(ne1);
+        }
+        if(!EH[ne2])
+        {
+          Q.push(ne2);
         }
-        assert(h == FHcount);
       }
     }
+    
+    {
+      vG[id].resize(FHcount,3);
+      vJ[id].resize(FHcount,1);
+      //nG += FHcount;
+      size_t h = 0;
+      assert(counts(id) == IM.rows());
+      for(int i = 0;i<counts(id);i++)
+      {
+        const size_t f = IM(i);
+        //if(f_flip)
+        //{
+        //  flip[f] = !flip[f];
+        //}
+        if(FH[f])
+        {
+          vG[id].row(h) = (flip[f]?F.row(f).reverse().eval():F.row(f));
+          vJ[id](h,0) = f;
+          h++;
+        }
+      }
+      assert(h == FHcount);
+    }
   }
+
+  // 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 MatrixXV & BC,
+    const MatrixXG & A,
+    const MatrixXJ & AJ,
+    const MatrixXG & B)->bool
+  {
+    const auto & bounding_box = [](
+      const Eigen::PlainObjectBase<DerivedV> & V,
+      const MatrixXG & F)->
+      MatrixXV
+    {
+      MatrixXV BB(2,3);
+      BB<<
+         1e26,1e26,1e26,
+        -1e26,-1e26,-1e26;
+      const size_t m = F.rows();
+      for(size_t f = 0;f<m;f++)
+      {
+        for(size_t c = 0;c<3;c++)
+        {
+          const auto & vfc = V.row(F(f,c));
+          BB.row(0) = BB.row(0).array().min(vfc.array()).eval();
+          BB.row(1) = BB.row(1).array().max(vfc.array()).eval();
+        }
+      }
+      return BB;
+    };
+    // 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);
+    if( (BBB.row(0)-ABB.row(1)).maxCoeff()>0  ||
+        (ABB.row(0)-BBB.row(1)).maxCoeff()>0 )
+    {
+      // bounding boxes do not overlap
+      return false;
+    }
+    ////////////////////////////////////////////////////////////////////////
+    // POTENTIAL ROBUSTNESS WEAK AREA
+    ////////////////////////////////////////////////////////////////////////
+    //
+    // q could be so close (<~1e-16) 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(1));
+    // In a perfect world, it's enough to test a single point.
+    double w;
+    winding_number_3(
+      V.data(),V.rows(),
+      B.data(),B.rows(),
+      q.data(),1,&w);
+    return fabs(w)>0.5;
+  };
+
+  // Reject components which are completely inside other components
+  vector<bool> keep(ncc,true);
+  size_t nG = 0;
+  // This is O( ncc * ncc * m)
+  for(size_t id = 0;id<ncc;id++)
+  {
+    for(size_t oid = 0;oid<ncc;oid++)
+    {
+      if(id == oid)
+      {
+        continue;
+      }
+      keep[id] = keep[id] && 
+        !is_component_inside_other(V,BC,vG[id],vJ[id],vG[oid]);
+    }
+    if(keep[id])
+    {
+      nG += vJ[id].rows();
+    }
+  }
+
   // collect G and J across components
   G.resize(nG,3);
   J.resize(nG,1);
   {
     size_t off = 0;
-    for(Index id = 0;id<ncc;id++)
+    for(Index id = 0;id<(Index)ncc;id++)
     {
-      assert(vG[id].rows() == vJ[id].rows());
-      G.block(off,0,vG[id].rows(),vG[id].cols()) = vG[id];
-      J.block(off,0,vJ[id].rows(),vJ[id].cols()) = vJ[id];
-      off += vG[id].rows();
+      if(keep[id])
+      {
+        assert(vG[id].rows() == vJ[id].rows());
+        G.block(off,0,vG[id].rows(),vG[id].cols()) = vG[id];
+        J.block(off,0,vJ[id].rows(),vJ[id].cols()) = vJ[id];
+        off += vG[id].rows();
+      }
     }
   }
 }
diff --git a/include/igl/sort.cpp b/include/igl/sort.cpp
index 0dce20af6..814381f91 100644
--- a/include/igl/sort.cpp
+++ b/include/igl/sort.cpp
@@ -221,4 +221,5 @@ template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<do
 template void igl::sort<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<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::sort<Eigen::Matrix<double, -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&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif
diff --git a/include/igl/unique_edge_map.h b/include/igl/unique_edge_map.h
index 5cbb17e20..e3583e310 100644
--- a/include/igl/unique_edge_map.h
+++ b/include/igl/unique_edge_map.h
@@ -1,5 +1,5 @@
-#ifndef IGL_EXTERIOR_EDGES_H
-#define IGL_EXTERIOR_EDGES_H
+#ifndef IGL_UNIQUE_EDGE_MAP_H
+#define IGL_UNIQUE_EDGE_MAP_H
 #include "igl_inline.h"
 #include <Eigen/Dense>
 #include <vector>
diff --git a/include/igl/winding_number.h b/include/igl/winding_number.h
index f3eb9e32f..7629f3b11 100644
--- a/include/igl/winding_number.h
+++ b/include/igl/winding_number.h
@@ -39,15 +39,15 @@ namespace igl
     const double * O,
     const int no,
     double * S);
-  // Only one evaluation origin
-  template <typename DerivedF>
-  IGL_INLINE void winding_number_3(
-    const double * V,
-    const int n,
-    const DerivedF * F,
-    const int m,
-    const double * O,
-    double * S);
+  //// Only one evaluation origin
+  //template <typename DerivedF>
+  //IGL_INLINE void winding_number_3(
+  //  const double * V,
+  //  const int n,
+  //  const DerivedF * F,
+  //  const int m,
+  //  const double * O,
+  //  double * S);
   // 2d
   template <typename DerivedF>
   IGL_INLINE void winding_number_2(