faster tri tri adj, unique simplices

Former-commit-id: c45b42bd184e5f191ddfefc8a5d9385648c57037

include/igl/cgal/SelfIntersectMesh.h

@@ -85,9 +85,11 @@ namespace igl
       // Number of self-intersecting triangle pairs
       typedef typename DerivedF::Index Index;
       Index count;
-      std::vector<std::list<CGAL::Object> > F_objects;
+      typedef std::vector<CGAL::Object> ObjectList;
+      std::vector<ObjectList > F_objects;
       Triangles T;
-      std::list<Index> lIF;
+      typedef std::vector<Index> IndexList;
+      IndexList lIF;
       std::vector<bool> offensive;
       std::vector<Index> offending_index;
       std::vector<Index> offending;
@@ -95,7 +97,7 @@ namespace igl
       typedef std::pair<Index,Index> EMK;
       // Make a short name for the type stored at each edge, the edge map's
       // value
-      typedef std::list<Index> EMV;
+      typedef std::vector<Index> EMV;
       // Make a short name for the edge map
       typedef std::map<EMK,EMV> EdgeMap;
       EdgeMap edge2faces;
@@ -205,11 +207,11 @@ namespace igl
       //   cdt  Contrained delaunay triangulation in projected 2D plane
       inline void projected_delaunay(
           const Triangle_3 & A,
-          const std::list<CGAL::Object> & A_objects_3,
+          const ObjectList & A_objects_3,
           CDT_plus_2 & cdt);
       // Getters:
     public:
-      //const std::list<int>& get_lIF() const{ return lIF;}
+      //const IndexList& get_lIF() const{ return lIF;}
       static inline void box_intersect(
         SelfIntersectMesh * SIM, 
         const SelfIntersectMesh::Box &a, 
@@ -222,6 +224,7 @@ namespace igl
 #include "mesh_to_cgal_triangle_list.h"
 
 #include <igl/REDRUM.h>
+#include <igl/get_seconds.h>
 #include <igl/C_STR.h>
 
 #include <boost/function.hpp>
@@ -322,8 +325,19 @@ inline igl::SelfIntersectMesh<
 {
   using namespace std;
   using namespace Eigen;
+
+  //const auto & tictoc = []()
+  //{
+  //  static double t_start = igl::get_seconds();
+  //  double diff = igl::get_seconds()-t_start;
+  //  t_start += diff;
+  //  return diff;
+  //};
+  //tictoc();
+
   // Compute and process self intersections
   mesh_to_cgal_triangle_list(V,F,T);
+  //cout<<"mesh_to_cgal_triangle_list: "<<tictoc()<<endl;
   // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
   // Create the corresponding vector of bounding boxes
   std::vector<Box> boxes;
@@ -338,6 +352,7 @@ inline igl::SelfIntersectMesh<
   // Leapfrog callback
   boost::function<void(const Box &a,const Box &b)> cb
     = boost::bind(&box_intersect, this, _1,_2);
+  //cout<<"boxes and bind: "<<tictoc()<<endl;
   // Run the self intersection algorithm with all defaults
   try{
     CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
@@ -350,6 +365,7 @@ inline igl::SelfIntersectMesh<
     }
     // Otherwise just fall through
   }
+  //cout<<"box_self_intersection_d: "<<tictoc()<<endl;
 
   // Convert lIF to Eigen matrix
   assert(lIF.size()%2 == 0);
@@ -357,7 +373,7 @@ inline igl::SelfIntersectMesh<
   {
     Index i=0;
     for(
-      typename list<Index>::const_iterator ifit = lIF.begin();
+      typename IndexList::const_iterator ifit = lIF.begin();
       ifit!=lIF.end();
       )
     {
@@ -368,6 +384,7 @@ inline igl::SelfIntersectMesh<
       i++;
     }
   }
+  //cout<<"IF: "<<tictoc()<<endl;
 
   if(params.detect_only)
   {
@@ -380,7 +397,8 @@ inline igl::SelfIntersectMesh<
     typename Eigen::Matrix<typename DerivedFF::Scalar,Dynamic,Dynamic>
     > NF(offending.size());
   // list of new vertices
-  list<Point_3> NV;
+  typedef vector<Point_3> Point_3List;
+  Point_3List NV;
   Index NV_count = 0;
   vector<CDT_plus_2> cdt(offending.size());
   vector<Plane_3> P(offending.size());
@@ -391,7 +409,9 @@ inline igl::SelfIntersectMesh<
   {
     // index in F
     const Index f = offending[o];
-    projected_delaunay(T[f],F_objects[f],cdt[o]);
+    {
+      projected_delaunay(T[f],F_objects[f],cdt[o]);
+    }
     // Q: Is this also delaunay in 3D?
     // A: No, because the projection is affine and delaunay is not affine
     // invariant.
@@ -439,7 +459,7 @@ inline igl::SelfIntersectMesh<
             assert(edge2faces.count(EMK(i,j))==1);
             // loop over neighbors
             for(
-              typename list<Index>::const_iterator nit = edge2faces[EMK(i,j)].begin();
+              typename IndexList::const_iterator nit = edge2faces[EMK(i,j)].begin();
               nit != edge2faces[EMK(i,j)].end() && !found;
               nit++)
             {
@@ -494,6 +514,7 @@ inline igl::SelfIntersectMesh<
       }
     }
   }
+  //cout<<"CDT: "<<tictoc()<<"  "<<t_proj_del<<endl;
   assert(NV_count == (Index)NV.size());
   // Build output
 #ifndef NDEBUG
@@ -535,7 +556,7 @@ inline igl::SelfIntersectMesh<
   {
     Index i = 0;
     for(
-      typename list<Point_3>::const_iterator nvit = NV.begin();
+      typename Point_3List::const_iterator nvit = NV.begin();
       nvit != NV.end();
       nvit++)
     {
@@ -571,7 +592,7 @@ inline igl::SelfIntersectMesh<
   {
     Index v = V.rows();
     for(
-      typename list<Point_3>::const_iterator nvit = NV.begin();
+      typename Point_3List::const_iterator nvit = NV.begin();
       nvit != NV.end();
       nvit++)
     {
@@ -585,6 +606,7 @@ inline igl::SelfIntersectMesh<
       v++;
     }
   }
+  //cout<<"Output + dupes: "<<tictoc()<<endl;
 
   // Q: Does this give the same result as TETGEN?
   // A: For the cow and beast, yes.
@@ -639,7 +661,7 @@ inline void igl::SelfIntersectMesh<
       // Create new list if there is no entry
       if(edge2faces.count(EMK(i,j))==0)
       {
-        edge2faces[EMK(i,j)] = list<Index>();
+        edge2faces[EMK(i,j)] = EMV();
       }
       // append to list
       edge2faces[EMK(i,j)].push_back(f);
@@ -1096,7 +1118,7 @@ inline void igl::SelfIntersectMesh<
   DerivedJ,
   DerivedIM>::projected_delaunay(
   const Triangle_3 & A,
-  const std::list<CGAL::Object> & A_objects_3,
+  const ObjectList & A_objects_3,
   CDT_plus_2 & cdt)
 {
   using namespace std;
@@ -1115,20 +1137,16 @@ inline void igl::SelfIntersectMesh<
     cdt.insert_constraint( corners[(i+1)%3], corners[(i+2)%3]);
   }
   // Insert constraints for intersection objects
-  for(
-    typename list<CGAL::Object>::const_iterator lit = A_objects_3.begin();
-    lit != A_objects_3.end();
-    lit++)
+  for( const auto & obj : A_objects_3)
   {
-    CGAL::Object obj = *lit;
-    if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
-    {
-      // Add point
-      cdt.insert(P.to_2d(*ipoint));
-    } else if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
+    if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
     {
       // Add segment constraint
       cdt.insert_constraint(P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
+    }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
+    {
+      // Add point
+      cdt.insert(P.to_2d(*ipoint));
     } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj))
     {
       // Add 3 segment constraints

include/igl/cgal/remesh_self_intersections.cpp

@@ -9,6 +9,7 @@
 #include "SelfIntersectMesh.h"
 #include <igl/C_STR.h>
 #include <list>
+#include <iostream>
 
 template <
   typename DerivedV,

include/igl/colon.h

@@ -13,6 +13,12 @@ namespace igl
 {
   // Note:
   // This should be potentially replaced with eigen's LinSpaced() function
+  //
+  // If step = 1, it's about 5 times faster to use:
+  //     X = Eigen::VectorXi::LinSpaced(n,0,n-1);
+  // than 
+  //     X = igl::colon<int>(0,n-1);
+  //
 
   // Colon operator like matlab's colon operator. Enumerats values between low
   // and hi with step step.

include/igl/draw_mesh.cpp

@@ -90,7 +90,7 @@ IGL_INLINE void igl::draw_mesh(
   {
     assert(F.maxCoeff() < V.rows());
     assert(V.cols() == 3);
-    assert(rC == rV || rC == rF || rC == rF*3 || C.size() == 0);
+    assert(rC == rV || rC == rF || rC == rF*3 || rC==1 || C.size() == 0);
     assert(C.cols() == 3 || C.size() == 0);
     assert(N.cols() == 3 || N.size() == 0);
     assert(TC.cols() == 2 || TC.size() == 0);

include/igl/exterior_edges.cpp

@@ -42,11 +42,11 @@ IGL_INLINE void igl::exterior_edges(
   all_edges(F,all_E);
   long int n = F.maxCoeff()+1;
   int m = F.minCoeff();
-  const auto & compress = [&n,&m](const int i, const int j)->long int
+  const auto & compress = [n,m](const int i, const int j)->long int
   {
     return n*(i-m)+(j-m);
   };
-  const auto & decompress = [&n,&m](const long int l,int & i, int & j)
+  const auto & decompress = [n,m](const long int l,int & i, int & j)
   {
     i = (l / n) + m;
     j = (l % n) + m;

include/igl/get_modifiers.cpp

@@ -36,7 +36,7 @@ IGL_INLINE int igl::get_modifiers()
   mod |= (carbon_is_keydown(kVK_Option)?GLUT_ACTIVE_ALT:0);
   mod |= (carbon_is_keydown(kVK_Control)?GLUT_ACTIVE_CTRL:0);
 #else
-#  error "Not supported.
+#  error "Not supported."
 #endif
   return mod;
 }

include/igl/get_seconds.h

@@ -12,6 +12,21 @@
 namespace igl
 {
   // Return the current time in seconds since program start
+  // 
+  // Example:
+  //    const auto & tictoc = []()
+  //    {
+  //      static double t_start = igl::get_seconds();
+  //      double diff = igl::get_seconds()-t_start;
+  //      t_start += diff;
+  //      return diff;
+  //    };
+  //    tictoc();
+  //    ... // part 1
+  //    cout<<"part 1: "<<tictoc()<<endl;
+  //    ... // part 2
+  //    cout<<"part 2: "<<tictoc()<<endl;
+  //    ... // etc
   IGL_INLINE double get_seconds();
 
 }

include/igl/matlab/prepare_lhs.cpp

@@ -11,6 +11,17 @@ IGL_INLINE void igl::prepare_lhs_double(
   copy(&V.data()[0],&V.data()[0]+V.size(),Vp);
 }
 
+template <typename DerivedV>
+IGL_INLINE void igl::prepare_lhs_logical(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  mxArray *plhs[])
+{
+  using namespace std;
+  plhs[0] = mxCreateLogicalMatrix(V.rows(),V.cols());
+  mxLogical * Vp = static_cast<mxLogical*>(mxGetData(plhs[0]));
+  copy(&V.data()[0],&V.data()[0]+V.size(),Vp);
+}
+
 template <typename DerivedV>
 IGL_INLINE void igl::prepare_lhs_index(
   const Eigen::PlainObjectBase<DerivedV> & V,

include/igl/matlab/prepare_lhs.h

@@ -15,6 +15,11 @@ namespace igl
   IGL_INLINE void prepare_lhs_double(
     const Eigen::PlainObjectBase<DerivedV> & V,
     mxArray *plhs[]);
+  // Casts to logical
+  template <typename DerivedV>
+  IGL_INLINE void prepare_lhs_logical(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    mxArray *plhs[]);
   // Writes out a matrix and adds 1
   template <typename DerivedV>
   IGL_INLINE void prepare_lhs_index(

include/igl/per_vertex_normals.cpp

@@ -76,7 +76,8 @@ IGL_INLINE void igl::per_vertex_normals(
     // throw normal at each corner
     for(int j = 0; j < 3;j++)
     {
-      // Does this need to be critical?
+      // Q: Does this need to be critical?
+      // A: Yes. Different (i,j)'s could produce the same F(i,j)
 //#pragma omp critical
       N.row(F(i,j)) += W(i,j)*FN.row(i);
     }

include/igl/project.cpp

@@ -90,6 +90,7 @@ IGL_INLINE int igl::project(
   const Eigen::PlainObjectBase<Derivedobj> & obj,
   Eigen::PlainObjectBase<Derivedwin> & win)
 {
+  assert(obj.size() >= 3);
   Eigen::Vector3d dobj(obj(0),obj(1),obj(2));
   Eigen::Vector3d dwin;
   int ret = project(dobj(0),dobj(1),dobj(2),

include/igl/triangle_triangle_adjacency.cpp

@@ -8,7 +8,8 @@
 #include "triangle_triangle_adjacency.h"
 #include "is_edge_manifold.h"
 #include "all_edges.h"
-#include <map>
+#include "unique_simplices.h"
+#include "unique_edge_map.h"
 #include <algorithm>
 
 template <typename Scalar, typename Index>
@@ -107,59 +108,105 @@ template <
     const Eigen::PlainObjectBase<DerivedF> & F,
     std::vector<std::vector<std::vector<TTIndex> > > & TT,
     std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
+{
+  return triangle_triangle_adjacency(F,true,TT,TTi);
+}
+
+template <
+  typename DerivedF, 
+  typename TTIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT)
+{
+  std::vector<std::vector<std::vector<TTIndex> > > not_used;
+  return triangle_triangle_adjacency(F,false,TT,not_used);
+}
+
+template <
+  typename DerivedF, 
+  typename TTIndex, 
+  typename TTiIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const bool construct_TTi,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT,
+    std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
 {
   using namespace Eigen;
   using namespace std;
-  using namespace igl;
   assert(F.cols() == 3 && "Faces must be triangles");
-  typedef typename DerivedF::Index Index;
   // number of faces
   const int m = F.rows();
-  // All occurances of directed edges
-  MatrixXi E;
-  all_edges(F,E);
-  assert(E.rows() == 3*m);
-  // uE2E[i] --> {j,k,...} means unique edge i corresponds to face edges j and
-  // k (where j-edge comes is the j/m edge of face j%m)
-  map<pair<Index,Index>,vector<Index> > uE2E;
-  for(int e = 0;e<E.rows();e++)
-  {
-    Index i = E(e,0);
-    Index j = E(e,1);
-    if(i<j)
-    {
-      uE2E[pair<Index,Index>(i,j)].push_back(e);
-    }else
-    {
-      uE2E[pair<Index,Index>(j,i)].push_back(e);
-    }
-  }
+  typedef typename DerivedF::Index Index;
+  typedef Matrix<typename DerivedF::Scalar,Dynamic,2> MatrixX2I;
+  typedef Matrix<typename DerivedF::Index,Dynamic,1> VectorXI;
+  MatrixX2I E,uE;
+  VectorXI EMAP;
+  vector<vector<Index> > uE2E;
+  unique_edge_map(F,E,uE,EMAP,uE2E);
+  return triangle_triangle_adjacency(E,EMAP,uE2E,construct_TTi,TT,TTi);
+}
+
+template <
+  typename DerivedE, 
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename TTIndex, 
+  typename TTiIndex>
+  IGL_INLINE void igl::triangle_triangle_adjacency(
+    const Eigen::PlainObjectBase<DerivedE> & E,
+    const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+    const std::vector<std::vector<uE2EType> > & uE2E,
+    const bool construct_TTi,
+    std::vector<std::vector<std::vector<TTIndex> > > & TT,
+    std::vector<std::vector<std::vector<TTiIndex> > > & TTi)
+{
+  using namespace std;
+  using namespace Eigen;
+  typedef typename DerivedE::Index Index;
+  const size_t m = E.rows()/3;
+  assert(E.rows() == m*3 && "E should come from list of triangles.");
   // E2E[i] --> {j,k,...} means face edge i corresponds to other faces edges j
   // and k
-  TT.resize (m,vector<vector<TTIndex> >(F.cols()));
-  TTi.resize(m,vector<vector<TTiIndex> >(F.cols()));
-  for(int e = 0;e<E.rows();e++)
+  TT.resize (m,vector<vector<TTIndex> >(3));
+  if(construct_TTi)
+  {
+    TTi.resize(m,vector<vector<TTiIndex> >(3));
+  }
+
+  // No race conditions because TT*[f][c]'s are in bijection with e's
+  // Minimum number of iterms per openmp thread
+  const size_t ne = E.rows();
+# ifndef IGL_OMP_MIN_VALUE
+#   define IGL_OMP_MIN_VALUE 1000
+# endif
+# pragma omp parallel for if (ne>IGL_OMP_MIN_VALUE)
+  // Slightly better memory access than loop over E
+  for(Index f = 0;f<m;f++)
   {
-    const Index i = E(e,0);
-    const Index j = E(e,1);
-    const Index f = e%m;
-    const Index c = e/m;
-    const vector<Index> & N = 
-      i<j ? uE2E[pair<Index,Index>(i,j)] : uE2E[pair<Index,Index>(j,i)];
-    for(const auto & ne : N)
+    for(Index c = 0;c<3;c++)
     {
-      const Index nf = ne%m;
-      const Index nc = ne/m;
-      TT[f][c].push_back(nf);
-      TTi[f][c].push_back(nc);
+      const Index e = f + m*c;
+      const Index i = E(e,0);
+      const Index j = E(e,1);
+      //const Index f = e%m;
+      //const Index c = e/m;
+      const vector<Index> & N = uE2E[EMAP(e)];
+      for(const auto & ne : N)
+      {
+        const Index nf = ne%m;
+        TT[f][c].push_back(nf);
+        if(construct_TTi)
+        {
+          const Index nc = ne/m;
+          TTi[f][c].push_back(nc);
+        }
+      }
     }
   }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
-template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template void igl::triangle_triangle_adjacency<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/triangle_triangle_adjacency.h

@@ -76,6 +76,40 @@ namespace igl
       const Eigen::PlainObjectBase<DerivedF> & F,
       std::vector<std::vector<std::vector<TTIndex> > > & TT,
       std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
+  template < typename DerivedF, typename TTIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT);
+  // Wrapper with bool to choose whether to compute TTi (this prototype should
+  // be "hidden").
+  template <
+    typename DerivedF, 
+    typename TTIndex, 
+    typename TTiIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      const bool construct_TTi,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT,
+      std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
+  // Inputs:
+  //   E  #F*3 by 2 list of all of directed edges in order (see `all_edges`)
+  //   EMAP #F*3 list of indices into uE, mapping each directed edge to unique
+  //     undirected edge
+  //   uE2E  #uE list of lists of indices into E of coexisting edges
+  // See also: unique_edge_map, all_edges
+  template <
+    typename DerivedE, 
+    typename DerivedEMAP,
+    typename uE2EType,
+    typename TTIndex, 
+    typename TTiIndex>
+    IGL_INLINE void triangle_triangle_adjacency(
+      const Eigen::PlainObjectBase<DerivedE> & E,
+      const Eigen::PlainObjectBase<DerivedEMAP> & EMAP,
+      const std::vector<std::vector<uE2EType > > & uE2E,
+      const bool construct_TTi,
+      std::vector<std::vector<std::vector<TTIndex> > > & TT,
+      std::vector<std::vector<std::vector<TTiIndex> > > & TTi);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/unique_simplices.cpp

@@ -24,7 +24,7 @@ IGL_INLINE void igl::unique_simplices(
   using namespace igl;
   // Sort each face
   MatrixXi sortF, unusedI;
-  igl::sort(F,2,1,sortF,unusedI);
+  igl::sort(F,2,true,sortF,unusedI);
   // Find unique faces
   MatrixXi C;
   igl::unique_rows(sortF,C,IA,IC);