Dynamic updates to AABB tree; intersection-blocking mesh decimation (#2301)

* working insertion and rotation b-b-b-but pointers aren't stuck to m_primitive

* insert now maintains primitive's pointers;+rotate is now getting heights in right ballpark

* sibling rotations working (and helping); dry run test working (and helping)

* working insertion, deletion (detach), and refit with padding.

* working; inexact

* working; inexact

* working; inexact

* note about poor assumptions

* working well after bug fixes. before refactor into functions

* simple self-intersection test function

* moved all functions to files

* docs

* blocking directly in qslim. better docs. aabb templates/tests;

* dont use size_t and fix namespace

* brute-force too fast on linux

* rm overloads in qslim/decimate; cgal template

* fix coplanar bug; factor out raytri.c

* fix and cgal debug

* refactor fast_find; fix bugs in fast_find; fix bugs in shared_vertex

* tutorials running again

* cleaned up aabb tutorials; templates

* format docs

* docs. arg names

* template name

* improve docs

* debugging test

* debugging test

* debugging test

* debugging test

* debugging test

* debugging test

* debugging test

* debugging test

* ebuggin test

* ebuggin test

* add epsilon to ray_triangle ifs

* erroneous includes

* rm leftover includes

* fix cmake bug

* uh actually fix cmake bug

* missing delete

* simple insert test

* don't pad all leaves F.rows() times

* fix pad bug

---------

Co-authored-by: Alec Jacobson <[email protected]>

cmake/igl/igl_copy_dll.cmake

@@ -50,6 +50,10 @@ function(igl_copy_dll target)
     if(NOT WIN32)
         return()
     endif()
+    if(NOT TARGET target)
+      message(STATUS "igl_copy_dll() was called with a non-target: ${target}")
+      return()
+    endif()
 
     # Sanity checks
     get_target_property(TYPE ${target} TYPE)

include/igl/AABB.h

@@ -29,6 +29,9 @@ namespace igl
     class AABB 
     {
 public:
+///////////////////////////////////////////////////////////////////////////////
+// Member variables
+///////////////////////////////////////////////////////////////////////////////
       /// Scalar type of vertex positions (e.g., `double`)
       typedef typename DerivedV::Scalar Scalar;
       /// Fixed-size (`DIM`) RowVector type using `Scalar`
@@ -42,26 +45,32 @@ public:
       AABB * m_left; 
       /// Pointer to "right" child node (`nullptr` if leaf)
       AABB * m_right;
+      /// Pointer to "parent" node (`nullptr` if root)
+      AABB * m_parent;
       /// Axis-Aligned Bounding Box containing this node
       Eigen::AlignedBox<Scalar,DIM> m_box;
       /// Index of single primitive in this node if full leaf, otherwise -1 for non-leaf
       int m_primitive;
+///////////////////////////////////////////////////////////////////////////////
+// Non-templated member functions
+///////////////////////////////////////////////////////////////////////////////
       //Scalar m_low_sqr_d;
       //int m_depth;
       /// @private
       AABB():
-        m_left(NULL), m_right(NULL),
+        m_left(nullptr), m_right(nullptr),m_parent(nullptr),
         m_box(), m_primitive(-1)
         //m_low_sqr_d(std::numeric_limits<double>::infinity()),
         //m_depth(0)
-    {
-      static_assert(DerivedV::ColsAtCompileTime == DIM || DerivedV::ColsAtCompileTime == Eigen::Dynamic,"DerivedV::ColsAtCompileTime == DIM || DerivedV::ColsAtCompileTime == Eigen::Dynamic");
-    }
+      {
+        static_assert(DerivedV::ColsAtCompileTime == DIM || DerivedV::ColsAtCompileTime == Eigen::Dynamic,"DerivedV::ColsAtCompileTime == DIM || DerivedV::ColsAtCompileTime == Eigen::Dynamic");
+      }
       /// @private
       // http://stackoverflow.com/a/3279550/148668
       AABB(const AABB& other):
-        m_left(other.m_left ? new AABB(*other.m_left) : NULL),
-        m_right(other.m_right ? new AABB(*other.m_right) : NULL),
+        m_left (other.m_left  ? new AABB(*other.m_left)  : nullptr),
+        m_right(other.m_right ? new AABB(*other.m_right) : nullptr),
+        m_parent(other.m_parent),
         m_box(other.m_box),
         m_primitive(other.m_primitive)
         //m_low_sqr_d(other.m_low_sqr_d),
@@ -69,6 +78,8 @@ public:
         //   m_left ? m_left->m_depth + 1 : 0,
         //   m_right ? m_right->m_depth + 1 : 0))
         {
+          if(m_left) { m_left->m_parent = this; }
+          if(m_right) { m_right->m_parent = this; }
         }
       /// @private
       // copy-swap idiom
@@ -78,6 +89,7 @@ public:
         using std::swap;
         swap(first.m_left,second.m_left);
         swap(first.m_right,second.m_right);
+        swap(first.m_parent,second.m_parent);
         swap(first.m_box,second.m_box);
         swap(first.m_primitive,second.m_primitive);
         //swap(first.m_low_sqr_d,second.m_low_sqr_d);
@@ -100,28 +112,365 @@ public:
       /// @private
       // Seems like there should have been an elegant solution to this using
       // the copy-swap idiom above:
-      IGL_INLINE void deinit()
+      IGL_INLINE void clear()
       {
         m_primitive = -1;
         m_box = Eigen::AlignedBox<Scalar,DIM>();
         delete m_left;
-        m_left = NULL;
+        m_left = nullptr;
         delete m_right;
-        m_right = NULL;
+        m_right = nullptr;
+        // Tell my parent I'm dead
+        if(m_parent)
+        {
+          if(m_parent->m_left == this)
+          {
+            m_parent->m_left = nullptr;
+          }else if(m_parent->m_right == this)
+          {
+            m_parent->m_right = nullptr;
+          }else
+          {
+            assert(false && "I'm not my parent's child");
+          }
+          auto * grandparent = m_parent->m_parent;
+          if(grandparent)
+          {
+            // Before
+            //     grandparent
+            //        /    \
+            //    parent   pibling
+            //      /   \
+            // sibling  this
+            //
+            // After 
+            //     grandparent
+            //        /    \
+            //    sibling®  pibling
+          }else
+          {
+            // Before
+            //    parent=root
+            //      /   \
+            // sibling  this
+            //
+            // After 
+            //     grandparent
+            //        /    \
+            //    sibling®  pibling
+          }
+        }
+        // Now my parent is dead to me.
+        m_parent = nullptr;
       }
       /// @private
       ~AABB()
       {
-        deinit();
+        clear();
       }
+      /// Return whether at leaf node
+      IGL_INLINE bool is_leaf() const;
+      /// Return whether at root node
+      IGL_INLINE bool is_root() const;
+      /// Return the root node of this node's tree by following its parent
+      IGL_INLINE AABB<DerivedV,DIM>* root() const;
+      IGL_INLINE AABB<DerivedV,DIM>* detach();
+      IGL_INLINE void refit_lineage();
+      /// Get a vector of leaves indexed by their m_primitive id (these better
+      /// be non-negative and tightly packed.
+      /// @param[in] m  number of leaves/elements (Ele.rows())
+      /// @returns leaves  m list of pointers to leaves
+      IGL_INLINE std::vector<AABB<DerivedV,DIM>*> gather_leaves(const int m);
+      /// \overload where m is the max m_primitive id in the tree.
+      IGL_INLINE std::vector<AABB<DerivedV,DIM>*> gather_leaves();
+      /// Pad leaves by `pad` in each dimension
+      /// @param[in] pad padding amount
+      /// @param[in] polish_rotate_passes number of passes to polish rotations
+      /// @returns pointer to (potentially new) root
+      IGL_INLINE AABB<DerivedV,DIM>* pad(
+        const std::vector<AABB<DerivedV,DIM>*> & leaves,
+        const Scalar pad, 
+        const int polish_rotate_passes=0);
+      /// @returns `this` if no update was needed, otherwise returns pointer to
+      /// (potentially new) root
+      ///
+      /// Example:
+      /// ```cpp
+      /// auto * up = leaf->update(new_box);
+      /// if(up != leaf)
+      /// {
+      ///   tree = up->root();
+      /// }else
+      /// {
+      ///   printf("no update occurred\n");
+      /// }
+      ///
+      /// // or simply
+      /// tree = leaf->update(new_box)->root();
+      /// ```
+      IGL_INLINE AABB<DerivedV,DIM>* update(
+          const Eigen::AlignedBox<Scalar,DIM> & new_box,
+          const Scalar pad=0);
+      /// Insert a (probably a leaf) AABB `other` into this AABB tree. If
+      /// `other`'s box is contained in this AABB's box then insert it as a child recursively.
+      ///
+      /// If `other`'s box is not contained in this AABB's box then insert it as a
+      /// sibling. 
+      ///
+      /// It's a very good idea to call either `rotate` (faster, less good) or `rotate_lineage` (slower, better)
+      /// after insertion. Rotating continues to improve the tree's quality so
+      /// after doing a bunch of insertions you might even consider calling
+      /// `rotate` on all nodes.
+      ///
+      /// `insert` attempts to minimize total internal surface area. Where as
+      /// `init` is top-down and splits boxes based on the median along the
+      /// longest dimension. When initializing a tree, `init` seems to result in
+      /// great trees (small height and small total internal surface area).
+      ///
+      /// @param[in] other pointer to another AABB node 
+      /// @returns pointer to the parent of `other` or `other` itself. This
+      /// could be == to a `new`ly created internal node or to `other` if
+      /// `this==other`. Calling ->root() on this returned node will give you
+      /// the root of the tree.
+      ///
+      /// ##### Example
+      /// 
+      /// ```cpp
+      /// // Create a tree (use pointer to track changes to root)
+      /// auto * tree = new igl::AABB<DerivedV,3>::AABB();
+      /// // Fill the tree (e.g., using ->init())
+      /// …
+      /// // Create a new leafe node
+      /// auto * leaf = new igl::AABB<DerivedV,3>::AABB();
+      /// // Fill the leaf node with a primitive and box
+      /// …
+      /// // Insert into the tree and find the possibly new root
+      /// tree = tree->insert(leaf)->root();
+      /// ```
+      IGL_INLINE AABB<DerivedV,DIM>* insert(AABB * other);
+      /// Insert `other` as a sibling to `this` by creating a new internal node
+      /// to be their shared parent.
+      ///
+      ///     Before
+      ///              parent
+      ///              /    \
+      ///          this(C)  sibling
+      ///            /  \
+      ///          left right
+      ///    
+      ///     After
+      ///              parent
+      ///              /    \
+      ///           newbie   sibling
+      ///            /   \
+      ///         this    other
+      ///         /    \
+      ///       left  right
+      ///    
+      ///
+      /// @param[in] other pointer to another AABB node 
+      /// @returns pointer to the new shared parent.
+      IGL_INLINE AABB<DerivedV,DIM>* insert_as_sibling(AABB * other);
+      /// Try to swap this node with its close relatives if it will decrease
+      /// total internal surface area.
+      ///    
+      ///
+      ///            grandparent
+      ///            /         \
+      ///         parent       pibling°
+      ///         /    \         /    \
+      ///     sibling  this   cuz1°  cuz2°
+      ///      /    \
+      ///     nib1° nib2°
+      ///
+      ///     °Swap Candidates
+      ///
+      /// @param[in] dry_run  if true then don't actually swap
+      /// @return[in] the change in total internal surface area, 0 if no
+      /// improvement and rotate won't be carried out.
+      IGL_INLINE Scalar rotate(const bool dry_run = false);
+      /// Try to swap this node with its cousins if it will decrease
+      /// total internal surface area.
+      ///    
+      /// @param[in] dry_run  if true then don't actually swap
+      /// @return[in] the change in total internal surface area, 0 if no
+      /// improvement and rotate won't be carried out.
+      ///    
+      ///     Before
+      ///            grandparent
+      ///            /         \
+      ///         parent       pibling
+      ///         /    \         /    \
+      ///     sibling  this    cuz1  cuz2
+      ///    
+      ///    
+      ///     Candidates
+      ///            grandparent
+      ///            /         \
+      ///         parent       pibling
+      ///         /    \         /    \
+      ///     sibling  cuz1   this   cuz2
+      ///    
+      ///     Or
+      ///            grandparent
+      ///            /         \
+      ///         parent       pibling
+      ///         /    \         /    \
+      ///     sibling  cuz2    cuz1  this
+      IGL_INLINE Scalar rotate_across(const bool dry_run = false);
+      /// Try to swap this node with its pibling if it will decrease
+      /// total internal surface area.
+      ///    
+      /// @param[in] dry_run  if true then don't actually swap
+      /// @return[in] the change in total internal surface area, 0 if no
+      /// improvement and rotate won't be carried out.
+      ///
+      ///     Before
+      ///        grandparent
+      ///           /    \
+      ///         other  parent
+      ///                /  \
+      ///             this  sibling
+      ///             
+      ///    
+      ///     Candidate
+      ///        grandparent
+      ///           /    \
+      ///        this    parent
+      ///                /  \
+      ///            other  sibling
+      IGL_INLINE Scalar rotate_up(const bool dry_run = false);
+      /// Try to swap this node with one of its niblings if it will decrease
+      /// total internal surface area.
+      ///
+      /// @param[in] dry_run  if true then don't actually swap
+      /// @return[in] the change in total internal surface area, 0 if no
+      /// improvement and rotate won't be carried out.
+      /// 
+      ///     Before
+      ///           parent
+      ///           /    \
+      ///         this   sibling
+      ///                /  \
+      ///             left  right
+      ///             
+      ///
+      ///     Candidates
+      ///           parent
+      ///           /    \
+      ///       left     sibling
+      ///                /  \
+      ///            this   right
+      ///
+      ///     Or
+      ///
+      ///           parent
+      ///           /    \
+      ///       right    sibling
+      ///                /  \
+      ///            left   this 
+      IGL_INLINE Scalar rotate_down(const bool dry_run = false);
+      /// "Rotate" (swap) `reining` with `challenger`.
+      ///
+      ///     Before
+      ///        grandparent
+      ///           /      \
+      ///      reining      parent
+      ///                   /    \
+      ///            challenger  sibling
+      ///             
+      ///    
+      ///     Candidate
+      ///        grandparent
+      ///           /      \
+      ///     challenger    parent
+      ///                   /    \
+      ///              reining   sibling
+      /// @param[in] reining  pointer to AABB node to be rotated
+      /// @param[in] grandparent  pointer to challenger's grandparent
+      /// @param[in] parent  pointer to challenger's parent
+      /// @param[in] challenger  pointer to AABB node to be rotated
+      /// @param[in] sibling  pointer to challenger's sibling
+      /// @returns true only if rotation was possible and successfully carried
+      /// out.
+      static IGL_INLINE Scalar rotate_up(
+        const bool dry_run,
+        AABB<DerivedV,DIM>* reining,
+        AABB<DerivedV,DIM>* grandparent,
+        AABB<DerivedV,DIM>* parent,
+        AABB<DerivedV,DIM>* challenger,
+        AABB<DerivedV,DIM>* sibling);
+      // Should this be a static function with an argument?
+      IGL_INLINE void rotate_lineage();
+      /// Number of nodes contained in subtree (is it?)
+      ///
+      /// \note At best, this function has a dubious name. This is really an
+      /// internal helper function for the serialization.
+      ///
+      /// \see size()
+      ///
+      /// @return Number of elements m then total tree size should be 2*h where h is
+      /// the deepest depth 2^ceil(log(#Ele*2-1))
+      IGL_INLINE int subtree_size() const;
+      /// @param[in]  box  query box
+      /// @param[in,out] leaves  list of leaves to append to
+      IGL_INLINE bool append_intersecting_leaves(
+        const Eigen::AlignedBox<Scalar,DIM> & box,
+        std::vector<const AABB<DerivedV,DIM>*> & leaves) const;
+      /// Compute sum of surface area of all internal (non-root, non-leaf) boxes
+      IGL_INLINE typename DerivedV::Scalar internal_surface_area() const;
+      /// Validate the subtree under this node by running a bunch of assertions.
+      /// Does nothing when not in debug mode
+      IGL_INLINE void validate() const;
+      /// print the memory addresses of the tree in a somewhat legible way
+      IGL_INLINE void print(const int depth = 0) const;
+      /// @returns Actual size of tree. Total number of nodes in tree. A singleton root
+      /// has size 1.
+      /// 
+      /// \see subtree_size
+      IGL_INLINE int size() const;
+      /// @returns Height of the tree. A singleton root has height 1.
+      IGL_INLINE int height() const;
+private:
+      /// If new distance (sqr_d_candidate) is less than current distance
+      /// (sqr_d), then update this distance and its associated values
+      /// _in-place_:
+      ///
+      /// @param[in] p  dim-long query point (only used in DEBUG mode)
+      /// @param[in] sqr_d  candidate minimum distance for this query, see
+      ///   output
+      /// @param[in] i  candidate index into Ele of closest point, see output
+      /// @param[in] c  dim-long candidate closest point, see output
+      /// @param[in] sqr_d  current minimum distance for this query, see output
+      /// @param[in] i  current index into Ele of closest point, see output
+      /// @param[in] c  dim-long current closest point, see output
+      /// @param[out] sqr_d   minimum of initial value and squared distance to
+      ///   this primitive
+      /// @param[out] i  possibly updated index into Ele of closest point
+      /// @param[out] c  dim-long possibly updated closest point
+      IGL_INLINE void set_min(
+        const RowVectorDIMS & p,
+        const Scalar sqr_d_candidate,
+        const int i_candidate,
+        const RowVectorDIMS & c_candidate,
+        Scalar & sqr_d,
+        int & i,
+        Eigen::PlainObjectBase<RowVectorDIMS> & c) const;
+public:
+///////////////////////////////////////////////////////////////////////////////
+// Templated member functions
+///////////////////////////////////////////////////////////////////////////////
       /// Build an Axis-Aligned Bounding Box tree for a given mesh and given
       /// serialization of a previous AABB tree.
       ///
       /// @param[in] V  #V by dim list of mesh vertex positions. 
       /// @param[in] Ele  #Ele by dim+1 list of mesh indices into #V. 
-      /// @param[in] bb_mins  max_tree by dim list of bounding box min corner positions
-      /// @param[in] bb_maxs  max_tree by dim list of bounding box max corner positions
-      /// @param[in] elements  max_tree list of element or (not leaf id) indices into Ele
+      /// @param[in] bb_mins  max_tree by dim list of bounding box min corner
+      ///   positions
+      /// @param[in] bb_maxs  max_tree by dim list of bounding box max corner
+      ///   positions
+      /// @param[in] elements  max_tree list of element or (not leaf id) indices
+      ///   into Ele
       /// @param[in] i  recursive call index {0}
       template <
         typename DerivedEle, 
@@ -148,12 +497,12 @@ public:
       ///
       /// @param[in] V  #V by dim list of mesh vertex positions. 
       /// @param[in] Ele  #Ele by dim+1 list of mesh indices into #V. 
-      /// @param[in] SI  #Ele by dim list revealing for each coordinate where Ele's
-      ///              barycenters would be sorted: SI(e,d) = i --> the dth coordinate of
-      ///              the barycenter of the eth element would be placed at position i in a
-      ///              sorted list.
-      /// @param[in] I  #I list of indices into Ele of elements to include (for recursive
-      ///     calls)
+      /// @param[in] SI  #Ele by dim list revealing for each coordinate where
+      ///   Ele's barycenters would be sorted: SI(e,d) = i --> the dth
+      ///   coordinate of the barycenter of the eth element would be placed at
+      ///   position i in a sorted list.
+      /// @param[in] I  #I list of indices into Ele of elements to include (for
+      ///   recursive calls)
       /// 
       template <typename DerivedEle, typename DerivedSI, typename DerivedI>
       IGL_INLINE void init(
@@ -161,15 +510,21 @@ public:
           const Eigen::MatrixBase<DerivedEle> & Ele, 
           const Eigen::MatrixBase<DerivedSI> & SI,
           const Eigen::MatrixBase<DerivedI>& I);
-      /// Return whether at leaf node
-      IGL_INLINE bool is_leaf() const;
+      /// @returns `this` if no update was needed, otherwise returns pointer to
+      /// (potentially new) root
+      template <typename DerivedEle>
+      IGL_INLINE AABB<DerivedV,DIM>* update_primitive(
+          const Eigen::MatrixBase<DerivedV> & V,
+          const Eigen::MatrixBase<DerivedEle> & Ele,
+          const Scalar pad=0);
+
       /// Find the indices of elements containing given point: this makes sense
       /// when Ele is a co-dimension 0 simplex (tets in 3D, triangles in 2D).
       ///
-      /// @param[in]  V  #V by dim list of mesh vertex positions. **Should be same as used to
-      ///     construct mesh.**
-      /// @param[in]  Ele  #Ele by dim+1 list of mesh indices into #V. **Should be same as used to
-      ///     construct mesh.**
+      /// @param[in]  V  #V by dim list of mesh vertex positions. **Should be
+      ///   same as used to construct mesh.**
+      /// @param[in]  Ele  #Ele by dim+1 list of mesh indices into #V. **Should
+      ///   be same as used to construct mesh.**
       /// @param[in]  q  dim row-vector query position
       /// @param[in]  first  whether to only return first element containing q
       /// @return  list of indices of elements containing q
@@ -180,17 +535,14 @@ public:
           const Eigen::MatrixBase<Derivedq> & q,
           const bool first=false) const;
 
-      /// Number of nodes contained in subtree
-      ///
-      /// @return Number of elements m then total tree size should be 2*h where h is
-      /// the deepest depth 2^ceil(log(#Ele*2-1))
-      IGL_INLINE int subtree_size() const;
-
       /// Serialize this class into 3 arrays (so we can pass it pack to matlab)
       ///
-      /// @param[out]  bb_mins  max_tree by dim list of bounding box min corner positions
-      /// @param[out]  bb_maxs  max_tree by dim list of bounding box max corner positions
-      /// @param[out]  elements  max_tree list of element or (not leaf id) indices into Ele
+      /// @param[out]  bb_mins  max_tree by dim list of bounding box min corner
+      ///   positions
+      /// @param[out]  bb_maxs  max_tree by dim list of bounding box max corner
+      ///   positions
+      /// @param[out]  elements  max_tree list of element or (not leaf id)
+      ///   indices into Ele
       /// @param[in]  i  recursive call index into these arrays {0}
       template <
         typename Derivedbb_mins, 
@@ -225,11 +577,11 @@ public:
       /// @param[in] V  #V by dim list of vertex positions
       /// @param[in] Ele  #Ele by dim list of simplex indices
       /// @param[in] p  dim-long query point 
-      /// @param[in] low_sqr_d  lower bound on squared distance, specified maximum squared
-      ///              distance 
-      /// @param[in] up_sqr_d  current upper bounded on squared distance, current minimum
-      ///              squared distance (only consider distances less than this), see
-      ///              output.
+      /// @param[in] low_sqr_d  lower bound on squared distance, specified
+      ///   maximum squared distance 
+      /// @param[in] up_sqr_d  current upper bounded on squared distance,
+      ///   current minimum squared distance (only consider distances less than
+      ///   this), see output.
       /// @param[out]  i  facet index corresponding to smallest distances
       /// @param[out]  c  closest point
       /// @return squared distance
@@ -250,9 +602,9 @@ public:
       /// @param[in] V  #V by dim list of vertex positions
       /// @param[in] Ele  #Ele by dim list of simplex indices
       /// @param[in] p  dim-long query point 
-      /// @param[in] up_sqr_d  current upper bounded on squared distance, current minimum
-      ///              squared distance (only consider distances less than this), see
-      ///              output.
+      /// @param[in] up_sqr_d  current upper bounded on squared distance,
+      ///   current minimum squared distance (only consider distances less than
+      ///   this), see output.
       /// @param[out]  i  facet index corresponding to smallest distances
       /// @param[out]  c  closest point
       /// @return squared distance
@@ -312,8 +664,6 @@ public:
         const RowVectorDIMS & dir,
         const Scalar min_t,
         igl::Hit & hit) const;
-
-public:
       /// Compute the squared distance from all query points in P to the
       /// _closest_ points on the primitives stored in the AABB hierarchy for
       /// the mesh (V,Ele).
@@ -337,7 +687,6 @@ public:
         Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
         Eigen::PlainObjectBase<DerivedI> & I,
         Eigen::PlainObjectBase<DerivedC> & C) const;
-
       /// Compute the squared distance from all query points in P already stored
       /// in its own AABB hierarchy to the _closest_ points on the primitives
       /// stored in the AABB hierarchy for the mesh (V,Ele).
@@ -418,32 +767,6 @@ private:
         Scalar & sqr_d,
         int & i,
         Eigen::PlainObjectBase<RowVectorDIMS> & c) const;
-      // If new distance (sqr_d_candidate) is less than current distance
-      // (sqr_d), then update this distance and its associated values
-      // _in-place_:
-      //
-      // Inputs:
-      //   p  dim-long query point (only used in DEBUG mode)
-      //   sqr_d  candidate minimum distance for this query, see output
-      //   i  candidate index into Ele of closest point, see output
-      //   c  dim-long candidate closest point, see output
-      //   sqr_d  current minimum distance for this query, see output
-      //   i  current index into Ele of closest point, see output
-      //   c  dim-long current closest point, see output
-      // Outputs:
-      //   sqr_d   minimum of initial value and squared distance to this
-      //     primitive
-      //   i  possibly updated index into Ele of closest point
-      //   c  dim-long possibly updated closest point
-      IGL_INLINE void set_min(
-        const RowVectorDIMS & p,
-        const Scalar sqr_d_candidate,
-        const int i_candidate,
-        const RowVectorDIMS & c_candidate,
-        Scalar & sqr_d,
-        int & i,
-        Eigen::PlainObjectBase<RowVectorDIMS> & c) const;
-
       /// Intersect a ray with the mesh return all hits
       ///
       /// @param[in]  V  #V by dim list of vertex positions

include/igl/box_faces.cpp

@@ -0,0 +1,94 @@
+#include "box_faces.h"
+#include "AABB.h"
+#include <vector>
+#include <utility>
+
+template <typename DerivedV, typename DerivedQ>
+IGL_INLINE void igl::box_faces(
+  const Eigen::AlignedBox<typename DerivedV::Scalar,3> & box,
+  const typename DerivedV::Scalar shrink,
+  Eigen::PlainObjectBase<DerivedV> & P,
+  Eigen::PlainObjectBase<DerivedQ> & Q)
+{
+  auto min_corner = box.min();
+  auto max_corner = box.max();
+  // shrink by 3%
+  min_corner = min_corner + shrink*(max_corner-min_corner);
+  max_corner = max_corner - shrink*(max_corner-min_corner);
+  P.resize(8,3);
+  Q.resize(6,4);
+  int p = 0;
+  int q = 0;
+  Q.row(q++) << p+0,p+1,p+2,p+3;
+  Q.row(q++) << p+0,p+1,p+5,p+4;
+  Q.row(q++) << p+1,p+2,p+6,p+5;
+  Q.row(q++) << p+2,p+3,p+7,p+6;
+  Q.row(q++) << p+3,p+0,p+4,p+7;
+  Q.row(q++) << p+4,p+5,p+6,p+7;
+  P.row(p++) = min_corner;
+  P.row(p++) = Eigen::RowVector3d(max_corner[0],min_corner[1],min_corner[2]);
+  P.row(p++) = Eigen::RowVector3d(max_corner[0],max_corner[1],min_corner[2]);
+  P.row(p++) = Eigen::RowVector3d(min_corner[0],max_corner[1],min_corner[2]);
+  P.row(p++) = Eigen::RowVector3d(min_corner[0],min_corner[1],max_corner[2]);
+  P.row(p++) = Eigen::RowVector3d(max_corner[0],min_corner[1],max_corner[2]);
+  P.row(p++) = max_corner;
+  P.row(p++) = Eigen::RowVector3d(min_corner[0],max_corner[1],max_corner[2]);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedP,
+  typename DerivedQ,
+  typename DerivedD >
+IGL_INLINE void igl::box_faces(
+  const igl::AABB<DerivedV,3> & tree,
+  Eigen::PlainObjectBase<DerivedP> & P,
+  Eigen::PlainObjectBase<DerivedQ> & Q,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  const int num_nodes = tree.size();
+  P.resize(8*num_nodes,3);
+  Q.resize(6*num_nodes,4);
+  D.resize(6*num_nodes);
+  int d = 0;
+  int p = 0;
+  int q = 0;
+  std::vector<std::pair<const igl::AABB<DerivedV,3> *,int> > stack;
+  stack.push_back({&tree,0});
+  while(!stack.empty())
+  {
+    const auto pair = stack.back();
+    const auto * node = pair.first;
+    const int depth = pair.second;
+    D(d++) = depth;
+    D(d++) = depth;
+    D(d++) = depth;
+    D(d++) = depth;
+    D(d++) = depth;
+    D(d++) = depth;
+    stack.pop_back();
+    const auto & box = node->m_box;
+
+    DerivedP Pi;
+    DerivedQ Qi;
+    box_faces(box,0.03,Pi,Qi);
+    P.block(p,0,8,3) = Pi;
+    Q.block(q,0,6,4) = Qi.array()+p;
+    p += 8;
+    q += 6;
+    if(node->m_left)
+    {
+      stack.push_back({node->m_left,depth+1});
+    }
+    if(node->m_right)
+    {
+      stack.push_back({node->m_right,depth+1});
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::box_faces<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::Matrix<int, -1, 1, 0, -1, 1> >(igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/box_faces.h

@@ -0,0 +1,46 @@
+#ifndef IGL_BOX_FACES_H
+#define IGL_BOX_FACES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+namespace igl
+{
+  /// Compute the quad faces of an axis-aligned bounding box) shrunk by a given
+  /// factor.
+  ///
+  /// @param[in] box Axis-aligned bounding box
+  /// @param[in] shrink Factor by which to shrink the box
+  /// @param[out] P  #P by 3 list of vertex positions
+  /// @param[out] Q  #Q by 4 list of triangle indices into rows of P
+  template <typename DerivedV, typename DerivedQ>
+  IGL_INLINE void box_faces(
+    const Eigen::AlignedBox<typename DerivedV::Scalar,3> & box,
+    const typename DerivedV::Scalar shrink,
+    Eigen::PlainObjectBase<DerivedV> & P,
+    Eigen::PlainObjectBase<DerivedQ> & Q);
+  // Forward declaration
+  template <typename DerivedV, int DIM> class AABB;
+  /// Compute the quad faces of a tree of axis-aligned bounding boxes.
+  ///
+  /// @param[in] tree Tree of axis-aligned bounding boxes
+  /// @param[out] P  #P by 3 list of vertex positions
+  /// @param[out] Q  #Q by 4 list of triangle indices into rows of P
+  /// @param[out] D  #Q list of tree depths (0==root)
+  template <
+    typename DerivedV,
+    typename DerivedP,
+    typename DerivedQ,
+    typename DerivedD >
+  IGL_INLINE void box_faces(
+    const igl::AABB<DerivedV,3> & tree,
+    Eigen::PlainObjectBase<DerivedP> & P,
+    Eigen::PlainObjectBase<DerivedQ> & Q,
+    Eigen::PlainObjectBase<DerivedD> & D);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "box_faces.cpp"
+#endif
+
+#endif

include/igl/box_surface_area.cpp

@@ -0,0 +1,35 @@
+#include "box_surface_area.h"
+
+
+template <typename DerivedCorner>
+IGL_INLINE typename DerivedCorner::Scalar igl::box_surface_area(
+  const Eigen::MatrixBase<DerivedCorner> & min_corner,
+  const Eigen::MatrixBase<DerivedCorner> & max_corner)
+{
+  using Scalar = typename DerivedCorner::Scalar;
+  const auto dimensions = (max_corner - min_corner).eval();
+  const auto num_dimensions = dimensions.size();
+    
+  Scalar surface_area = 0;
+  for (int i = 0; i < num_dimensions; ++i) {
+      for (int j = i + 1; j < num_dimensions; ++j) {
+          surface_area += 2 * dimensions[i] * dimensions[j];
+      }
+  }
+    
+  return surface_area;
+}
+
+template <typename Scalar, int AmbientDim>
+IGL_INLINE Scalar igl::box_surface_area(
+  const Eigen::AlignedBox<Scalar,AmbientDim> & box)
+{
+  return igl::box_surface_area(box.min(),box.max());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template double igl::box_surface_area<double, 2>(Eigen::AlignedBox<double, 2> const&);
+template double igl::box_surface_area<double, 3>(Eigen::AlignedBox<double, 3> const&);
+#endif

include/igl/box_surface_area.h

@@ -0,0 +1,30 @@
+#ifndef IGL_BOX_SURFACE_AREA_H
+#define IGL_BOX_SURFACE_AREA_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+namespace igl
+{
+  /// Compute the surface area of a box given its min and max corners
+  /// 
+  /// @param[in] min_corner  #d vector of min corner position
+  /// @param[in] max_corner  #d vector of max corner position
+  /// @return            surface area of box
+  template <typename DerivedCorner>
+  IGL_INLINE typename DerivedCorner::Scalar box_surface_area(
+    const Eigen::MatrixBase<DerivedCorner> & min_corner,
+    const Eigen::MatrixBase<DerivedCorner> & max_corner);
+  /// \overload
+  /// @param[in] box  axis-aligned bounding box
+  template <typename Scalar, int AmbientDim>
+  IGL_INLINE Scalar box_surface_area(
+    const Eigen::AlignedBox<Scalar,AmbientDim> & box);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "box_surface_area.cpp"
+#endif
+
+#endif
+

include/igl/collapse_edge.cpp

@@ -82,10 +82,10 @@ IGL_INLINE bool igl::collapse_edge(
   {
     E(e,0) = IGL_COLLAPSE_EDGE_NULL;
     E(e,1) = IGL_COLLAPSE_EDGE_NULL;
-    EF(e,0) = IGL_COLLAPSE_EDGE_NULL;
-    EF(e,1) = IGL_COLLAPSE_EDGE_NULL;
-    EI(e,0) = IGL_COLLAPSE_EDGE_NULL;
-    EI(e,1) = IGL_COLLAPSE_EDGE_NULL;
+    // Don't clear EF, EI in case post_collapse would like to access previous
+    // connectivity. It's questionable whether E should be cleared at all, but
+    // if the user is drawing things with E then it's convenient that it's
+    // mapped to "null" edges (similar to F).
   };
 
   // update edge info

include/igl/collapse_edge.h

@@ -48,6 +48,66 @@ namespace igl
   /// @param[out] f1  index into F of face collpased on left
   /// @param[out] f2  index into F of face collpased on right
   /// @return true if edge was collapsed
+  ///
+  ///
+  /// Define [s,d] = sort(E(e,:)) so that s<d, then d is "detached" from
+  /// connectivity meaning all faces/edges incident on d will now be incident on
+  /// s. (This reduces fragmentation by preferring to collapse toward the start
+  /// of V)¹. If E(e,1)==s then we say the edge is "flipped" (`eflip` true in
+  /// the implementation).
+  ///
+  /// f1 is set to EF(e,0) and f2 is set to EF(e,1). Let v1 be EI(e,0) the
+  /// corner of F(f1,:) opposite e. _If_ (s<d) then e1 will be the edge after e
+  /// within f1:
+  ///
+  ///         s<d
+  ///     ✅s----e-----d☠️
+  ///        \   ←    /
+  ///         \ ↘f₁↗ /
+  ///         e₁    /
+  ///           \  /
+  ///            \/
+  ///
+  /// _If_ (s>d) then e1 will be the edge after e within f1:
+  ///    
+  ///         s>d
+  ///     ✅s----e-----d☠️
+  ///        \   ←    /
+  ///         \ ↘f₁↗ /
+  ///          \    e₁
+  ///           \  /
+  ///            \/
+  ///
+  ///
+  /// ¹Or at least it would if we templated these functions to allow using
+  /// RowMajor V.
+  ///
+  /// It really seems that this callback should provide a meaningful edge on the
+  /// _new_ mesh. Meanwhile – Oof – You can use this gross mechanism to find the faces incident on the
+  /// collapsed vertex:
+  ///
+  /// ```cpp
+  /// const auto survivors = 
+  ///   [&F,&e,&EMAP](const int f1, const int e1, int & d1)
+  /// {
+  ///   for(int c=0;c<3;c++)
+  ///   {
+  ///     d1 = EMAP(f1+c*F.rows());
+  ///     if((d1 != e) && (d1 != e1)) { break; }
+  ///   }
+  /// };
+  /// int d1,d2;
+  /// survivors(f1,e1,d1);
+  /// survivors(f2,e2,d2);
+  /// // Will circulating by continuing in the CCW direction of E(d1,:)
+  /// // encircle the common edge? That is, is E(d1,1) the common vertex?
+  /// const bool ccw = E(d1,1) == E(d2,0) || E(d1,1) == E(d2,1);
+  /// std::vector<int> Nf;
+  /// {
+  ///   std::vector<int> Nv;
+  ///   igl::circulation(d1,ccw,F,EMAP,EF,EI,Nv,Nf);
+  /// }
+  /// ```
   IGL_INLINE bool collapse_edge(
     const int e,
     const Eigen::RowVectorXd & p,

include/igl/collapse_edge_would_create_intersections.cpp

@@ -0,0 +1,163 @@
+#include "collapse_edge_would_create_intersections.h"
+#include "AABB.h"
+#include "circulation.h"
+#include "writePLY.h"
+#include "triangle_triangle_intersect.h"
+#include <Eigen/Geometry>
+#include <vector>
+#include <iostream>
+#include <algorithm>
+#include <cassert>
+
+IGL_INLINE bool igl::collapse_edge_would_create_intersections(
+  const int e,
+  const Eigen::RowVectorXd & p,
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXi & E,
+  const Eigen::VectorXi & EMAP,
+  const Eigen::MatrixXi & EF,
+  const Eigen::MatrixXi & EI,
+  const igl::AABB<Eigen::MatrixXd,3> & tree,
+  const int inf_face_id)
+{
+  // Merge two lists of integers
+  const auto merge = [&](
+    const std::vector<int> & A, const std::vector<int> & B)->
+    std::vector<int>
+  {
+    std::vector<int> C;
+    C.reserve( A.size() + B.size() ); // preallocate memory
+    C.insert( C.end(), A.begin(), A.end() );
+    C.insert( C.end(), B.begin(), B.end() );
+    // https://stackoverflow.com/a/1041939/148668
+    std::sort( C.begin(), C.end() );
+    C.erase( std::unique( C.begin(), C.end() ), C.end() );
+    return C;
+  };
+
+  std::vector<int> old_one_ring;
+  {
+    std::vector<int> Nsv,Nsf,Ndv,Ndf;
+    igl::circulation(e, true,F,EMAP,EF,EI,Nsv,Nsf);
+    igl::circulation(e,false,F,EMAP,EF,EI,Ndv,Ndf);
+    old_one_ring = merge(Nsf,Ndf);
+  }
+  int f1 = EF(e,0);
+  int f2 = EF(e,1);
+  std::vector<int> new_one_ring = old_one_ring;
+  // erase if ==f1 or ==f2
+  new_one_ring.erase(
+    std::remove(new_one_ring.begin(), new_one_ring.end(), f1), 
+    new_one_ring.end());
+  new_one_ring.erase(
+    std::remove(new_one_ring.begin(), new_one_ring.end(), f2), 
+    new_one_ring.end());
+
+
+  // big box containing new_one_ring
+  Eigen::AlignedBox<double,3> big_box;
+  // Extend box by placement point
+  big_box.extend(p.transpose());
+  // Extend box by all other corners (skipping old edge vertices)
+  for(const auto f : new_one_ring)
+  {
+    Eigen::RowVector3d corners[3];
+    for(int c = 0;c<3;c++)
+    {
+      if(F(f,c) == E(e,0) || F(f,c) == E(e,1))
+      {
+        corners[c] = p;
+      }else
+      {
+        corners[c] = V.row(F(f,c));
+        big_box.extend(V.row(F(f,c)).transpose());
+      }
+    }
+    // Degenerate triangles are considered intersections
+    if((corners[0]-corners[1]).cross(corners[0]-corners[2]).squaredNorm() < 1e-16)
+    {
+      return true;
+    }
+  }
+  
+
+  std::vector<const igl::AABB<Eigen::MatrixXd,3>*> candidates;
+  tree.append_intersecting_leaves(big_box,candidates);
+  
+
+  // Exclude any candidates that are in old_one_ring.
+  // consider using unordered_set above so that this is O(n+m) rather than O(nm)
+  candidates.erase(
+    std::remove_if(candidates.begin(), candidates.end(),
+        [&](const igl::AABB<Eigen::MatrixXd,3>* candidate) {
+            return std::find(old_one_ring.begin(), old_one_ring.end(), candidate->m_primitive) != old_one_ring.end();
+        }),
+    candidates.end());
+  // print candidates
+  //const bool stinker = e==2581;
+  constexpr bool stinker = false;
+  if(stinker)
+  {
+    igl::writePLY("before.ply",V,F);
+    std::cout<<"Ee = ["<<E(e,0)<<" "<<E(e,1)<<"]+1;"<<std::endl;
+    std::cout<<"p = ["<<p<<"];"<<std::endl;
+    // print new_one_ring as matlab vector of indices
+    std::cout<<"new_one_ring = [";
+    for(const auto f : new_one_ring)
+    {
+      std::cout<<f<<" ";
+    }
+    std::cout<<"]+1;"<<std::endl;
+    // print candidates as matlab vector of indices
+    std::cout<<"candidates = [";
+    for(const auto * candidate : candidates)
+    {
+      std::cout<<candidate->m_primitive<<" ";
+    }
+    std::cout<<"]+1;"<<std::endl;
+  }
+  
+  // For each pair of candidate and new_one_ring, check if they intersect
+  bool found_intersection = false;
+  for(const int & f : new_one_ring)
+  {
+    if(inf_face_id >= 0 && f >= inf_face_id) { continue; }
+
+    Eigen::AlignedBox<double,3> small_box;
+    small_box.extend(p.transpose());
+    for(int c = 0;c<3;c++)
+    {
+      if(F(f,c) != E(e,0) && F(f,c) != E(e,1))
+      {
+        small_box.extend(V.row(F(f,c)).transpose());
+      }
+    }
+    for(const auto * candidate : candidates)
+    {
+      const int g = candidate->m_primitive;
+      //constexpr bool inner_stinker = false;
+      const bool inner_stinker = stinker && (f==1492 && g==1554);
+      if(inner_stinker){ printf("  f: %d g: %d\n",f,g); }
+      if(!small_box.intersects(candidate->m_box))
+      {
+        if(inner_stinker){ printf("  ✅ boxes don't overlap\n"); }
+        continue;
+      }
+      // Corner replaced by p
+      int c;
+      for(c = 0;c<3;c++)
+      {
+        if(F(f,c) == E(e,0) || F(f,c) == E(e,1))
+        {
+          break;
+        }
+      }
+      assert(c<3);
+      found_intersection = triangle_triangle_intersect(V,F,E,EMAP,EF,f,c,p,g);
+      if(found_intersection) { break; }
+    }
+    if(found_intersection) { break; }
+  }
+  return found_intersection;
+}

include/igl/collapse_edge_would_create_intersections.h

@@ -0,0 +1,44 @@
+#ifndef IGL_COLLAPSE_EDGE_WOULD_CREATE_INTERSECTIONS_H
+#define IGL_COLLAPSE_EDGE_WOULD_CREATE_INTERSECTIONS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  /// Determine if collapse the edge `e` would create new intersections.
+  ///
+  /// @param[in] e  index into E of edge to try to collapse. E(e,:) = [s d] or [d s] so
+  ///     that s<d, then d is collapsed to s.
+  /// @param[in] p  dim list of vertex position where to place merged vertex
+  /// [mesh inputs]
+  /// @param[in,out] V  #V by dim list of vertex positions, lesser index of E(e,:) will be set
+  ///     to midpoint of edge.
+  /// @param[in,out] F  #F by 3 list of face indices into V.
+  /// @param[in,out] E  #E by 2 list of edge indices into V.
+  /// @param[in,out] EMAP #F*3 list of indices into E, mapping each directed edge to unique
+  ///     unique edge in E
+  /// @param[in,out] EF  #E by 2 list of edge flaps, EF(e,0)=f means e=(i-->j) is the edge of
+  ///     F(f,:) opposite the vth corner, where EI(e,0)=v. Similarly EF(e,1) "
+  ///     e=(j->i)
+  /// @param[in,out] EI  #E by 2 list of edge flap corners (see above).
+  /// [mesh inputs]
+  /// @param[in] tree AABB tree whose leaves correspond to the current
+  ///   (non-null) faces in (V,F)
+  ///
+  /// \see collapse_edge
+  template <typename DerivedV, int DIM> class AABB;
+  IGL_INLINE bool collapse_edge_would_create_intersections(
+    const int e,
+    const Eigen::RowVectorXd & p,
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXi & E,
+    const Eigen::VectorXi & EMAP,
+    const Eigen::MatrixXi & EF,
+    const Eigen::MatrixXi & EI,
+    const igl::AABB<Eigen::MatrixXd,3> & tree,
+    const int inf_face_id = -1);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "collapse_edge_would_create_intersections.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/is_self_intersecting.cpp

@@ -0,0 +1,63 @@
+#include "is_self_intersecting.h"
+#include "../../find.h"
+#include "../../doublearea.h"
+#include "../../remove_unreferenced.h"
+#include "RemeshSelfIntersectionsParam.h"
+#include "remesh_self_intersections.h"
+#include "../../collapse_edge.h"
+
+template <
+  typename DerivedV,
+  typename DerivedF>
+bool igl::copyleft::cgal::is_self_intersecting(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F)
+{
+  assert(V.cols() == 3);
+  assert(F.cols() == 3);
+  using MatrixV = Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, 3>;
+  using MatrixF = Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, 3>;
+  const auto valid = 
+    igl::find((F.array() != IGL_COLLAPSE_EDGE_NULL).rowwise().any().eval());
+  // Extract only the valid faces
+  MatrixF FF = F(valid, Eigen::all);
+  // Remove unreferneced vertices
+  MatrixV VV;
+  {
+    Eigen::VectorXi I;
+    igl::remove_unreferenced(V,MatrixF(FF),VV,FF,I);
+  }
+  Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,1> A;
+  igl::doublearea(VV,FF,A);
+  if(A.minCoeff() <= 0)
+  {
+    return true;
+  }
+  if(
+       (FF.array().col(0) == FF.array().col(1)).any() ||
+       (FF.array().col(1) == FF.array().col(2)).any() ||
+       (FF.array().col(2) == FF.array().col(0)).any())
+
+  {
+    return true;
+  }
+
+  // check for self-intersections VV,FF
+  igl::copyleft::cgal::RemeshSelfIntersectionsParam params;
+  params.detect_only = true;
+  params.first_only = true;
+  Eigen::MatrixXi IF;
+  Eigen::VectorXi J,IM;
+  {
+    MatrixV tempV;
+    MatrixF tempF;
+    igl::copyleft::cgal::remesh_self_intersections(
+      V,F,params,tempV,tempF,IF,J,IM);
+  }
+  return IF.rows() > 0;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::copyleft::cgal::is_self_intersecting<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&);
+#endif

include/igl/copyleft/cgal/is_self_intersecting.h

@@ -0,0 +1,36 @@
+#ifndef IGL_COPYLEFT_CGAL_IS_SELF_INTERSECTING_H
+#define IGL_COPYLEFT_CGAL_IS_SELF_INTERSECTING_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      /// Determine if a mesh has _any_ self-intersections. Skips any
+      /// `(IGL_COLLAPSE_EDGE_NULL,IGL_COLLAPSE_EDGE_NULL,IGL_COLLAPSE_EDGE_NULL)`
+      /// faces and returns true if any faces have zero area.
+      /// 
+      /// @param[in] V  #V by 3 list of vertex positions
+      /// @param[in] F  #F by 3 list of triangle indices into V
+      /// @return true if any faces intersect
+      ///
+      /// \see remesh_self_intersections, SelfIntersectMesh
+      template <
+        typename DerivedV,
+        typename DerivedF>
+      bool is_self_intersecting(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "is_self_intersecting.cpp"
+#endif
+
+#endif

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -508,28 +508,19 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick>>> const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>>>, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>>>>>> const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1>>&, 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>>&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3>> const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick>>> const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>>>, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>>>>>> const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>>&, 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>>&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3>> const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick>>> const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>>>, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>>>>>> const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3>>&, 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>>&);
-//MARK
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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> >&);
-// generated by autoexplicit.sh
-template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3> >&, 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> >&);
-// generated by autoexplicit.sh
 template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> >&, 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> >&);
-// Alec
+template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> >&, 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> >&);
-// generated by autoexplicit.sh
 template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, -1, 1, -1, -1> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3>, 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::MatrixBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 1, -1, 3> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epeck, 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, CGAL::Epick, 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, -1, 0, -1, -1>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+template void igl::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epeck, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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::copyleft::cgal::remesh_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, CGAL::Epick, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, std::map<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > >, std::less<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index const, std::vector<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object>, std::allocator<std::pair<Eigen::Matrix<int, -1, 3, 0, -1, 3>::Index, CGAL::Object> > > > > > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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> >&);
 //
 // To-do: we should really stop using Eigen::Index or anything but `int` for
 // index types, it causes this bloating of templates below and chaos on windows.

include/igl/copyleft/cgal/remesh_self_intersections.cpp

@@ -134,6 +134,8 @@ IGL_INLINE void igl::copyleft::cgal::remesh_self_intersections(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::copyleft::cgal::remesh_self_intersections<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::Matrix<int, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3>, 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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -1, 3, 0, -1, 3> >&, 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> >&);
 // generated by autoexplicit.sh
 template void igl::copyleft::cgal::remesh_self_intersections<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<CGAL::Epeck::FT, -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::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, igl::copyleft::cgal::RemeshSelfIntersectionsParam const&, Eigen::PlainObjectBase<Eigen::Matrix<CGAL::Epeck::FT, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/decimate.cpp

@@ -9,6 +9,8 @@
 #include "collapse_edge.h"
 #include "edge_flaps.h"
 #include "decimate_trivial_callbacks.h"
+#include "AABB.h"
+#include "intersection_blocking_collapse_edge_callbacks.h"
 #include "is_edge_manifold.h"
 #include "remove_unreferenced.h"
 #include "find.h"
@@ -20,12 +22,19 @@
 IGL_INLINE bool igl::decimate(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
-  const size_t max_m,
+  const int max_m,
+  const bool block_intersections,
   Eigen::MatrixXd & U,
   Eigen::MatrixXi & G,
   Eigen::VectorXi & J,
   Eigen::VectorXi & I)
 {
+  igl::AABB<Eigen::MatrixXd, 3> * tree = nullptr;
+  if(block_intersections)
+  {
+    tree = new igl::AABB<Eigen::MatrixXd, 3>();
+    tree->init(V,F);
+  }
   // Original number of faces
   const int orig_m = F.rows();
   // Tracking number of faces
@@ -49,16 +58,23 @@ IGL_INLINE bool igl::decimate(
       return false;
     }
   }
-  decimate_pre_collapse_callback always_try;
-  decimate_post_collapse_callback never_care;
-  decimate_trivial_callbacks(always_try,never_care);
+  decimate_pre_collapse_callback pre_collapse;
+  decimate_post_collapse_callback post_collapse;
+  decimate_trivial_callbacks(pre_collapse,post_collapse);
+  if(block_intersections)
+  {
+    igl::intersection_blocking_collapse_edge_callbacks(
+      pre_collapse, post_collapse, // These will get copied as needed
+      tree,
+      pre_collapse, post_collapse);
+  }
   bool ret = decimate(
     VO,
     FO,
     shortest_edge_and_midpoint,
     max_faces_stopping_condition(m,orig_m,max_m),
-    always_try,
-    never_care,
+    pre_collapse,
+    post_collapse,
     E,
     EMAP,
     EF,
@@ -73,21 +89,11 @@ IGL_INLINE bool igl::decimate(
   Eigen::VectorXi _1,I2;
   igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1,I2);
   I = I(I2).eval();
+  assert(tree == nullptr || tree == tree->root());
+  delete tree;
   return ret;
 }
 
-IGL_INLINE bool igl::decimate(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  const size_t max_m,
-  Eigen::MatrixXd & U,
-  Eigen::MatrixXi & G,
-  Eigen::VectorXi & J)
-{
-  Eigen::VectorXi I;
-  return igl::decimate(V,F,max_m,U,G,J,I);
-}
-
 IGL_INLINE bool igl::decimate(
   const Eigen::MatrixXd & OV,
   const Eigen::MatrixXi & OF,

include/igl/decimate.h

@@ -82,6 +82,8 @@ namespace igl
   /// @param[in] V  #V by dim list of vertex positions
   /// @param[in] F  #F by 3 list of face indices into V.
   /// @param[in] max_m  desired number of output faces
+  /// @param[in] block_intersections  whether to block intersections (see
+  ///   intersection_blocking_collapse_edge_callbacks)
   /// @param[out] U  #U by dim list of output vertex posistions (can be same ref as V)
   /// @param[out] G  #G by 3 list of output face indices into U (can be same ref as G)
   /// @param[out] J  #G list of indices into F of birth face
@@ -90,19 +92,12 @@ namespace igl
   IGL_INLINE bool decimate(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
-    const size_t max_m,
+    const int max_m,
+    const bool block_intersections,
     Eigen::MatrixXd & U,
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J,
     Eigen::VectorXi & I);
-  /// \overload
-  IGL_INLINE bool decimate(
-    const Eigen::MatrixXd & V,
-    const Eigen::MatrixXi & F,
-    const size_t max_m,
-    Eigen::MatrixXd & U,
-    Eigen::MatrixXi & G,
-    Eigen::VectorXi & J);
 
   /// Collapses edges of a **closed manifold mesh** (V,F) using user defined
   /// callbacks in a priority queue. Functions control the cost and placement of each collapse the

include/igl/fast_find_intersections.cpp

@@ -8,118 +8,246 @@
 #include "fast_find_intersections.h"
 #include "AABB.h"
 #include "tri_tri_intersect.h"
-
-
+#include "triangle_triangle_intersect_shared_edge.h"
+#include "triangle_triangle_intersect_shared_vertex.h"
+#include <stdio.h>
 
 template <
   typename DerivedV1,
   typename DerivedF1,
   typename DerivedV2,
   typename DerivedF2,
-  typename DerivedI,
-  typename DerivedE>
-IGL_INLINE void igl::fast_find_intersections(
-  const Eigen::MatrixBase<DerivedV1>& V1,
-  const Eigen::MatrixBase<DerivedF1>& F1,
-  const Eigen::MatrixBase<DerivedV2>& V2,
-  const Eigen::MatrixBase<DerivedF2>& F2,
-        Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
-        Eigen::PlainObjectBase<DerivedE>& edges )
+  typename DerivedIF,
+  typename DerivedEV,
+  typename DerivedEE,
+  typename DerivedEI>
+IGL_INLINE bool igl::fast_find_intersections(
+  const igl::AABB<DerivedV1,3> & tree,
+  const Eigen::MatrixBase<DerivedV1> & V1,
+  const Eigen::MatrixBase<DerivedF1> & F1,
+  const Eigen::MatrixBase<DerivedV2> & V2,
+  const Eigen::MatrixBase<DerivedF2> & F2,
+  const bool detect_only,
+  const bool first_only,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedEV> & EV,
+  Eigen::PlainObjectBase<DerivedEE> & EE,
+  Eigen::PlainObjectBase<DerivedEI> & EI)
 {
+  constexpr bool stinker = false;
   using AABBTree=igl::AABB<DerivedV1,3>;
-  AABBTree tree;
-  tree.init(V1,F1);
-
-  fast_find_intersections(tree, V1, F1, V2, F2, intersect_pairs, edges);
-}
-
-
-template <
-  typename DerivedV1,
-  typename DerivedF1,
-  typename DerivedV2,
-  typename DerivedF2,
-  typename DerivedI,
-  typename DerivedE>
-IGL_INLINE void igl::fast_find_intersections(
-  const AABB<DerivedV1,3>           & tree,
-  const Eigen::MatrixBase<DerivedV1>& V1,
-  const Eigen::MatrixBase<DerivedF1>& F1,
-  const Eigen::MatrixBase<DerivedV2>& V2,
-  const Eigen::MatrixBase<DerivedF2>& F2,
-        Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
-        Eigen::PlainObjectBase<DerivedE>& edges )
-{
-    using AABBTree=igl::AABB<DerivedV1,3>;
-    using Scalar=typename DerivedV1::Scalar;
-    using BBOX=Eigen::AlignedBox<Scalar,3>;
-    using VERTEX=Eigen::Matrix<typename DerivedE::Scalar,1,3,Eigen::RowMajor>;
-
-    std::vector<VERTEX> _edges;
-    std::vector<int>    _intersect_pairs;
+  using Scalar=typename DerivedV1::Scalar;
+  static_assert(
+    std::is_same<Scalar,typename DerivedV2::Scalar>::value,
+    "V1 and V2 must have the same scalar type");
+  static_assert(
+    std::is_same<Scalar,typename DerivedEV::Scalar>::value,
+    "V1 and EV must have the same scalar type");
+  using RowVector3S=Eigen::Matrix<Scalar,1,3>;
+  if(first_only){ assert(detect_only && "first_only must imply detect_only"); }
+
+  // Determine if V1,F1 and V2,F2 point to the same data
+  const bool self_test = (&V1 == &V2) && (&F1 == &F2);
+  if(stinker){ printf("%s\n",self_test?"🍎&(V1,F1) == 🍎&(V2,F2)":"🍎≠🍊"); }
+
+  int num_if = 0;
+  int num_ee = 0;
+  const auto append_intersection = 
+    [&IF,&EV,&EE,&EI,&detect_only,&num_if,&num_ee](
+      const int f1, 
+      const int f2,
+      const bool coplanar,
+      const RowVector3S & v1,
+      const RowVector3S & v2)
+  {
+    if(!coplanar && !detect_only)
+    {
+      if(num_ee >= EE.rows())
+      {
+        EE.conservativeResize(2*EE.rows()+1,2);
+        EI.conservativeResize(EE.rows());
+        EV.conservativeResize(2*EE.rows(),3);
+      }
+      EI(num_ee) = num_if;
+      EV.row(2*num_ee+0) = v1;
+      EV.row(2*num_ee+1) = v2;
+      EE.row(num_ee) << 2*num_ee+0, 2*num_ee+1;
+      num_ee++;
+    }
 
-    for(int i=0; i<F2.rows(); ++i)
+    if(num_if >= IF.rows())
     {
-      BBOX tri_box;
+      IF.conservativeResize(2*IF.rows()+1,2);
+    }
+    IF.row(num_if++) << f1,f2;
+  };
 
-      for(int j=0;j<3;++j)
-        tri_box.extend( V2.row( F2(i,j) ).transpose() );
-      
-      // find leaf nodes containing intersecting tri_box
-      // need to specify exact type instead of auto to allow recursion
-      std::function<void(const AABBTree &,int)> check_intersect = 
-        [&](const AABBTree &t,int d) -> void
+  // Returns corner in ith face opposite of shared edge; -1 otherwise
+  const auto shared_edge = [&F1](const int f, const int g)->int
+  {
+    for(int c = 0;c<3;c++)
+    {
+      int s = F1(f,(c+1)%3);
+      int d = F1(f,(c+2)%3);
+      for(int e = 0;e<3;e++)
       {
-        if(t.m_primitive != -1) //check for the actual intersection //t.is_leaf()
+        // Find in opposite direction on jth face
+        if(F1(g,e) == d && F1(g,(e+1)%3) == s)
         {
-          bool coplanar=false;
-          VERTEX edge1,edge2;
-
-          if(igl::tri_tri_intersection_test_3d(
-            V2.row(F2(i,0)),            V2.row(F2(i,1)),            V2.row(F2(i,2)),
-            V1.row(F1(t.m_primitive,0)),V1.row(F1(t.m_primitive,1)),V1.row(F1(t.m_primitive,2)),
-            coplanar,
-            edge1,edge2))
-          {
-            if(!coplanar)
-            {
-              _intersect_pairs.push_back(t.m_primitive);
-              _intersect_pairs.push_back(i);
-
-              _edges.push_back(edge1);
-              _edges.push_back(edge2);
-            }
-          }
-        } else {
-          if(t.m_box.intersects( tri_box )) 
-          { // need to check all branches 
-            check_intersect(*t.m_left, d+1);
-            check_intersect(*t.m_right,d+1);
-          }
+          return c;
         }
-      };
-
-      // run actual search
-      check_intersect(tree, 0);
+      }
     }
-    edges.resize(_edges.size(), 3);
-
-    for(int i=0; i!=_edges.size(); ++i)
+    return -1;
+  };
+  // Returns corner of shared vertex in ith face; -1 otherwise
+  const auto shared_vertex = [&F1](const int f, const int g, int & sf, int & sg)->bool
+  {
+    for(sf = 0;sf<3;sf++)
     {
-      edges.row(i) = _edges[i];
+      for(sg = 0;sg<3;sg++)
+      {
+        if(F1(g,sg) == F1(f,sf))
+        {
+          return true;
+        }
+      }
     }
-
-    intersect_pairs.resize(_intersect_pairs.size()/2,2);
-    for(int i=0; i!=_intersect_pairs.size(); ++i)
+    return false;
+  };
+
+  RowVector3S dummy;
+  for(int f2 = 0; f2<F2.rows(); ++f2)
+  {
+    if(stinker){ printf("f2: %d\n",f2); }
+    Eigen::AlignedBox<Scalar,3> box;
+    box.extend( V2.row( F2(f2,0) ).transpose() );
+    box.extend( V2.row( F2(f2,1) ).transpose() );
+    box.extend( V2.row( F2(f2,2) ).transpose() );
+    std::vector<const AABBTree*> candidates;
+    tree.append_intersecting_leaves(box, candidates);
+    for(const auto * candidate : candidates)
     {
-      intersect_pairs(i/2, i%2) = _intersect_pairs[i];
+      const int f1 = candidate->m_primitive;
+      if(stinker){ printf("  f1: %d\n",f1); }
+      bool found_intersection = false;
+      bool yes_shared_verted = false;
+      bool yes_shared_edge = false;
+      if(self_test)
+      {
+        // Skip self-test and direction f2>=f1 (assumes by symmetry we'll find
+        // the other direction since we're testing all pairs)
+        if(f1 >= f2)
+        {
+          if(stinker){ printf("    ⏭\n"); }
+          continue;
+        }
+        const int c = shared_edge(f1,f2);
+        yes_shared_edge = c != -1;
+        if(yes_shared_edge)
+        {
+          if(stinker){ printf("    ⚠️  shared edge\n"); }
+          if(stinker)
+          {
+            printf("    %d: %d %d %d\n",f1,F1(f1,0),F1(f1,1),F1(f1,2));
+            printf("    %d: %d %d %d\n",f2,F1(f2,0),F1(f2,1),F1(f2,2));
+            printf("   edge: %d %d\n",F1(f1,(c+1)%3),F1(f1,(c+2)%3));
+          }
+          found_intersection = igl::triangle_triangle_intersect_shared_edge(
+            V1,F1,f1,c,V1.row(F1(f1,c)),f2,1e-8);
+          if(found_intersection)
+          {
+            append_intersection(f1,f2,true,dummy,dummy);
+          }
+        }else
+        {
+          int sf,sg;
+          yes_shared_verted = shared_vertex(f1,f2,sf,sg);
+          if(yes_shared_verted)
+          {
+            if(stinker){ printf("    ⚠️  shared vertex\n"); }
+            // Just to be sure. c≠sf
+            const int c = (sf+1)%3;
+            assert(F1(f1,sf) == F1(f2,sg));
+            found_intersection = igl::triangle_triangle_intersect_shared_vertex(
+              V1,F1,f1,sf,c,V1.row(F1(f1,c)),f2,sg,1e-14);
+            if(found_intersection && detect_only)
+            {
+              append_intersection(f1,f2,true,dummy,dummy);
+            }
+          }
+          
+        }
+      }
+      if(
+        !self_test || 
+        (!yes_shared_verted && !yes_shared_edge) || 
+        (yes_shared_verted && found_intersection && !detect_only))
+      {
+        bool coplanar = false;
+        RowVector3S v1,v2;
+        const bool tt_found_intersection = 
+          igl::tri_tri_intersection_test_3d(
+            V2.row(F2(f2,0)),V2.row(F2(f2,1)),V2.row(F2(f2,2)),
+            V1.row(F1(f1,0)),V1.row(F1(f1,1)),V1.row(F1(f1,2)),
+            coplanar,
+            v1,v2);
+        if(found_intersection && !tt_found_intersection)
+        {
+          // We failed to find the edge. Mark it as an intersection but don't
+          // include edge.
+          append_intersection(f1,f2,true,dummy,dummy);
+        }else if(tt_found_intersection)
+        {
+          found_intersection = true;
+          append_intersection(f1,f2,false,v1,v2);
+        }
+      }
+      if(stinker) { printf("    %s\n",found_intersection? "☠️":"❌"); }
+      if(num_if && first_only) { break; }
     }
+    if(num_if && first_only) { break; }
+  }
+  if(!detect_only)
+  {
+    EV.conservativeResize(2*num_ee,3);
+    EE.conservativeResize(num_ee,2);
+    EI.conservativeResize(num_ee,1);
+  }
+  IF.conservativeResize(num_if,2);
+  return IF.rows();
 }
 
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedIF,
+  typename DerivedEV,
+  typename DerivedEE,
+  typename DerivedEI>
+IGL_INLINE bool igl::fast_find_intersections(
+  const Eigen::MatrixBase<DerivedV1> & V1,
+  const Eigen::MatrixBase<DerivedF1> & F1,
+  const Eigen::MatrixBase<DerivedV2> & V2,
+  const Eigen::MatrixBase<DerivedF2> & F2,
+  const bool detect_only,
+  const bool first_only,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedEV> & EV,
+  Eigen::PlainObjectBase<DerivedEE> & EE,
+  Eigen::PlainObjectBase<DerivedEI> & EI)
+{
+  igl::AABB<DerivedV1,3> tree1;
+  tree1.init(V1,F1);
+  return fast_find_intersections(
+    tree1,V1,F1,V2,F2,detect_only,first_only,IF,EV,EE,EI);
+}
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template void igl::fast_find_intersections<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<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template void igl::fast_find_intersections<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<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template bool igl::fast_find_intersections<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 
 #endif

include/igl/fast_find_intersections.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2022 Vladimir S. FONOV <[email protected]>
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
 // 
 // 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 
@@ -10,20 +11,25 @@
 #define FAST_FIND_MESH_INTERSECT_H
 
 #include "igl_inline.h"
-#include "AABB.h"
 #include <Eigen/Core>
 
-namespace igl {
+namespace igl 
+{
+  template <typename DerivedV, int DIM> class AABB;
   /// Identify triangles where two meshes interesect 
   /// using AABBTree and tri_tri_intersection_test_3d.
   ///
-  /// @param[in] V1  #V by 3 list representing vertices on the first mesh
-  /// @param[in] F1  #F by 3 list representing triangles on the first mesh
-  /// @param[in] V2  #V by 3 list representing vertices on the second mesh
-  /// @param[in] F2  #F by 3 list representing triangles on the second mesh
-  /// @param[out] intersect_pairs  correspondance list of intersecting triangles
-  ///                    column 0 - mesh 1, column 1 - mesh2  
-  /// @param[out] edges      list of pairs of intersection edges
+  /// @param[in] V1  #V1 by 3 list representing vertices on the first mesh
+  /// @param[in] F1  #F1 by 3 list representing triangles on the first mesh
+  /// @param[in] V2  #V2 by 3 list representing vertices on the second mesh
+  /// @param[in] F2  #F2 by 3 list representing triangles on the second mesh
+  /// @param[out] IF #IF by 2 list of intersecting triangle pairs, so that 
+  ///   F1(IF(i,0),:) intersects F2(IF(i,1),:)
+  /// @param[out] EV #EV by 3 list of vertices definining intersection segments
+  /// for non-coplanar intersections
+  /// @param[out] EE #EE by 2 list of edges indices into rows of EV
+  /// @param[out] EI #EI by 1 list of indices into rows IF indicating source of
+  ///   intersection.
   ///
   /// \see copyleft::cgal::intersect_other
   template <
@@ -31,32 +37,44 @@ namespace igl {
     typename DerivedF1,
     typename DerivedV2,
     typename DerivedF2,
-    typename DerivedI,
-    typename DerivedE>
-  IGL_INLINE void fast_find_intersections(
-    const Eigen::MatrixBase<DerivedV1>& V1,
-    const Eigen::MatrixBase<DerivedF1>& F1,
-    const Eigen::MatrixBase<DerivedV2>& V2,
-    const Eigen::MatrixBase<DerivedF2>& F2,
-          Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
-          Eigen::PlainObjectBase<DerivedE>& edges );
+    typename DerivedIF,
+    typename DerivedEV,
+    typename DerivedEE,
+    typename DerivedEI>
+  IGL_INLINE bool fast_find_intersections(
+    const AABB<DerivedV1,3> & tree,
+    const Eigen::MatrixBase<DerivedV1> & V1,
+    const Eigen::MatrixBase<DerivedF1> & F1,
+    const Eigen::MatrixBase<DerivedV2> & V2,
+    const Eigen::MatrixBase<DerivedF2> & F2,
+    const bool detect_only,
+    const bool first_only,
+    Eigen::PlainObjectBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedEV> & EV,
+    Eigen::PlainObjectBase<DerivedEE> & EE,
+    Eigen::PlainObjectBase<DerivedEI> & EI);
   /// \overload
-  /// @param[in] tree - AABB tree bult from the first mesh
+  /// \brief Tree built internally.
   template <
     typename DerivedV1,
     typename DerivedF1,
     typename DerivedV2,
     typename DerivedF2,
-    typename DerivedI,
-    typename DerivedE>
-  IGL_INLINE void fast_find_intersections(
-    const AABB<DerivedV1,3>           & tree,
-    const Eigen::MatrixBase<DerivedV1>& V1,
-    const Eigen::MatrixBase<DerivedF1>& F1,
-    const Eigen::MatrixBase<DerivedV2>& V2,
-    const Eigen::MatrixBase<DerivedF2>& F2,
-          Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
-          Eigen::PlainObjectBase<DerivedE>& edges );
+    typename DerivedIF,
+    typename DerivedEV,
+    typename DerivedEE,
+    typename DerivedEI>
+  IGL_INLINE bool fast_find_intersections(
+    const Eigen::MatrixBase<DerivedV1> & V1,
+    const Eigen::MatrixBase<DerivedF1> & F1,
+    const Eigen::MatrixBase<DerivedV2> & V2,
+    const Eigen::MatrixBase<DerivedF2> & F2,
+    const bool detect_only,
+    const bool first_only,
+    Eigen::PlainObjectBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedEV> & EV,
+    Eigen::PlainObjectBase<DerivedEE> & EE,
+    Eigen::PlainObjectBase<DerivedEI> & EI);
 };
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/fast_find_self_intersections.cpp

@@ -1,201 +1,39 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2022 Vladimir S. FONOV <[email protected]>
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
 // 
 // 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 "fast_find_self_intersections.h"
-#include "AABB.h"
-#include "tri_tri_intersect.h"
-
-namespace igl{
-namespace internal {
-
-    // helper function, to check if two faces have shared vertices
-    template <typename Derived>
-    bool adjacent_faces(const Eigen::MatrixBase<Derived>& A,
-                        const Eigen::MatrixBase<Derived>& B)
-    {
-        for(int i=0;i<3;++i)
-          for(int j=0;j<3;j++)
-          {
-            if(A(i)==B(j))
-              return true;
-          }
-        return false;
-    }
-
-}
-}
+#include "fast_find_intersections.h"
 
 template <
   typename DerivedV,
   typename DerivedF,
-  typename DerivedI>
+  typename DerivedIF,
+  typename DerivedEV,
+  typename DerivedEE,
+  typename DerivedEI>
 IGL_INLINE bool igl::fast_find_self_intersections(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedF>& F,
-        Eigen::PlainObjectBase<DerivedI>& intersect)
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const bool detect_only,
+  const bool first_only,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedEV> & EV,
+  Eigen::PlainObjectBase<DerivedEE> & EE,
+  Eigen::PlainObjectBase<DerivedEI> & EI)
 {
-    using Scalar=typename DerivedV::Scalar;
-    using BBOX=Eigen::AlignedBox<Scalar,3>;
-    using AABBTree=igl::AABB<DerivedV,3>;
-    AABBTree tree;
-
-    tree.init(V,F);
-    bool _intersects=false;
-
-    intersect.resize(F.rows(),1);
-    intersect.setConstant(0);
-
-    for(int i=0; i<F.rows(); ++i)
-    {
-      if( intersect(i) ) continue;
-
-      BBOX tri_box;
-
-      for(int j=0;j<3;++j)
-        tri_box.extend( V.row( F(i,j) ).transpose() );
-      
-      // find leaf nodes containing intersecting tri_box
-      // need to declare full type to enable recursion
-      std::function<bool(const AABBTree &,int)> check_intersect = 
-        [&](const AABBTree &t,int d) -> bool
-      {
-        if(t.m_primitive != -1) //check for the actual intersection (is_leaf)
-        {
-          if(t.m_primitive==i) //itself
-              return false;
-          if(igl::internal::adjacent_faces(F.row(i), F.row(t.m_primitive)) )
-              return false;
-
-          bool coplanar=false;
-          Eigen::Matrix<Scalar,1,3,Eigen::RowMajor> edge1,edge2;
-
-          if(igl::tri_tri_intersection_test_3d(
-            V.row(F(i,0)),V.row(F(i,1)),V.row(F(i,2)),
-            V.row(F(t.m_primitive,0)),V.row(F(t.m_primitive,1)),V.row(F(t.m_primitive,2)),
-            coplanar,
-            edge1,edge2))
-          {
-            if(!coplanar)
-            {
-              intersect(i)=1;
-              intersect(t.m_primitive)=1;
-              return true;
-            }
-          }
-        } else {
-          if(t.m_box.intersects(tri_box)) {
-            // need to check both subtrees
-            bool r1=check_intersect(*t.m_left ,d+1);
-            bool r2=check_intersect(*t.m_right,d+1);
-            return r1||r2;
-          }
-        }
-        return false;
-      };
-
-      bool r=check_intersect(tree,0);
-      _intersects = _intersects || r;
-    }
-    return _intersects;
+  // This is really just a wrapper around fast_find_intersections which will
+  // internally detect that V,F are the second set
+  return igl::fast_find_intersections(
+    V,F,V,F,detect_only,first_only,IF,EV,EE,EI);
 }
 
-
-template <
-  typename DerivedV,
-  typename DerivedF,
-  typename DerivedI,
-  typename DerivedE>
-IGL_INLINE bool igl::fast_find_self_intersections(
-  const Eigen::MatrixBase<DerivedV>& V,
-  const Eigen::MatrixBase<DerivedF>& F,
-        Eigen::PlainObjectBase<DerivedI>& intersect,
-        Eigen::PlainObjectBase<DerivedE>& edges )
-{
-    using Scalar=typename DerivedV::Scalar;
-    using BBOX=Eigen::AlignedBox<Scalar,3>;
-    using AABBTree=igl::AABB<DerivedV,3>;
-    using EDGE=Eigen::Matrix<typename DerivedE::Scalar,1,3,Eigen::RowMajor>;
-
-    std::vector<EDGE> _edges;
-    AABBTree tree;
-
-    tree.init(V,F);
-    bool _intersects=false;
-
-    intersect.resize(F.rows(),1);
-    intersect.setConstant(0);
-
-    for(int i=0; i<F.rows(); ++i)
-    {
-      if( intersect(i) ) continue;
-
-      BBOX tri_box;
-
-      for(int j=0;j<3;++j)
-        tri_box.extend( V.row( F(i,j) ).transpose() );
-      
-      // find leaf nodes containing intersecting tri_box
-      std::function<bool(const AABBTree &,int)> check_intersect = 
-        [&](const AABBTree &t,int d) -> bool
-      {
-        if(t.m_primitive != -1) //check for the actual intersection //t.is_leaf()
-        {
-          if(t.m_primitive==i) //itself
-              return false;
-          if(igl::internal::adjacent_faces(F.row(i), F.row(t.m_primitive)) )
-              return false;
-
-          bool coplanar=false;
-          EDGE edge1,edge2;
-
-          if(igl::tri_tri_intersection_test_3d(
-            V.row(F(i,0)),            V.row(F(i,1)),            V.row(F(i,2)),
-            V.row(F(t.m_primitive,0)),V.row(F(t.m_primitive,1)),V.row(F(t.m_primitive,2)),
-            coplanar,
-            edge1,edge2))
-          {
-            if(!coplanar)
-            {
-              intersect(i)=1;
-              intersect(t.m_primitive)=1;
-              _edges.push_back(edge1);
-              _edges.push_back(edge2);
-              return true;
-            }
-          }
-        } else {
-          if(t.m_box.intersects( tri_box ))
-          {
-            bool r1=check_intersect(*t.m_left, d+1);
-            bool r2=check_intersect(*t.m_right,d+1);
-            return r1||r2;
-          }
-        }
-        return false;
-      };
-
-      // run actual search
-      bool r=check_intersect(tree, 0);
-      _intersects = _intersects || r;
-    }
-
-    edges.resize(_edges.size(),3);
-    int i=0;
-
-    for(auto e=std::begin(_edges);e!=std::end(_edges);++e,++i)
-    {
-      edges.row(i) = *e;
-    }
-    return _intersects;
-}
-
-
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-template bool igl::fast_find_self_intersections<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<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-template bool igl::fast_find_self_intersections<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
+template bool igl::fast_find_self_intersections<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<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, bool, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/fast_find_self_intersections.h

@@ -12,38 +12,38 @@
 #include "igl_inline.h"
 #include <Eigen/Core>
 
-namespace igl {
-
-  /// Identify triangles where mesh intersects itself
-  /// using AABBTree and tri_tri_intersection_test_3d
+namespace igl
+{
+  /// Identify triangles where two meshes interesect 
+  /// using AABBTree and tri_tri_intersection_test_3d.
   ///
-  /// @param[in] V  #V by 3 list representing vertices
-  /// @param[in] F  #F by 3 list representing triangles.
-  /// @param[out] intersect  #F by 1 indicator that triangle intersects anothe triangle
-  /// @param[out] edges      list of pairs of intersection edges
-  /// @return whether any self-interections were found
+  /// @param[in] V  #V by 3 list representing vertices on the first mesh
+  /// @param[in] F  #F by 3 list representing triangles on the first mesh
+  /// @param[out] IF #IF by 2 list of intersecting triangle pairs, so that 
+  ///   F1(IF(i,0),:) intersects F2(IF(i,1),:)
+  /// @param[out] EV #EV by 3 list of vertices definining intersection segments
+  /// for non-coplanar intersections
+  /// @param[out] EE #EE by 2 list of edges indices into rows of EV
+  /// @param[out] EI #EI by 1 list of indices into rows IF indicating source of
+  ///   intersection.
   ///
-  /// \see copyleft::cgal::remesh_self_intersections
-  template <
-    typename DerivedV,
-    typename DerivedF,
-    typename DerivedI,
-    typename DerivedE>
-  IGL_INLINE bool fast_find_self_intersections(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedF>& F,
-    Eigen::PlainObjectBase<DerivedI>& intersect,
-    Eigen::PlainObjectBase<DerivedE>& edges );
-  /// \overload
+  /// \see copyleft::cgal::SelfIntersectMesh
   template <
     typename DerivedV,
     typename DerivedF,
-    typename DerivedI>
+    typename DerivedIF,
+    typename DerivedEV,
+    typename DerivedEE,
+    typename DerivedEI>
   IGL_INLINE bool fast_find_self_intersections(
-    const Eigen::MatrixBase<DerivedV>& V,
-    const Eigen::MatrixBase<DerivedF>& F,
-    Eigen::PlainObjectBase<DerivedI>& intersect);
-
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const bool detect_only,
+    const bool first_only,
+    Eigen::PlainObjectBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedEV> & EV,
+    Eigen::PlainObjectBase<DerivedEE> & EE,
+    Eigen::PlainObjectBase<DerivedEI> & EI);
 };
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/intersection_blocking_collapse_edge_callbacks.cpp

@@ -0,0 +1,176 @@
+#include "intersection_blocking_collapse_edge_callbacks.h"
+#include "AABB.h"
+#include "circulation.h"
+#include "decimate_trivial_callbacks.h"
+#include "collapse_edge_would_create_intersections.h"
+
+// If debugging, it's a good idea to run this once in release mode to collect
+// the edge id then adjust below.
+//#define IGL_INTERSECTION_BLOCKING_COLLAPSE_EDGE_CALLBACKS_DEBUG
+#ifdef IGL_INTERSECTION_BLOCKING_COLLAPSE_EDGE_CALLBACKS_DEBUG
+#include "copyleft/cgal/is_self_intersecting.h"
+#include "writePLY.h"
+#endif
+
+void igl::intersection_blocking_collapse_edge_callbacks(
+  const igl::decimate_pre_collapse_callback  & orig_pre_collapse,
+  const igl::decimate_post_collapse_callback & orig_post_collapse,
+  const std::vector<igl::AABB<Eigen::MatrixXd,3> *> & leaves,
+        igl::AABB<Eigen::MatrixXd,3> * & tree,
+        igl::decimate_pre_collapse_callback  & pre_collapse,
+        igl::decimate_post_collapse_callback & post_collapse
+        )
+{
+  // Not clear padding would help much but could try it.
+  const double pad = 0;
+  const int leaves_size = (int)leaves.size();
+  // Capture the original callbacks by value so that caller can destory their
+  // copy.
+  pre_collapse = 
+    [orig_pre_collapse,leaves_size,&tree](
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      const Eigen::MatrixXi & E,
+      const Eigen::VectorXi & EMAP,
+      const Eigen::MatrixXi & EF,
+      const Eigen::MatrixXi & EI,
+      const igl::min_heap< std::tuple<double,int,int> > & Q,
+      const Eigen::VectorXi & EQ,
+      const Eigen::MatrixXd & C,
+      const int e)->bool
+    {
+      if(!orig_pre_collapse(V,F,E,EMAP,EF,EI,Q,EQ,C,e))
+      {
+        return false;
+      }
+      
+      // Check if there would be (new) intersections
+      return 
+        !igl::collapse_edge_would_create_intersections(
+          e,C.row(e).eval(),V,F,E,EMAP,EF,EI,*tree,leaves_size);
+    };
+  // leaves could also be captured by value
+  post_collapse =
+    [orig_post_collapse,leaves,pad,leaves_size,&tree](
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      const Eigen::MatrixXi & E,
+      const Eigen::VectorXi & EMAP,
+      const Eigen::MatrixXi & EF,
+      const Eigen::MatrixXi & EI,
+      const igl::min_heap< std::tuple<double,int,int> > & Q,
+      const Eigen::VectorXi & EQ,
+      const Eigen::MatrixXd & C,
+      const int e,
+      const int e1,
+      const int e2,
+      const int f1,
+      const int f2,
+      const bool collapsed)
+    {
+      if(collapsed)
+      {
+        // detach leaves of deleted faces
+        for(int f : {f1,f2})
+        {
+          // Skip faces that aren't in leaves (e.g., infinite faces)
+          if(f >= leaves_size) { continue; }
+          auto * sibling = leaves[f]->detach();
+          sibling->refit_lineage();
+          tree = sibling->root();
+          delete leaves[f];
+        }
+        // If finding `Nf` becomes a bottleneck we could remember it via
+        // `pre_collapse` the same way that
+        // `qslim_optimal_collapse_edge_callbacks` remembers `v1` and `v2`
+        const int m = F.rows();
+        const auto survivors = 
+          [&m,&e,&EF,&EI,&EMAP](const int f1, const int e1, int & d1)
+        {
+          int c;
+          for(c=0;c<3;c++)
+          {
+            d1 = EMAP(f1+c*m);
+            if((d1 != e) && (d1 != e1)) { break; }
+          }
+          assert(c<3);
+        };
+        int d1,d2;
+        survivors(f1,e1,d1);
+        survivors(f2,e2,d2);
+        // Will circulating by continuing in the CCW direction of E(d1,:)
+        // encircle the common edge? That is, is E(d1,1) the common vertex?
+        const bool ccw = E(d1,1) == E(d2,0) || E(d1,1) == E(d2,1);
+#ifndef NDEBUG
+        // Common vertex.
+        const int s = E(d1,ccw?1:0);
+        assert(s == E(d2,0) || s == E(d2,1));
+#endif
+        std::vector<int> Nf;
+        {
+          std::vector<int> Nv;
+          igl::circulation(d1,ccw,F,EMAP,EF,EI,Nv,Nf);
+        }
+        for(const int & f : Nf)
+        {
+          // Skip faces that aren't in leaves (e.g., infinite faces)
+          if(f >= leaves_size) { continue; }
+          Eigen::AlignedBox<double,3> box;
+          box
+            .extend(V.row(F(f,0)).transpose())
+            .extend(V.row(F(f,1)).transpose())
+            .extend(V.row(F(f,2)).transpose());
+          // Always grab root (returns self if no update)
+          tree = leaves[f]->update(box,pad)->root();
+          assert(tree == tree->root());
+        }
+          assert(tree == tree->root());
+      }
+#ifdef IGL_INTERSECTION_BLOCKING_COLLAPSE_EDGE_CALLBACKS_DEBUG
+#warning "🐌🐌🐌🐌🐌🐌🐌🐌🐌🐌🐌 Slow intersection checking..."
+      constexpr bool stinker = false;
+      //const bool stinker = e==2581;
+      if(stinker && igl::copyleft::cgal::is_self_intersecting(V,F))
+      {
+        igl::writePLY("after.ply",V,F);
+        printf("💩🛌@e=%d \n",e);
+        exit(1);
+      }
+#endif
+      // Finally. Run callback.
+      return orig_post_collapse(
+        V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2,collapsed);
+    };
+}
+
+IGL_INLINE void igl::intersection_blocking_collapse_edge_callbacks(
+  const igl::decimate_pre_collapse_callback  & orig_pre_collapse,
+  const igl::decimate_post_collapse_callback & orig_post_collapse,
+        igl::AABB<Eigen::MatrixXd,3> * & tree,
+        igl::decimate_pre_collapse_callback  & pre_collapse,
+        igl::decimate_post_collapse_callback & post_collapse)
+{
+  return intersection_blocking_collapse_edge_callbacks(
+    orig_pre_collapse,
+    orig_post_collapse,
+    tree->gather_leaves(),
+    tree,
+    pre_collapse,
+    post_collapse);
+}
+
+IGL_INLINE void igl::intersection_blocking_collapse_edge_callbacks(
+  igl::AABB<Eigen::MatrixXd,3> * & tree,
+  igl::decimate_pre_collapse_callback  & pre_collapse,
+  igl::decimate_post_collapse_callback & post_collapse)
+{
+  igl::decimate_pre_collapse_callback  always_try;
+  igl::decimate_post_collapse_callback never_care;
+  igl::decimate_trivial_callbacks(always_try,never_care);
+  intersection_blocking_collapse_edge_callbacks(
+    always_try,
+    never_care,
+    tree,
+    pre_collapse,
+    post_collapse);
+}

include/igl/intersection_blocking_collapse_edge_callbacks.h

@@ -0,0 +1,94 @@
+#ifndef IGL_INTERSECTION_BLOCKING_COLLAPSE_EDGE_CALLBACKS_H
+#define IGL_INTERSECTION_BLOCKING_COLLAPSE_EDGE_CALLBACKS_H
+#include "igl_inline.h"
+#include "decimate_callback_types.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Forward declaration
+  template <typename DerivedV, int DIM> class AABB;
+  /// Wrap around callbacks for decimate to prevent intersections from being
+  /// created.
+  ///
+  /// @param[in] orig_pre_collapse  original pre_collapse callback (callers copy
+  ///   can be destroyed)
+  /// @param[in] orig_post_collapse original post_collapse callback (ditto)
+  /// @param[in] leaves  list of pointers to leaves of tree (ditto)
+  /// @param[in,out] tree  pointer to AABB tree whose leaves will correspond to
+  ///   the current (non-null) faces in (V,F)  (callers copy can _not_ be
+  ///   destroyed – not the tree and not this pointer – until use of
+  ///   pre_/post_collapse is done)
+  /// @param[out] pre_collapse  new pre_collapse callback
+  /// @param[out] post_collapse  new post_collapse callback
+  ///
+  /// #### Example
+  ///
+  /// ```cpp
+  ///  // Build tree around mesh (must be edge-manifold, may have boundaries)
+  ///  igl::AABB<Eigen::MatrixXd, 3> * tree = new igl::AABB<Eigen::MatrixXd, 3>();
+  ///  tree->init(V,F);
+  ///  // Connect boundaries to dummy infinite vertex (after tree building)
+  ///  Eigen::MatrixXd VO;
+  ///  Eigen::MatrixXi FO;
+  ///  igl::connect_boundary_to_infinity(V,F,VO,FO);
+  ///  // prepare edge structures
+  ///  Eigen::VectorXi EMAP;
+  ///  Eigen::MatrixXi E,EF,EI;
+  ///  igl::edge_flaps(FO,E,EMAP,EF,EI);
+  ///  // prepare callbacks
+  ///  igl::decimate_cost_and_placement_callback cost_and_placement;
+  ///  igl::decimate_pre_collapse_callback  pre_collapse;
+  ///  igl::decimate_post_collapse_callback post_collapse;
+  ///  cost_and_placement = igl::shortest_edge_and_midpoint;
+  ///  igl::decimate_trivial_callbacks(pre_collapse,post_collapse);
+  ///  igl::intersection_blocking_collapse_edge_callbacks(
+  ///      pre_collapse, post_collapse, // These will get copied as needed
+  ///      tree,
+  ///      pre_collapse, post_collapse);
+  ///  int m = F.rows();
+  ///  const int orig_m = m;
+  ///  Eigen::MatrixXd U;
+  ///  Eigen::MatrixXi G;
+  ///  Eigen::VectorXi J,I;
+  ///  const bool ret = igl::decimate(
+  ///    VO, FO,
+  ///    cost_and_placement,
+  ///    igl::max_faces_stopping_condition(m,orig_m,target_m),
+  ///    pre_collapse,
+  ///    post_collapse,
+  ///    E, EMAP, EF, EI,
+  ///    U, G, J, I);
+  ///  G = G(igl::find((J.array()<orig_m).eval()), Eigen::all).eval();
+  ///  {
+  ///    Eigen::VectorXi _;
+  ///    igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_);
+  ///  }
+  /// ```
+  ///
+  /// \see decimate.h
+  IGL_INLINE void intersection_blocking_collapse_edge_callbacks(
+    const decimate_pre_collapse_callback  & orig_pre_collapse,
+    const decimate_post_collapse_callback & orig_post_collapse,
+    const std::vector<AABB<Eigen::MatrixXd,3> *> & leaves,
+          AABB<Eigen::MatrixXd,3> * & tree,
+          decimate_pre_collapse_callback  & pre_collapse,
+          decimate_post_collapse_callback & post_collapse);
+  /// \overload Same as above but computes leaves from tree
+  IGL_INLINE void intersection_blocking_collapse_edge_callbacks(
+    const decimate_pre_collapse_callback  & orig_pre_collapse,
+    const decimate_post_collapse_callback & orig_post_collapse,
+          AABB<Eigen::MatrixXd,3> * & tree,
+          decimate_pre_collapse_callback  & pre_collapse,
+          decimate_post_collapse_callback & post_collapse);
+  /// \overload Same as above but uses trivial callbacks
+  IGL_INLINE void intersection_blocking_collapse_edge_callbacks(
+    AABB<Eigen::MatrixXd,3> * & tree,
+    decimate_pre_collapse_callback  & pre_collapse,
+    decimate_post_collapse_callback & post_collapse);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intersection_blocking_collapse_edge_callbacks.cpp"
+#endif
+#endif

include/igl/pad_box.cpp

@@ -0,0 +1,17 @@
+#include "pad_box.h"
+
+template <typename Scalar, int DIM>
+IGL_INLINE void igl::pad_box(
+  const Scalar pad,
+  Eigen::AlignedBox<Scalar,DIM> & box)
+{
+  box.min().array() -= pad;
+  box.max().array() += pad;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::pad_box<double, 2>(double, Eigen::AlignedBox<double, 2>&);
+template void igl::pad_box<double, 3>(double, Eigen::AlignedBox<double, 3>&);
+#endif

include/igl/pad_box.h

@@ -0,0 +1,23 @@
+#ifndef IGL_PAD_BOX_H
+#define IGL_PAD_BOX_H
+
+#include <Eigen/Geometry>
+#include "igl_inline.h"
+
+namespace igl
+{
+  /// Pads a box by a given amount
+  ///
+  /// @param[in] pad: amount to pad each side of the box
+  /// @param[in,out] box  box to be padded.
+  template <typename Scalar, int DIM>
+  IGL_INLINE void pad_box(
+    const Scalar pad,
+    Eigen::AlignedBox<Scalar,DIM> & box);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "pad_box.cpp"
+#endif
+
+#endif

include/igl/per_vertex_point_to_plane_quadrics.h

@@ -43,7 +43,10 @@ namespace igl
     const Eigen::MatrixXi & EF,
     const Eigen::MatrixXi & EI,
     std::vector<
-      std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> > & quadrics);
+      std::tuple<
+        Eigen::MatrixXd,
+        Eigen::RowVectorXd,
+        double> > & quadrics);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "per_vertex_point_to_plane_quadrics.cpp"

include/igl/qslim.cpp

@@ -18,17 +18,26 @@
 #include "qslim_optimal_collapse_edge_callbacks.h"
 #include "quadric_binary_plus_operator.h"
 #include "remove_unreferenced.h"
+#include "intersection_blocking_collapse_edge_callbacks.h"
+#include "AABB.h"
 
 IGL_INLINE bool igl::qslim(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,
-  const size_t max_m,
+  const int max_m,
+  const bool block_intersections,
   Eigen::MatrixXd & U,
   Eigen::MatrixXi & G,
   Eigen::VectorXi & J,
   Eigen::VectorXi & I)
 {
   using namespace igl;
+  igl::AABB<Eigen::MatrixXd, 3> * tree = nullptr;
+  if(block_intersections)
+  {
+    tree = new igl::AABB<Eigen::MatrixXd, 3>();
+    tree->init(V,F);
+  }
 
   // Original number of faces
   const int orig_m = F.rows();
@@ -62,6 +71,13 @@ IGL_INLINE bool igl::qslim(
   decimate_post_collapse_callback      post_collapse;
   qslim_optimal_collapse_edge_callbacks(
     E,quadrics,v1,v2, cost_and_placement, pre_collapse,post_collapse);
+  if(block_intersections)
+  {
+    igl::intersection_blocking_collapse_edge_callbacks(
+      pre_collapse, post_collapse, // These will get copied as needed
+      tree,
+      pre_collapse, post_collapse);
+  }
   // Call to greedy decimator
   bool ret = decimate(
     VO, FO,
@@ -80,5 +96,7 @@ IGL_INLINE bool igl::qslim(
   igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_1,I2);
   I = I(I2).eval();
 
+  assert(tree == nullptr || tree == tree->root());
+  delete tree;
   return ret;
 }

include/igl/qslim.h

@@ -23,6 +23,8 @@ namespace igl
   ///     infinity" in a single triangle.
   /// @param[in] F  #F by 3 list of triangle indices into V
   /// @param[in] max_m  desired number of output faces
+  /// @param[in] block_intersections  whether to block intersections (see
+  ///   intersection_blocking_collapse_edge_callbacks)
   /// @param[out] U  #U by dim list of output vertex posistions (can be same ref as V)
   /// @param[out] G  #G by 3 list of output face indices into U (can be same ref as F)
   /// @param[out] J  #G list of indices into F of birth face
@@ -30,7 +32,8 @@ namespace igl
   IGL_INLINE bool qslim(
     const Eigen::MatrixXd & V,
     const Eigen::MatrixXi & F,
-    const size_t max_m,
+    const int max_m,
+    const bool block_intersections,
     Eigen::MatrixXd & U,
     Eigen::MatrixXi & G,
     Eigen::VectorXi & J,

include/igl/qslim_optimal_collapse_edge_callbacks.h

@@ -35,7 +35,11 @@ namespace igl
   /// \see collapse_edge
   IGL_INLINE void qslim_optimal_collapse_edge_callbacks(
     Eigen::MatrixXi & E,
-    std::vector<std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> > & 
+    std::vector<
+      std::tuple<
+        Eigen::MatrixXd,
+        Eigen::RowVectorXd,
+        double> > & 
       quadrics,
     int & v1,
     int & v2,

include/igl/quad_edges.cpp

@@ -0,0 +1,24 @@
+#include "quad_edges.h"
+#include "unique_simplices.h"
+
+template <
+  typename DerivedQ,
+  typename DerivedE >
+IGL_INLINE void igl::quad_edges(
+  const Eigen::MatrixBase<DerivedQ> & Q,
+  Eigen::PlainObjectBase<DerivedE> & E)
+{
+  E.resize(4*Q.rows(),2);
+  E <<
+    Q.col(0), Q.col(1),
+    Q.col(1), Q.col(2),
+    Q.col(2), Q.col(3),
+    Q.col(3), Q.col(0);
+  igl::unique_simplices(Eigen::MatrixXi(E),E);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::quad_edges<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/quad_edges.h

@@ -0,0 +1,25 @@
+#ifndef IGL_QUAD_EDGES_H
+#define IGL_QUAD_EDGES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Compute the edges of a quad mesh.
+  ///
+  /// @param[in] Q  #Q by 4 list of quad indices into rows of some vertex list V
+  /// @param[out] E  #E by 2 list of edge indices into rows of V
+  template <
+    typename DerivedQ,
+    typename DerivedE >
+  IGL_INLINE void quad_edges(
+    const Eigen::MatrixBase<DerivedQ> & Q,
+    Eigen::PlainObjectBase<DerivedE> & E);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "quad_edges.cpp"
+#endif
+
+#endif
+

include/igl/ray_mesh_intersect.cpp

@@ -119,6 +119,8 @@ IGL_INLINE bool ray_mesh_intersect(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template bool igl::ray_mesh_intersect<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1> const, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+// generated by autoexplicit.sh
 template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, 1, 3, 1, 1, 3>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1> const, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<int, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
 // generated by autoexplicit.sh
 template bool igl::ray_mesh_intersect<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);

include/igl/ray_triangle_intersect.cpp

@@ -0,0 +1,97 @@
+#include "ray_triangle_intersect.h"
+#include <Eigen/Geometry>
+#include <type_traits>
+
+namespace
+{
+  // It'd be nice to have a varadic version of this and put in the igl namespace
+  template <typename Derived>
+  void assert_3_vector(const Eigen::MatrixBase<Derived>& mat) {
+      static_assert(
+          (Derived::RowsAtCompileTime == 3 && Derived::ColsAtCompileTime == 1) ||
+          (Derived::RowsAtCompileTime == 1 && Derived::ColsAtCompileTime == 3) ||
+          (Derived::RowsAtCompileTime == Eigen::Dynamic && Derived::ColsAtCompileTime == 1) ||
+          (Derived::RowsAtCompileTime == 1 && Derived::ColsAtCompileTime == Eigen::Dynamic),
+          "Matrix must be 3x1, 1x3, Eigen::Dynamic×1, or 1×Eigen::Dynamic.");
+  
+      // Add runtime check for dynamic sizes
+      assert((mat.rows() == 3 && mat.cols() == 1) ||
+             (mat.rows() == 1 && mat.cols() == 3));
+  }
+}
+
+template <
+  typename DerivedO,
+  typename DerivedD,
+  typename DerivedV0,
+  typename DerivedV1,
+  typename DerivedV2>
+IGL_INLINE bool igl::ray_triangle_intersect(
+   const Eigen::MatrixBase<DerivedO>  & _O,
+   const Eigen::MatrixBase<DerivedD>  & _D,
+   const Eigen::MatrixBase<DerivedV0> & _V0,
+   const Eigen::MatrixBase<DerivedV1> & _V1,
+   const Eigen::MatrixBase<DerivedV2> & _V2,
+   const typename DerivedO::Scalar epsilon,
+   typename DerivedO::Scalar & t,
+   typename DerivedO::Scalar & u,
+   typename DerivedO::Scalar & v,
+   bool & parallel)
+{
+  // Eigen grossness to ensure we have compile-time 3 vectors below.
+  assert_3_vector(_O);
+  assert_3_vector(_D);
+  assert_3_vector(_V0);
+  assert_3_vector(_V1);
+  assert_3_vector(_V2);
+  using Scalar = typename DerivedO::Scalar;
+  static_assert(std::is_same<Scalar, typename DerivedD::Scalar>::value,
+    "O and D must have the same scalar type");
+  static_assert(std::is_same<Scalar, typename DerivedV0::Scalar>::value,
+    "O and V0 must have the same scalar type");
+  static_assert(std::is_same<Scalar, typename DerivedV1::Scalar>::value,
+    "O and V1 must have the same scalar type");
+  static_assert(std::is_same<Scalar, typename DerivedV2::Scalar>::value,
+    "O and V2 must have the same scalar type");
+  const auto  O = _O.template head<3>();
+  const auto  D = _D.template head<3>();
+  const auto V0 = _V0.template head<3>();
+  const auto V1 = _V1.template head<3>();
+  const auto V2 = _V2.template head<3>();
+  
+  // This Eigen rewrite of raytri.c seems… way faster? like 8× over
+  // intersect_triangle and 50× over intersect_triangle1,2,3. It's not clear to
+  // me why.
+
+  const auto edge1 = (V1-V0).eval();
+  const auto edge2 = (V2-V0).eval();
+  const auto pvec = D.cross(edge2).eval();
+  const Scalar det = edge1.dot(pvec);
+  if( det > -epsilon && det < epsilon )
+  {
+    parallel = true;
+    return false;
+  }
+  parallel = false;
+  const Scalar inv_det = 1.0 / det;
+  const auto tvec = (O - V0).eval();
+  u = tvec.dot(pvec) * inv_det;
+  if( u < 0.0-epsilon || u > 1.0+epsilon )
+  {
+    return false;
+  }
+  const auto qvec = tvec.cross(edge1).eval();
+  v = D.dot(qvec) * inv_det;
+  if( v < 0.0-epsilon || u + v > 1.0+epsilon )
+  {
+    return false;
+  }
+  t = edge2.dot(qvec) * inv_det;
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::ray_triangle_intersect<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&, bool&);
+#endif

include/igl/ray_triangle_intersect.h

@@ -0,0 +1,50 @@
+#ifndef IGL_RAY_TRIANGLE_INTERSECT_H
+#define IGL_RAY_TRIANGLE_INTERSECT_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Determine whether (and if so where) a ray intersects a triangle at a
+  /// single point.
+  ///
+  /// @param[in] O  3d origin of ray
+  /// @param[in] D  3d direction of ray
+  /// @param[in] V0 3d position of first triangle vertex
+  /// @param[in] V1 3d position of second triangle vertex
+  /// @param[in] V2 3d position of third triangle vertex
+  /// @param[in] epsilon  epsilon for determining whether ray is parallel to
+  ///  triangle
+  /// @param[out] t  distance along ray to intersection (if any)
+  /// @param[out] u  barycentric coordinate of V1 triangle vertex
+  /// @param[out] v  barycentric coordinate of V2 triangle vertex
+  /// @param[out] parallel whether ray was considered parallel to triangle (and
+  ///   if so then will return false)
+  /// @returns true if ray intersects triangle
+  ///
+  ///
+  template <
+    typename DerivedO,
+    typename DerivedD,
+    typename DerivedV0,
+    typename DerivedV1,
+    typename DerivedV2>
+  IGL_INLINE bool ray_triangle_intersect(
+     const Eigen::MatrixBase<DerivedO> & O,
+     const Eigen::MatrixBase<DerivedD> & D,
+     const Eigen::MatrixBase<DerivedV0> & V0,
+     const Eigen::MatrixBase<DerivedV1> & V1,
+     const Eigen::MatrixBase<DerivedV2> & V2,
+     const typename DerivedO::Scalar epsilon,
+     typename DerivedO::Scalar & t,
+     typename DerivedO::Scalar & u,
+     typename DerivedO::Scalar & v,
+     bool & parallel);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "ray_triangle_intersect.cpp"
+#endif
+
+#endif

include/igl/remove_unreferenced.cpp

@@ -98,6 +98,8 @@ IGL_INLINE void igl::remove_unreferenced(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::remove_unreferenced<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -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::remove_unreferenced<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::remove_unreferenced<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/tri_tri_intersect.cpp

@@ -888,6 +888,10 @@ IGL_INLINE bool igl::tri_tri_overlap_test_2d(
 
 
 #ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+template bool igl::tri_tri_intersection_test_3d<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template bool igl::tri_tri_intersection_test_3d<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, -1, 1, 1, -1>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template bool igl::tri_tri_intersection_test_3d<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template bool igl::tri_tri_intersection_test_3d<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 template bool igl::tri_tri_intersection_test_3d<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 #endif

include/igl/triangle_triangle_intersect.cpp

@@ -0,0 +1,135 @@
+#include "triangle_triangle_intersect.h"
+#include "triangle_triangle_intersect_shared_edge.h"
+#include "triangle_triangle_intersect_shared_vertex.h"
+#include "tri_tri_intersect.h"
+#include <Eigen/Geometry>
+#include <stdio.h>
+
+//#define IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
+// CGAL::Epeck
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#warning "🐌🐌🐌🐌🐌🐌🐌🐌 Slow debug mode for igl::triangle_triangle_intersect"
+#endif
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedE,
+  typename DerivedEMAP,
+  typename DerivedEF,
+  typename Derivedp>
+IGL_INLINE bool igl::triangle_triangle_intersect(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  const Eigen::MatrixBase<DerivedEF> & EF,
+  const int f,
+  const int c,
+  const Eigen::MatrixBase<Derivedp> & p,
+  const int g)
+{
+  static_assert(
+    std::is_same<typename DerivedV::Scalar,typename Derivedp::Scalar>::value, 
+    "V and p should have same Scalar type");
+  assert(V.cols() == 3);
+  assert(p.cols() == 3);
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
+  using Kernel = CGAL::Epeck;
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Segment_3<Kernel>  Segment_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3;
+    bool cgal_found_intersection = false;
+    Point_3 Vg[3];
+    Point_3 Vf[3];
+    for(int i = 0;i<3;i++)
+    {
+      Vg[i] = Point_3(V(F(g,i),0),V(F(g,i),1),V(F(g,i),2));
+      if(i == c)
+      {
+        Vf[i] = Point_3(p(0),p(1),p(2));
+      }else
+      {
+        Vf[i] = Point_3(V(F(f,i),0),V(F(f,i),1),V(F(f,i),2));
+      }
+    }
+    Triangle_3 Tg(Vg[0],Vg[1],Vg[2]);
+    Triangle_3 Tf(Vf[0],Vf[1],Vf[2]);
+#endif
+
+  // I'm leaving this debug printing stuff in for a bit until I trust this
+  // better. 
+  constexpr bool stinker = false;
+  //const bool stinker = (f==1492 && g==1554);
+  if(stinker) { printf("👀\n"); }
+  bool found_intersection = false;
+  // So edge opposite F(f,c) is the outer edge.
+  const bool share_edge_opposite_c = [&]()
+  {
+    const int o = EMAP(f + c*F.rows());
+    return (EF(o,0) == f && EF(o,1) == g) || (EF(o,1) == f && EF(o,0) == g);
+  }();
+  const int o = EMAP(f + c*F.rows());
+  // Do they share the edge opposite c?
+  if(share_edge_opposite_c)
+  {
+    if(stinker) { printf("⚠️ shares an edge\n"); }
+    found_intersection = triangle_triangle_intersect_shared_edge(V,F,f,c,p,g,1e-8);
+  }else
+  {
+    if(stinker) { printf("does not share an edge\n"); }
+    // Do they share a vertex?
+    int sf,sg;
+    bool found_shared_vertex = false;
+    for(sf = 0;sf<3;sf++)
+    {
+      if(sf == c){ continue;}
+      for(sg = 0;sg<3;sg++)
+      {
+        if(F(f,sf) == F(g,sg))
+        {
+          found_shared_vertex = true;
+          break;
+        }
+      }
+      if(found_shared_vertex) { break;} 
+    }
+    if(found_shared_vertex)
+    {
+      if(stinker) { printf("⚠️ shares a vertex\n"); }
+      found_intersection = 
+        triangle_triangle_intersect_shared_vertex(V,F,f,sf,c,p,g,sg,1e-14);
+    }else
+    {
+      bool coplanar;
+      Eigen::RowVector3d i1,i2;
+      found_intersection = igl::tri_tri_intersection_test_3d(
+                V.row(F(g,0)).template cast<double>(), 
+                V.row(F(g,1)).template cast<double>(), 
+                V.row(F(g,2)).template cast<double>(),
+                            p.template cast<double>(),
+          V.row(F(f,(c+1)%3)).template cast<double>(),
+          V.row(F(f,(c+2)%3)).template cast<double>(),
+          coplanar,
+          i1,i2);
+      if(stinker) { printf("tri_tri_intersection_test_3d says %s\n",found_intersection?"☠️":"✅"); }
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
+      if(CGAL::do_intersect(Tg,Tf))
+      {
+        cgal_found_intersection = true;
+        printf("  ✅ sure it's anything\n");
+      }
+      assert(found_intersection == cgal_found_intersection);
+#endif
+    }
+  }
+  if(stinker) { printf("%s\n",found_intersection?"☠️":"✅"); }
+  return found_intersection;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template bool igl::triangle_triangle_intersect<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::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, int);
+template bool igl::triangle_triangle_intersect<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::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, int);
+#endif

include/igl/triangle_triangle_intersect.h

@@ -0,0 +1,68 @@
+#ifndef IGL_TRIANGLE_TRIANGLE_INTERSECT_H
+#define IGL_TRIANGLE_TRIANGLE_INTERSECT_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Determine whether two triangles intersect. We consider the `f`th and `g`th
+  /// triangles in `F` indexing rows of `V` for 3D positions, but the `c`th corner
+  /// of the `f`th triangle is replaced by `p`. In matlab, this would be
+  ///
+  /// ```matlab
+  /// Tf = V(F(f,:),:);
+  /// Tf(c,:) = p;
+  /// ```
+  ///
+  /// and 
+  ///
+  /// ```matlab
+  /// Tg = V(F(g,:),:);
+  /// ```
+  ///
+  /// Triangles can share an edge, but only if it's the one opposite the replaced
+  /// corner. 
+  ///
+  /// @param[in] V  #V by 3 list of vertex positions
+  /// @param[in] F  #F by 3 list of triangle indices into rows of V
+  /// @param[in] E #E by 2 list of unique undirected edge indices into rows of V
+  /// @param[in] EMAP  #F*3 list of indices into F, mapping each directed edge to
+  ///   unique edge in {1,...,E}
+  /// @param[in] EF  #E by 2 list of edge indices into F
+  /// @param[in] f  index into F of first triangle
+  /// @param[in] c  index into F of corner of first triangle to replace with `p`
+  /// @param[in] p  3D position to replace corner of first triangle
+  /// @param[in] g  index into F of second triangle
+  /// @returns  true if triangles intersect
+  ///
+  /// \note This very specialized function should be complemented with the more
+  /// general functions in igl::tri_tri_intersect (which should be renamed and
+  /// split into appropriate files). There's a reasonably amount of shared code
+  /// with igl::copyleft::cgal::remesh_self_intersections too.
+  ///
+  /// \see edge_flaps, tri_tri_intersect
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedE,
+    typename DerivedEMAP,
+    typename DerivedEF,
+    typename Derivedp>
+  IGL_INLINE bool triangle_triangle_intersect(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+    const Eigen::MatrixBase<DerivedEF> & EF,
+    const int f,
+    const int c,
+    const Eigen::MatrixBase<Derivedp> & p,
+    const int g);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangle_triangle_intersect.cpp"
+#endif
+
+#endif

include/igl/triangle_triangle_intersect_shared_edge.cpp

@@ -0,0 +1,116 @@
+#include "triangle_triangle_intersect_shared_edge.h"
+#include "PI.h"
+#include <Eigen/Geometry>
+#include <type_traits>
+
+//#define IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_EDGE_DEBUG
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_EDGE_DEBUG
+// CGAL::Epeck
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#warning "🐌🐌🐌🐌🐌🐌🐌🐌 Slow debug mode for igl::triangle_triangle_intersect"
+#endif
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedp>
+IGL_INLINE bool igl::triangle_triangle_intersect_shared_edge(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const int f,
+  const int c,
+  const Eigen::MatrixBase<Derivedp> & p,
+  const int g,
+  const typename DerivedV::Scalar epsilon)
+{
+  constexpr bool stinker = false;
+  static_assert(
+    std::is_same<typename DerivedV::Scalar,typename Derivedp::Scalar>::value, 
+    "V and p should have same Scalar type");
+  assert(V.cols() == 3);
+  assert(p.cols() == 3);
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
+  using Kernel = CGAL::Epeck;
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Segment_3<Kernel>  Segment_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3;
+    bool cgal_found_intersection = false;
+    Point_3 Vg[3];
+    Point_3 Vf[3];
+    for(int i = 0;i<3;i++)
+    {
+      Vg[i] = Point_3(V(F(g,i),0),V(F(g,i),1),V(F(g,i),2));
+      if(i == c)
+      {
+        Vf[i] = Point_3(p(0),p(1),p(2));
+      }else
+      {
+        Vf[i] = Point_3(V(F(f,i),0),V(F(f,i),1),V(F(f,i),2));
+      }
+    }
+    Triangle_3 Tg(Vg[0],Vg[1],Vg[2]);
+    Triangle_3 Tf(Vf[0],Vf[1],Vf[2]);
+#endif
+  bool found_intersection = false;
+  // Only intersects if the dihedral angle is zero (precondition: no zero
+  // area triangles before or after collapse)
+
+  // normal of g
+  const auto vg10 = (V.row(F(g,1))-V.row(F(g,0))).template head<3>();
+  const auto vg20 = (V.row(F(g,2))-V.row(F(g,0))).template head<3>();
+  const auto ng = vg10.cross(vg20);
+
+  // normal of f (with p replaced)
+  const auto vf1p = (V.row(F(f,(c+1)%3))-p).template head<3>();
+  const auto vf2p = (V.row(F(f,(c+2)%3))-p).template head<3>();
+  const auto nf = vf1p.cross(vf2p);
+
+  // normalized edge vector
+  const auto o_vec_un = (V.row(F(f,(c+2)%3))-V.row(F(f,(c+1)%3))).template head<3>();
+  const auto o_vec = o_vec_un.stableNormalized();
+  const auto theta = std::abs(std::atan2(o_vec.dot(ng.cross(nf)),ng.dot(nf)));
+  if(stinker){ printf(" theta: %g\n",theta); }
+  // This magic number should be a parameter
+  if(theta < igl::PI - epsilon)
+  {
+    return false;
+  }
+  // Triangles really really might intersect.
+  found_intersection = true;
+  // Find intersection
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_EDGE_DEBUG
+  if(CGAL::do_intersect(Tg,Tf))
+  {
+    CGAL::Object obj = CGAL::intersection(Tg,Tf);
+    if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
+    {
+      printf("  ❌ segment is just the edge.\n");
+    }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
+    {
+      printf("  🤔 how just a single point?\n");
+    } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj))
+    {
+      cgal_found_intersection = true;
+      printf("  ✅ sure it's a triangle\n");
+    } else if(const std::vector<Point_3 > *polyp =
+        CGAL::object_cast< std::vector<Point_3 > >(&obj))
+    {
+      printf("  ✅ polygon\n");
+    }else {
+      printf("  🤔 da fuke?\n");
+    }
+  }
+  assert(found_intersection == cgal_found_intersection);
+#endif
+  return found_intersection;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::triangle_triangle_intersect_shared_edge<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar);
+// generated by autoexplicit.sh
+template bool igl::triangle_triangle_intersect_shared_edge<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar);
+// generated by autoexplicit.sh
+template bool igl::triangle_triangle_intersect_shared_edge<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar);
+#endif

include/igl/triangle_triangle_intersect_shared_edge.h

@@ -0,0 +1,44 @@
+#ifndef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_EDGE_H
+#define IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_EDGE_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Determine whether two triangles --- which share an edge--- intersect. We
+  /// consider the `f`th and `g`th triangles in `F` indexing rows of `V` for 3D
+  /// positions, but the `c`th corner (opposite the shared edge) of the `f`th
+  /// triangle is replaced by `p`. 
+  ///
+  /// @param[in] V  #V by 3 list of vertex positions
+  /// @param[in] F  #F by 3 list of triangle indices into rows of V
+  /// @param[in] f  index into F of first triangle
+  /// @param[in] c  index into F of corner opposite shared edge
+  /// @param[in] p  3D position to replace cth corner of first triangle
+  /// @param[in] g  index into F of second triangle
+  /// @param[in] epsilon  tolerance used to determine coplanarity
+  /// @returns  true if triangles intersect
+  ///
+  /// \see edge_flaps, tri_tri_intersect, triangle_triangle_intersect
+  ///
+  /// \pre both faces are assumed to have non-trivial area
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Derivedp>
+  IGL_INLINE bool triangle_triangle_intersect_shared_edge(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const int f,
+    const int c,
+    const Eigen::MatrixBase<Derivedp> & p,
+    const int g,
+    const typename DerivedV::Scalar epsilon);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangle_triangle_intersect_shared_edge.cpp"
+#endif
+
+#endif

include/igl/triangle_triangle_intersect_shared_vertex.cpp

@@ -0,0 +1,335 @@
+#include "triangle_triangle_intersect_shared_vertex.h"
+#include "ray_triangle_intersect.h"
+#include "barycentric_coordinates.h"
+#include "matlab_format.h"
+#include <Eigen/Geometry>
+#include <iostream>
+#include <iomanip>
+#include <stdio.h>
+// std::signbit
+#include <cmath>
+
+//#define IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_VERTEX_DEBUG
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_VERTEX_DEBUG
+// CGAL::Epeck
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#warning "🐌🐌🐌🐌🐌🐌🐌🐌 Slow debug mode for igl::triangle_triangle_intersect"
+#endif
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedp>
+IGL_INLINE bool igl::triangle_triangle_intersect_shared_vertex(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const int f,
+  const int sf,
+  const int c,
+  const Eigen::MatrixBase<Derivedp> & p,
+  const int g,
+  const int sg,
+  const typename DerivedV::Scalar epsilon)
+{
+  static_assert(
+    std::is_same<typename DerivedV::Scalar,typename Derivedp::Scalar>::value, 
+    "V and p should have same Scalar type");
+  assert(V.cols() == 3);
+  assert(p.cols() == 3);
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
+  using Kernel = CGAL::Epeck;
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Segment_3<Kernel>  Segment_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3;
+    bool cgal_found_intersection = false;
+    Point_3 Vg[3];
+    Point_3 Vf[3];
+    for(int i = 0;i<3;i++)
+    {
+      Vg[i] = Point_3(V(F(g,i),0),V(F(g,i),1),V(F(g,i),2));
+      if(i == c)
+      {
+        Vf[i] = Point_3(p(0),p(1),p(2));
+      }else
+      {
+        Vf[i] = Point_3(V(F(f,i),0),V(F(f,i),1),V(F(f,i),2));
+      }
+    }
+    Triangle_3 Tg(Vg[0],Vg[1],Vg[2]);
+    Triangle_3 Tf(Vf[0],Vf[1],Vf[2]);
+#endif
+  constexpr bool stinker = false;
+  bool found_intersection = false;
+  // If they share a vertex and intersect, then an opposite edge must
+  // stab through the other triangle.
+
+  // Using ray_triangle_intersect now, not sure we need this copy or cast
+  Eigen::RowVector3d g0 = V.row(F(g,0)).template cast<double>();
+  Eigen::RowVector3d g1 = V.row(F(g,1)).template cast<double>();
+  Eigen::RowVector3d g2 = V.row(F(g,2)).template cast<double>();
+  Eigen::RowVector3d fs;
+  if(((sf+1)%3)==c)
+  {
+    fs = p;
+  }else
+  {
+    fs = V.row(F(f,(sf+1)%3));
+  }
+  Eigen::RowVector3d fd;
+  if( ((sf+2)%3)==c )
+  {
+    fd = p.template cast<double>();
+  }else
+  {
+    fd = V.row(F(f,(sf+2)%3)).template cast<double>();
+  }
+  Eigen::RowVector3d fdir = fd - fs;
+  double t,u,v;
+
+  if(stinker)
+  {
+    std::cout<<"T = ["<<g0<<";" <<g1<<";"<<g2<<"];"<<std::endl;
+    std::cout<<"src = [" <<fs<<"];"<<std::endl;
+    std::cout<<"dir = [" <<fdir<<"];"<<std::endl;
+  }
+  // p = (1-u-v)*a + u*b + v*c
+  const auto bary = [](
+      const Eigen::RowVector3d & p,
+      const Eigen::RowVector3d & a,
+      const Eigen::RowVector3d & b,
+      const Eigen::RowVector3d & c,
+      double & u,
+      double & v)
+  {
+    const auto v0 = (b-a).eval();
+    const auto v1 = (c-a).eval();
+    const auto v2 = (p-a).eval();
+    const double d00 = v0.dot(v0);
+    const double d01 = v0.dot(v1);
+    const double d11 = v1.dot(v1);
+    const double d20 = v2.dot(v0);
+    const double d21 = v2.dot(v1);
+    const double denom = d00 * d11 - d01 * d01;
+    u = (d11 * d20 - d01 * d21) / denom;
+    v = (d00 * d21 - d01 * d20) / denom;
+    // Equivalent:
+    //Eigen::RowVector3d l;
+    //igl::barycentric_coordinates(p,a,b,c,l);
+    //u = l(1); v = l(2);
+  };
+
+  // Does the segment (A,B) intersect the triangle (0,0),(1,0),(0,1)?
+  const auto intersect_unit_helper = [](
+      const Eigen::RowVector2d & A,
+      const Eigen::RowVector2d & B) -> bool
+  {
+    // Check if P is inside (0,0),(1,0),(0,1) triangle
+    const auto inside_unit = []( const Eigen::RowVector2d & P) -> bool
+    {
+      return P(0) >= 0 && P(1) >= 0 && P(0) + P(1) <= 1;
+    };
+    if(inside_unit(A) || inside_unit(B))
+    { 
+      return true; 
+    }
+
+    const auto open_interval_contains_zero = [](
+        const double a, const double b) -> bool
+    {
+      // handle case where either is 0.0 or -0.0
+      if(a==0 || b==0) { return false; }
+      return std::signbit(a) != std::signbit(b);
+    };
+
+    // Now check if the segment intersects any of the edges.
+    // Does A-B intesect X-axis?
+    if(open_interval_contains_zero(A(1),B(1)))
+    {
+      assert((A(1) - B(1)) != 0);
+      // A and B are on opposite sides of the X-axis
+      const double t = A(1) / (A(1) - B(1));
+      const double x = A(0) + t * (B(0) - A(0));
+      if(x >= 0 && x <= 1)
+      {
+        return true;
+      }
+    }
+    // Does A-B intesect Y-axis?
+    if(open_interval_contains_zero(A(0),B(0)))
+    {
+      assert((A(0) - B(0)) != 0);
+      // A and B are on opposite sides of the Y-axis
+      const double t = A(0) / (A(0) - B(0));
+      const double y = A(1) + t * (B(1) - A(1));
+      if(y >= 0 && y <= 1)
+      {
+        return true;
+      }
+    }
+    // Does A-B intersect the line x+y=1?
+    if(open_interval_contains_zero(A(0) + A(1) - 1.0,B(0) + B(1) - 1.0))
+    {
+      assert((A(0) + A(1) - 1.0) - (B(0) + B(1) - 1.0) != 0);
+      // A and B are on opposite sides of the line x+y=1
+      // x+y=1
+      // A(0) + t * (B(0) - A(0)) + A(1) + t * (B(1) - A(1)) = 1
+      // t * (B(0) - A(0) + B(1) - A(1)) = 1 - A(0) - A(1)
+      const double  t = (1 - A(0) - A(1)) / (B(0) - A(0) + B(1) - A(1));
+      const double y = A(1) + t * (B(1) - A(1));
+      if(y >= 0 && y <= 1)
+      {
+        return true;
+      }
+    }
+    return false;
+  };
+  const auto intersect_unit = [&intersect_unit_helper](
+      const Eigen::RowVector2d & A,
+      const Eigen::RowVector2d & B) -> bool
+  {
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_VERTEX_DEBUG
+    printf("A=[%g,%g];B=[%g,%g];\n",
+        A(0),A(1),B(0),B(1));
+#endif
+    const bool ret = intersect_unit_helper(A,B);
+    return ret;
+  };
+  const auto point_on_plane = [&epsilon](
+      const Eigen::RowVector3d & p,
+      const Eigen::RowVector3d & a,
+      const Eigen::RowVector3d & b,
+      const Eigen::RowVector3d & c) -> bool
+  {
+    const auto n = (b-a).cross(c-a);
+    const auto d = n.dot(p-a);
+    return std::abs(d) < epsilon*n.stableNorm();
+  };
+
+
+  //if(intersect_triangle1(
+  //      fs.data(),fdir.data(),
+  //      g0.data(),g1.data(),g2.data(),
+  //      &t,&u,&v))
+  bool parallel = false;
+  if(ray_triangle_intersect(
+        fs,fdir,
+        g0,g1,g2,
+        epsilon,
+        t,u,v,parallel))
+  {
+    found_intersection = t > 0.0 && t<1.0+epsilon;
+  }else if(parallel)
+  {
+    if(stinker){ printf("    parallel\n"); }
+    if(point_on_plane(fs,g0,g1,g2))
+    {
+      if(stinker){ printf("    coplanar\n"); }
+      // deal with parallel
+      Eigen::RowVector2d s2,d2;
+      bary(fs,g0,g1,g2,s2(0),s2(1));
+      bary(fd,g0,g1,g2,d2(0),d2(1));
+      found_intersection = intersect_unit(s2,d2);
+    }
+  }
+  if(stinker){ printf("    found_intersection: %d\n",found_intersection); }
+
+  if(!found_intersection)
+  {
+    Eigen::RowVector3d fv[3];
+    fv[0] = V.row(F(f,0)).template cast<double>();
+    fv[1] = V.row(F(f,1)).template cast<double>();
+    fv[2] = V.row(F(f,2)).template cast<double>();
+    fv[c] = p.template cast<double>();
+    Eigen::RowVector3d gs = V.row(F(g,(sg+1)%3)).template cast<double>();
+    Eigen::RowVector3d gd = V.row(F(g,(sg+2)%3)).template cast<double>();
+    Eigen::RowVector3d gdir = gd - gs;
+
+    if(stinker)
+    {
+      std::cout<<"T = ["<<fv[0]<<";"<<fv[1]<<";"<<fv[2]<<"];"<<std::endl;
+      std::cout<<"src = [" <<gs<<"];"<<std::endl;
+      std::cout<<"dir = [" <<gdir<<"];"<<std::endl;
+    }
+    if(ray_triangle_intersect(
+          gs,gdir,
+          fv[0],fv[1],fv[2],
+          epsilon,
+          t,u,v,parallel))
+    {
+      found_intersection = t > 0 && t<1+epsilon;
+    }else if(parallel)
+    {
+      if(stinker){ printf("    parallel2\n"); }
+      if(point_on_plane(gs,fv[0],fv[1],fv[2]))
+      {
+        if(stinker){ printf("    coplanar\n"); }
+        // deal with parallel
+        //assert(false);
+        Eigen::RowVector2d s2,d2;
+        bary(gs,fv[0],fv[1],fv[2],s2(0),s2(1));
+        bary(gd,fv[0],fv[1],fv[2],d2(0),d2(1));
+        found_intersection = intersect_unit(s2,d2);
+      }
+    }
+  }
+  if(stinker){ printf("    found_intersection2: %d\n",found_intersection); }
+#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_VERTEX_DEBUG
+  if(CGAL::do_intersect(Tg,Tf))
+  {
+    CGAL::Object obj = CGAL::intersection(Tg,Tf);
+    if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
+    {
+      printf("  ✅ sure it's a segment\n");
+      cgal_found_intersection = true;
+    }else if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
+    {
+      printf("  ❌ it's just the point.\n");
+    } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj))
+    {
+      cgal_found_intersection = true;
+      printf("  ✅ sure it's a triangle\n");
+    } else if(const std::vector<Point_3 > *polyp =
+        CGAL::object_cast< std::vector<Point_3 > >(&obj))
+    {
+      cgal_found_intersection = true;
+      printf("  ✅ polygon\n");
+    }else {
+      printf("  🤔 da fuke?\n");
+    }
+  }
+  printf("%d,%d  %s vs %s\n",f,c,found_intersection?"☠️":"✅",cgal_found_intersection?"☠️":"✅");
+  if(found_intersection != cgal_found_intersection)
+  {
+    printf("Tg = [[%g,%g,%g];[%g,%g,%g];[%g,%g,%g]];\n",
+        CGAL::to_double(Tg.vertex(0).x()),
+        CGAL::to_double(Tg.vertex(0).y()),
+        CGAL::to_double(Tg.vertex(0).z()),
+        CGAL::to_double(Tg.vertex(1).x()),
+        CGAL::to_double(Tg.vertex(1).y()),
+        CGAL::to_double(Tg.vertex(1).z()),
+        CGAL::to_double(Tg.vertex(2).x()),
+        CGAL::to_double(Tg.vertex(2).y()),
+        CGAL::to_double(Tg.vertex(2).z()));
+    printf("Tf = [[%g,%g,%g];[%g,%g,%g];[%g,%g,%g]];\n",
+        CGAL::to_double(Tf.vertex(0).x()),
+        CGAL::to_double(Tf.vertex(0).y()),
+        CGAL::to_double(Tf.vertex(0).z()),
+        CGAL::to_double(Tf.vertex(1).x()),
+        CGAL::to_double(Tf.vertex(1).y()),
+        CGAL::to_double(Tf.vertex(1).z()),
+        CGAL::to_double(Tf.vertex(2).x()),
+        CGAL::to_double(Tf.vertex(2).y()),
+        CGAL::to_double(Tf.vertex(2).z()));
+  }
+  assert(found_intersection == cgal_found_intersection);
+#endif
+  return found_intersection;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::triangle_triangle_intersect_shared_vertex<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, int, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, int, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar);
+template bool igl::triangle_triangle_intersect_shared_vertex<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, int, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, int, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar);
+template bool igl::triangle_triangle_intersect_shared_vertex<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, int, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, int, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar);
+#endif

include/igl/triangle_triangle_intersect_shared_vertex.h

@@ -0,0 +1,49 @@
+#ifndef IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_VERTEX_H
+#define IGL_TRIANGLE_TRIANGLE_INTERSECT_SHARED_VERTEX_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Determine whether two triangles --- which share a vertex F(f,sf) ==
+  /// F(g,sg) --- intersect. We consider the `f`th and `g`th triangles in `F`
+  /// indexing rows of `V` for 3D positions, but the `c`th corner (opposite the
+  /// shared edge) of the `f`th triangle is replaced by `p`. 
+  ///
+  /// @param[in] V  #V by 3 list of vertex positions
+  /// @param[in] F  #F by 3 list of triangle indices into rows of V
+  /// @param[in] f  index into F of first triangle
+  /// @param[in] sf  corner index of shared vertex in first triangle
+  /// @param[in] c  index into F of corner opposite shared edge (assumed c≠sf)
+  /// @param[in] p  3D position to replace cth corner of first triangle
+  /// @param[in] g  index into F of second triangle
+  /// @param[in] sg  corner index of shared vertex in second triangle
+  /// @param[in] epsilon  tolerance used to determine intersection
+  /// @returns  true if triangles intersect
+  ///
+  /// \see edge_flaps, tri_tri_intersect, triangle_triangle_intersect,
+  /// triangle_triangle_intersect_shared_edge
+  ///
+  /// \pre both faces are assumed to have non-trivial area
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Derivedp>
+  IGL_INLINE bool triangle_triangle_intersect_shared_vertex(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const int f,
+    const int sf,
+    const int c,
+    const Eigen::MatrixBase<Derivedp> & p,
+    const int g,
+    const int sg,
+    const typename DerivedV::Scalar epsilon);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangle_triangle_intersect_shared_vertex.cpp"
+#endif
+
+#endif

include/igl/unproject_onto_mesh.cpp

@@ -73,6 +73,8 @@ IGL_INLINE bool igl::unproject_onto_mesh(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template bool igl::unproject_onto_mesh<Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, std::function<bool (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, igl::Hit&)> const&, int&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);
+// generated by autoexplicit.sh
 template bool igl::unproject_onto_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
 // generated by autoexplicit.sh
 template bool igl::unproject_onto_mesh<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, int&, Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> >&);

tests/include/igl/AABB.cpp

@@ -1,5 +1,15 @@
-#include <igl/AABB.h>
 #include <test_common.h>
+#include <igl/AABB.h>
+#include <igl/EPS.h>
+#include <igl/avg_edge_length.h>
+#include <igl/barycenter.h>
+#include <igl/colon.h>
+#include <igl/get_seconds.h>
+#include <igl/point_mesh_squared_distance.h>
+#include <igl/point_simplex_squared_distance.h>
+#include <igl/randperm.h>
+#include <igl/read_triangle_mesh.h>
+#include <iostream>
 
 TEST_CASE("AABB: find_2d", "[igl]")
 {
@@ -52,3 +62,171 @@ TEST_CASE("AABB: find_3d", "[igl]")
   REQUIRE(r[0] == 0);
   REQUIRE(r[1] == 1);
 }
+
+TEST_CASE("AABB: insert", "[igl]")
+{
+  igl::AABB<Eigen::MatrixXd, 3> * tree = new igl::AABB<Eigen::MatrixXd, 3>();
+  for(int i = 0;i<3;i++)
+  {
+    igl::AABB<Eigen::MatrixXd, 3> * leaf = new igl::AABB<Eigen::MatrixXd, 3>();
+    leaf->m_box = Eigen::AlignedBox<double, 3>(
+      Eigen::RowVector3d(-i,-i,-i),
+      Eigen::RowVector3d(i,i,i));
+    auto * ret = tree->insert(leaf);
+    tree = ret->root();
+    REQUIRE(leaf->is_leaf());
+  }
+  REQUIRE(tree->is_root());
+  delete tree;
+}
+
+TEST_CASE("AABB: dynamic", "[igl]")
+{
+  const int MAX_RUNS = 10;
+  const auto test_case = [&MAX_RUNS](const std::string &param)
+  {
+    IGL_TICTOC_LAMBDA;
+    Eigen::MatrixXd V;
+    Eigen::MatrixXi F;
+    // Load example mesh: GetParam() will be name of mesh file
+    igl::read_triangle_mesh(test_common::data_path(param), V, F);
+    // Make into soup
+    V = V(Eigen::Map<Eigen::VectorXi>(F.data(),F.size()), Eigen::all).eval();
+    F = Eigen::Map<Eigen::MatrixXi>(igl::colon<int>(0,V.rows()-1).data(),V.rows()/3,3).eval();
+    Eigen::MatrixXd BC;
+    igl::barycenter(V,F,BC);
+    Eigen::VectorXi I_gt = igl::colon<int>(0,F.rows()-1);
+    Eigen::VectorXd sqrD_gt = Eigen::VectorXd::Zero(F.rows());
+    Eigen::VectorXd sqrD;
+    Eigen::VectorXi I;
+    Eigen::MatrixXd C;
+    igl::AABB<Eigen::MatrixXd, 3> * tree = new igl::AABB<Eigen::MatrixXd, 3>();
+    tree->init(V,F);
+    tictoc();
+    tree->squared_distance(V,F,BC,sqrD,I,C);
+    const double t0 = tictoc();
+    test_common::assert_eq(I,I_gt);
+    test_common::assert_near(sqrD,sqrD_gt,1e-15);
+    REQUIRE(tree->size() == 2*F.rows()-1);
+
+    const double pad = igl::EPS<double>();
+    const double h = igl::avg_edge_length(V,F);
+    const auto leaves = tree->gather_leaves(F.rows());
+    for(auto * leaf : leaves)
+    {
+      REQUIRE(leaf);
+      REQUIRE(leaf->is_leaf());
+    }
+    // This is slow. Only test on small meshes
+    if(F.rows() < 4000)
+    {
+      // Gather list of pointers to leaves and pad by ε to force rebuild
+      tree = tree->pad(leaves,pad,2);
+      REQUIRE(tree->is_root());
+      REQUIRE(tree == tree->root());
+      for(auto * leaf : leaves)
+      {
+        REQUIRE(leaf);
+        REQUIRE(leaf->is_leaf());
+      }
+      tictoc();
+      tree->squared_distance(V,F,BC,sqrD,I,C);
+      const double t1 = tictoc();
+      REQUIRE(t1 < 10*t0);
+      test_common::assert_eq(I,I_gt);
+      test_common::assert_near(sqrD,sqrD_gt,1e-15);
+      REQUIRE(tree->size() == 2*F.rows()-1);
+    }
+#ifndef NDEBUG
+    if(F.rows()>10000)
+    {
+      std::cerr<<"#ifndef NDEBUG: Skipping timing test."<<std::endl;
+      delete tree;
+      return;
+    }
+#endif
+
+
+    Eigen::VectorXi RV;
+    Eigen::VectorXi RF;
+    // Perturb a small subset of the triangles
+    {
+      {
+        igl::randperm(F.rows(),RF);
+        RF = RF.topRows(std::min(12,(int)F.rows())).eval();
+        RV.resize(RF.size()*3);
+        RV << RF, RF.array()+F.rows(), RF.array()+2*F.rows();
+      }
+      Eigen::MatrixXd TF = 0.1*h*Eigen::MatrixXd::Random(RF.size(),3);
+      Eigen::MatrixXd TV(RV.rows(),3);
+      TV<<TF,TF,TF;
+      V(RV,Eigen::all) += TV;
+      igl::barycenter(V,F,BC);
+    }
+    const int qi = RF(0);
+
+    double t_dynamic_tree;
+    {
+      tictoc();
+      // update those perturbed triangles
+      for(int i = 0;i<RF.size();i++)
+      {
+        tree = leaves[RF(i)]->update_primitive(V,F,pad)->root();
+      }
+      const double t_refit = tictoc();
+      for(int r = 0;r<MAX_RUNS;r++)
+      {
+        tree->squared_distance(V,F,Eigen::MatrixXd(BC.row(qi)),sqrD,I,C);
+      }
+      const double t_query = tictoc()/MAX_RUNS;
+      REQUIRE(I(0) == qi);
+      t_dynamic_tree = t_refit+t_query;
+    }
+    double t_static_tree;
+    {
+      tictoc();
+      for(int r = 0;r<MAX_RUNS;r++)
+      {
+        igl::point_mesh_squared_distance(
+          Eigen::MatrixXd(BC.row(qi)),V,F,sqrD,I,C);
+      }
+      t_static_tree = tictoc()/MAX_RUNS;
+      REQUIRE(I(0) == qi);
+    }
+    double t_brute_force;
+    {
+      tictoc();
+      double min_dist;
+      int min_i;
+      for(int r = 0;r<MAX_RUNS;r++)
+      {
+        min_dist = std::numeric_limits<double>::infinity();
+        min_i = -1;
+        Eigen::RowVector3d q = Eigen::RowVector3d(BC.row(qi));
+        for(int i = 0;i<F.rows();i++)
+        {
+          Eigen::RowVector3d c;
+          double d;
+          igl::point_simplex_squared_distance<3>(q,V,F,i,d,c);
+          if(d < min_dist)
+          {
+            min_dist = d;
+            min_i = i;
+          }
+        }
+      }
+      t_brute_force = tictoc()/MAX_RUNS;
+      REQUIRE(min_i == qi);
+    }
+
+    // Only compare speeds on large meshes
+    if(F.rows() > 300)
+    { 
+      REQUIRE(t_dynamic_tree < t_static_tree);
+      REQUIRE(t_dynamic_tree < t_brute_force);
+    }
+    delete tree;
+  };
+
+  test_common::run_test_cases(test_common::all_meshes(), test_case);
+}

tests/include/igl/decimate.cpp

@@ -24,7 +24,7 @@ TEST_CASE("decimate: hemisphere", "[igl]")
   // Perfect normals from positions
   Eigen::MatrixXd NV = V.rowwise().normalized();
   // Remove half of the faces
-  igl::decimate(V,F,F.rows()/2,U,G,J,I);
+  igl::decimate(V,F,F.rows()/2,false,U,G,J,I);
   // Expect that all normals still point in same direction as original
   Eigen::MatrixXd NU = U.rowwise().normalized();
   Eigen::MatrixXd NVI = NV(I,Eigen::all);
@@ -43,19 +43,19 @@ TEST_CASE("decimate: closed", "[igl]")
   {
     Eigen::MatrixXd V,U;
     Eigen::MatrixXi F,G;
-    Eigen::VectorXi J;
+    Eigen::VectorXi I,J;
     // Load example mesh: GetParam() will be name of mesh file
     igl::read_triangle_mesh(test_common::data_path(param), V, F);
-    igl::decimate(V,F,0,U,G,J);
+    igl::decimate(V,F,0,false,U,G,I,J);
     REQUIRE (4 == U.rows());
     REQUIRE (4 == G.rows());
     {
-      Eigen::MatrixXi I;
-      igl::sort(Eigen::MatrixXi(G),2,true,G,I);
+      Eigen::MatrixXi _;
+      igl::sort(Eigen::MatrixXi(G),2,true,G,_);
     }
     {
-      Eigen::VectorXi I;
-      igl::sortrows(Eigen::MatrixXi(G),true,G,I);
+      Eigen::VectorXi _;
+      igl::sortrows(Eigen::MatrixXi(G),true,G,_);
     }
     // Tet with sorted faces
     Eigen::MatrixXi T(4,3);

tests/include/igl/fast_find_intersections.cpp

@@ -1,78 +1,14 @@
-#include <test_common.h>
-#include <igl/tri_tri_intersect.h>
+#include "test_common.h"
 #include <igl/fast_find_intersections.h>
 #include <igl/fast_find_self_intersections.h>
 #include <igl/combine.h>
+#include <igl/unique.h>
+#include <igl/matlab_format.h>
+#include <iostream>
 
-TEST_CASE("tri_tri_intersection_test_3d intersect", "[igl]")
-{
-  // try with oblique interecion plane
-  for(double shift=-1000;shift<=1000;shift+=100.0)
-  {
-    Eigen::RowVector3d p1(0,0,0),    q1(1,0,0),r1(0,1,0);
-    Eigen::RowVector3d p2(shift,0,1),q2(1,1,0),r2(-shift,0,-1);
 
-    // should intersect along the vector (0,0,0) -> (0.5,0.5,0)
-    Eigen::RowVector3d s,t;
-    bool coplanar;
-    REQUIRE( igl::tri_tri_intersection_test_3d(p1,q1,r1, p2,q2,r2, coplanar, s, t) );
-    REQUIRE( !coplanar);
-    
-    if(s.norm()<1e-5)
-    {
-      Eigen::RowVector3d t_correct(0.5,0.5,0);
-      Eigen::RowVector3d s_correct(0,0,0);
-      test_common::assert_near( t, t_correct, 1e-10);
-      test_common::assert_near( s, s_correct, 1e-10);
-    } else {
-      Eigen::RowVector3d s_correct(0.5,0.5,0);
-      Eigen::RowVector3d t_correct(0,0,0);
-      test_common::assert_near( t, t_correct, 1e-10);
-      test_common::assert_near( s, s_correct, 1e-10);
-    }
-  }
-}
 
-TEST_CASE("tri_tri_intersection_test_3d not intersect", "[igl]")
-{
-  // positive test that two triangles intersect
-  Eigen::RowVector3d p1(0,0,0),q1(1,0,0),r1(0,1,0);
-  Eigen::RowVector3d p2(0,0,1),q2(1,0,1),r2(0,1,1);
-
-  // should intersect along the vector (0,0,0) -> (sqrt(2),sqrt(2),0)
-  Eigen::RowVector3d s,t;
-  bool coplanar;
-  REQUIRE( !igl::tri_tri_intersection_test_3d(p1,q1,r1,p2,q2,r2, coplanar, s, t) );
-}
-
-
-TEST_CASE("tri_tri_intersection_test_3d coplanar", "[igl]")
-{
-  // positive test that two triangles intersect
-  Eigen::RowVector3d p1(0,0,0),q1(1,0,0),r1(0,1,0);
-  Eigen::RowVector3d p2(0,0,0),q2(0.5,0,0),r2(0,0.5,0);
-
-  // should intersect along the vector (0,0,0) -> (sqrt(2),sqrt(2),0)
-  Eigen::RowVector3d s,t;
-  bool coplanar;
-  REQUIRE( igl::tri_tri_intersection_test_3d(p1,q1,r1,p2,q2,r2, coplanar, s, t) );
-  REQUIRE(coplanar);
-}
-
-
-TEST_CASE("fast_find_self_intersections: negative", "[igl]")
-{
-  Eigen::MatrixXd V;
-  Eigen::MatrixXi F;
-  Eigen::VectorXi I;
-
-  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
-  REQUIRE (! igl::fast_find_self_intersections(V,F,I) );
-
-  REQUIRE ( I.sum()==0);
-}
-
-TEST_CASE("fast_find_self_intersections: positive", "[igl]")
+TEST_CASE("fast_find_intersections: cube-triangle", "[igl]")
 {
   Eigen::MatrixXd V;
   Eigen::MatrixXi F;
@@ -90,53 +26,28 @@ TEST_CASE("fast_find_self_intersections: positive", "[igl]")
   Eigen::MatrixXi Fp(1,3);
   Fp << 0,1,2;
   
-  Eigen::MatrixXd V_;
-  Eigen::MatrixXi F_;
-  igl::combine<Eigen::MatrixXd,Eigen::MatrixXi>({V,Vp},{F,Fp}, V_,F_);
-
-  Eigen::VectorXi I;
-  Eigen::MatrixXd edges;
-
-  REQUIRE ( igl::fast_find_self_intersections(V_,F_,I,edges) );
-  // should have 9 triangles that are intersecting each other
-  REQUIRE ( I.sum()==9);
-  // and 16 line edges
-  REQUIRE ( edges.rows()==16);
-  // plane that we added intersects other triangles
-  REQUIRE ( I(F_.rows()-1)!=0 );
-  // TODO: check if the edges are at the right place (?) 
-}
 
-TEST_CASE("fast_find_intersections:", "[igl]")
-{
-  Eigen::MatrixXd V;
-  Eigen::MatrixXi F;
-
-  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
-
-  auto dims=(V.colwise().maxCoeff()-V.colwise().minCoeff())*2;
-  auto cz=(V.col(2).minCoeff()+V.col(2).maxCoeff())/2;
-
-  // add a triangle intersecting cube 
-  Eigen::MatrixXd Vp(3,3);
-  Vp << V.col(0).minCoeff()-dims(0),   V.col(1).minCoeff()-dims(1),   cz,
-        V.col(0).minCoeff()-dims(0),   V.col(1).maxCoeff()+dims(1)*3, cz,
-        V.col(0).maxCoeff()+dims(0)*3, V.col(1).maxCoeff()-dims(1),   cz;
-  Eigen::MatrixXi Fp(1,3);
-  Fp << 0,1,2;
-  
-  Eigen::MatrixXi I;
-  Eigen::MatrixXd edges;
-
-  igl::fast_find_intersections(V,F,Vp,Fp,I,edges);
-
-  // should have 8 triangles that are intersecting plane
-  REQUIRE ( I.rows()==8);
+  Eigen::MatrixXi IF,EE;
+  Eigen::MatrixXd EV;
+  Eigen::VectorXi EI;
+  igl::fast_find_intersections(V,F,Vp,Fp,false,false,IF,EV,EE,EI);
+  {
+    Eigen::VectorXi I;
+    igl::unique(IF,I);
+    // should have 8 triangles that are intersecting plane (so 9 unique)
+    REQUIRE( I.size()== 9 );
+  }
 
   // all triangles from the cube intersect the same plane
-  REQUIRE( (I.col(1).array()==0).all());
+  REQUIRE( (IF.col(1).array()==0).all() );
   
-  // and 16 line edges
-  REQUIRE ( edges.rows()==16);
+  // and 8 line edges
+  REQUIRE( EE.rows()==8 );
+  REQUIRE( EV.rows()==2*8 );
+  REQUIRE( EI.size()==8 );
   // TODO: check if the edges are at the right place
 }
+
+TEST_CASE("fast_find_intersections: coplanar", "[igl]")
+{
+}

tests/include/igl/fast_find_self_intersections.cpp

@@ -0,0 +1,129 @@
+#include "test_common.h"
+#include <igl/fast_find_self_intersections.h>
+#include <igl/combine.h>
+#include <igl/sortrows.h>
+#include <igl/unique.h>
+#include <igl/matlab_format.h>
+#include <igl/PI.h>
+#include <iostream>
+#include <iomanip>
+
+
+TEST_CASE("fast_find_self_intersections: cube", "[igl]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+
+  Eigen::MatrixXi IF,EE;
+  Eigen::MatrixXd EV;
+  Eigen::VectorXi EI;
+  REQUIRE( !igl::fast_find_self_intersections(V,F,true,true,IF,EV,EE,EI) );
+
+  REQUIRE(IF.rows()==0);
+  REQUIRE(EV.rows()==0);
+  REQUIRE(EE.rows()==0);
+  REQUIRE(EI.size()==0);
+}
+
+TEST_CASE("fast_find_self_intersections: cube-triangle", "[igl]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+
+  auto dims=(V.colwise().maxCoeff()-V.colwise().minCoeff())*2;
+  auto cz=(V.col(2).minCoeff()+V.col(2).maxCoeff())/2;
+
+  // add a triangle intersecting cube 
+  Eigen::MatrixXd Vp(3,3);
+  Vp << V.col(0).minCoeff()-dims(0),   V.col(1).minCoeff()-dims(1),   cz,
+        V.col(0).minCoeff()-dims(0),   V.col(1).maxCoeff()+dims(1)*3, cz,
+        V.col(0).maxCoeff()+dims(0)*3, V.col(1).maxCoeff()-dims(1),   cz;
+  Eigen::MatrixXi Fp(1,3);
+  Fp << 0,1,2;
+  
+  Eigen::MatrixXd V_;
+  Eigen::MatrixXi F_;
+  igl::combine<Eigen::MatrixXd,Eigen::MatrixXi>({V,Vp},{F,Fp}, V_,F_);
+
+
+  Eigen::MatrixXi IF,EE;
+  Eigen::MatrixXd EV;
+  Eigen::VectorXi EI;
+  REQUIRE( igl::fast_find_self_intersections(V_,F_,false,false,IF,EV,EE,EI) );
+  {
+    Eigen::VectorXi I;
+    igl::unique(IF,I);
+    // should have 9 triangles that are intersecting each other
+    REQUIRE( I.size()==9 );
+    // plane that we added intersects other triangles
+    REQUIRE((IF.col(1).array()==(F_.rows()-1)).all() );
+  }
+  // and 8 edges
+  REQUIRE( EE.rows()==8 );
+  REQUIRE( EV.rows()==2*8 );
+  REQUIRE( EI.size()==8 );
+}
+
+
+TEST_CASE("fast_find_self_intersections: rose", "[igl]")
+{
+  Eigen::MatrixXd V(10,3);
+  Eigen::MatrixXi F(9,3);
+  for(int i=0;i<9;i++)
+  {
+    const double theta_i = 4.0*igl::PI*double(i)/9.0;
+    V.row(i) << std::cos(theta_i), std::sin(theta_i), 1;
+    F.row(i)<<9,i,(i+1)%9;
+  }
+  V.row(9) << 0,0,0;
+
+  Eigen::MatrixXi IF,EE;
+  Eigen::MatrixXd EV;
+  Eigen::VectorXi EI;
+  REQUIRE( igl::fast_find_self_intersections(V,F,false,false,IF,EV,EE,EI) );
+  Eigen::MatrixXi IF_gt(9,2);
+  IF_gt<<
+    0,4,
+    0,5,
+    1,5,
+    1,6,
+    2,6,
+    2,7,
+    3,7,
+    3,8,
+    4,8;
+  {
+    Eigen::MatrixXi sIF;
+    Eigen::VectorXi _;
+    igl::sortrows(Eigen::MatrixXi(IF),true,sIF,_);
+    test_common::assert_eq(IF_gt,sIF);
+  }
+}
+
+TEST_CASE("fast_find_self_intersections: shared-edge", "[igl]")
+{
+  Eigen::MatrixXd V(4,3);
+  V<<
+    0,0,0,
+    1,0,0,
+    0,1,0,
+    -1,0,0;
+  Eigen::MatrixXi F(2,3);
+  F <<
+    0,1,2,
+    0,2,3;
+  Eigen::MatrixXi IF,EE;
+  Eigen::MatrixXd EV;
+  Eigen::VectorXi EI;
+  REQUIRE(!igl::fast_find_self_intersections(V,F,false,false,IF,EV,EE,EI));
+  V.row(3) << 2.0,0,0;
+  REQUIRE( igl::fast_find_self_intersections(V,F,false,false,IF,EV,EE,EI));
+  REQUIRE( IF.rows()==1 );
+  REQUIRE( EE.rows()==0 );
+  REQUIRE( EV.rows()==0 );
+  REQUIRE( EI.rows()==0 );
+}

tests/include/igl/intersection_blocking_collapse_edge_callbacks.cpp

@@ -0,0 +1,143 @@
+#include <test_common.h>
+#include <igl/intersection_blocking_collapse_edge_callbacks.h>
+#include <igl/edge_flaps.h>
+#include <igl/decimate_callback_types.h>
+#include <igl/AABB.h>
+#include <igl/collapse_edge.h>
+#include <igl/min_heap.h>
+#include <tuple>
+
+TEST_CASE("intersection_blocking_collapse_edge_callbacks: simple", "[igl]")
+{
+  // A mesh where collapsing edge e1 will create a self-intersection and
+  // collapsed e2 will not
+  Eigen::MatrixXd V(21,3);
+  V<< 
+    0,20,0,
+    0,30,0,
+    10,10,0,
+    10,20,-30,
+    10,30,0,
+    20,0,0,
+    20,10,30,
+    20,20,0,
+    30,0,0,
+    30,10,0,
+    10,17,-15,
+    7,23,15,
+    13,23,15,
+    10,30,30,
+    20,30,0,
+    0,10,3,
+    20,5,15,
+    25,0,0,
+    25,5,15,
+    25,10,15,
+    30,5,0;
+  Eigen::MatrixXi F(22,3);
+  F<< 
+    5,17,16,
+    8,20,18,
+    8,18,17,
+    9,19,20,
+    6,16,18,
+    6,18,19,
+    18,16,17,
+    19,18,20,
+    2,5,16,
+    7,6,19,
+    2,16,6,
+    7,19,9,
+    2,3,0,
+    0,3,1,
+    3,4,1,
+    2,6,3,
+    6,7,3,
+    3,7,4,
+    10,11,12,
+    12,11,13,
+    4,7,14,
+    2,0,15;
+  igl::AABB<Eigen::MatrixXd, 3> * tree = new igl::AABB<Eigen::MatrixXd, 3>();
+  tree->init(V,F);
+
+
+  igl::decimate_pre_collapse_callback  pre_collapse;
+  igl::decimate_post_collapse_callback post_collapse;
+  igl::intersection_blocking_collapse_edge_callbacks(
+    tree,
+    pre_collapse,
+    post_collapse);
+  
+  Eigen::VectorXi EMAP;
+  Eigen::MatrixXi E,EF,EI;
+  igl::edge_flaps(F,E,EMAP,EF,EI);
+  igl::min_heap< std::tuple<double,int,int> > Q;
+  Eigen::VectorXi EQ;
+  Eigen::MatrixXd C = Eigen::MatrixXd::Zero(E.rows(),3);
+  // Try to do the collapses.
+  std::vector<int> edges_to_collapse;
+  for(const auto edge : {Eigen::RowVector2i(3,6),Eigen::RowVector2i(6,18)})
+  {
+    int e = -1;
+    const bool found = 
+      ( (E.array().col(0)==edge[0] && E.array().col(1)==edge[1])||
+        (E.array().col(0)==edge[1] && E.array().col(1)==edge[0])).maxCoeff(&e);
+    REQUIRE(found);
+    edges_to_collapse.push_back(e);
+  }
+  std::vector<bool> should_collapse = {false,true};
+  for(int i=0;i<edges_to_collapse.size();i++)
+  {
+    const auto e = edges_to_collapse[i];
+    bool collapsed = true;
+    int e1,e2,f1,f2;
+    // pre_collapse only has access to C not p
+    C.row(e) = 0.5*(V.row(E(e,0))+V.row(E(e,1)));
+    if(pre_collapse(V,F,E,EMAP,EF,EI,Q,EQ,C,e))
+    {
+      collapsed = igl::collapse_edge(
+        e,C.row(e).eval(),
+        /*Nsv,Nsf,Ndv,Ndf,*/
+        V,F,E,EMAP,EF,EI,e1,e2,f1,f2);
+    }else
+    {
+      // Aborted by pre collapse callback
+      collapsed = false;
+    }
+    REQUIRE(collapsed == should_collapse[i]);
+    post_collapse(V,F,E,EMAP,EF,EI,Q,EQ,C,e,e1,e2,f1,f2,collapsed);
+  }
+
+  {
+    auto * root = tree->root();
+    const auto leaves = tree->gather_leaves();
+    for(int f = 0;f<F.rows();f++)
+    {
+      if((F.row(f).array() == IGL_COLLAPSE_EDGE_NULL).all())
+      {
+        continue;
+      }
+      REQUIRE(f < leaves.size());
+      auto * leaf = leaves[f];
+      REQUIRE(root == leaf->root());
+      // check containment all the way up to
+      auto * node = leaf;
+      while(node)
+      {
+        REQUIRE(node->m_box.contains(V.row(F(f,0)).transpose()));
+        REQUIRE(node->m_box.contains(V.row(F(f,1)).transpose()));
+        REQUIRE(node->m_box.contains(V.row(F(f,2)).transpose()));
+        auto * parent = node->m_parent;
+        if(parent)
+        {
+          REQUIRE(parent->m_box.contains(node->m_box));
+        }
+        node = parent;
+      }
+    }
+  }
+
+  tree = tree->root();
+  delete tree;
+}

tests/include/igl/qslim.cpp

@@ -17,7 +17,7 @@ TEST_CASE("qslim: cylinder", "[igl]" "[slow]")
   Eigen::MatrixXd U;
   Eigen::MatrixXi G;
   Eigen::VectorXi I,J;
-  qslim(V,F,2*axis_devisions,U,G,I,J);
+  qslim(V,F,2*axis_devisions,false,U,G,I,J);
   REQUIRE (U.rows() == axis_devisions*2);
   //double l,u;
   igl::writePLY("qslim-cylinder-vf.ply",V,F);

tests/include/igl/tri_tri_intersect.cpp

@@ -0,0 +1,57 @@
+#include "test_common.h"
+#include <igl/tri_tri_intersect.h>
+
+TEST_CASE("tri_tri_intersection_test_3d intersect", "[igl]")
+{
+  // try with oblique interecion plane
+  for(double shift=-1000;shift<=1000;shift+=100.0)
+  {
+    Eigen::RowVector3d p1(0,0,0),    q1(1,0,0),r1(0,1,0);
+    Eigen::RowVector3d p2(shift,0,1),q2(1,1,0),r2(-shift,0,-1);
+
+    // should intersect along the vector (0,0,0) -> (0.5,0.5,0)
+    Eigen::RowVector3d s,t;
+    bool coplanar;
+    REQUIRE( igl::tri_tri_intersection_test_3d(p1,q1,r1, p2,q2,r2, coplanar, s, t) );
+    REQUIRE( !coplanar);
+    
+    if(s.norm()<1e-5)
+    {
+      Eigen::RowVector3d t_correct(0.5,0.5,0);
+      Eigen::RowVector3d s_correct(0,0,0);
+      test_common::assert_near( t, t_correct, 1e-10);
+      test_common::assert_near( s, s_correct, 1e-10);
+    } else {
+      Eigen::RowVector3d s_correct(0.5,0.5,0);
+      Eigen::RowVector3d t_correct(0,0,0);
+      test_common::assert_near( t, t_correct, 1e-10);
+      test_common::assert_near( s, s_correct, 1e-10);
+    }
+  }
+}
+
+TEST_CASE("tri_tri_intersection_test_3d not intersect", "[igl]")
+{
+  // positive test that two triangles intersect
+  Eigen::RowVector3d p1(0,0,0),q1(1,0,0),r1(0,1,0);
+  Eigen::RowVector3d p2(0,0,1),q2(1,0,1),r2(0,1,1);
+
+  // should intersect along the vector (0,0,0) -> (sqrt(2),sqrt(2),0)
+  Eigen::RowVector3d s,t;
+  bool coplanar;
+  REQUIRE( !igl::tri_tri_intersection_test_3d(p1,q1,r1,p2,q2,r2, coplanar, s, t) );
+}
+
+
+TEST_CASE("tri_tri_intersection_test_3d coplanar", "[igl]")
+{
+  // positive test that two triangles intersect
+  Eigen::RowVector3d p1(0,0,0),q1(1,0,0),r1(0,1,0);
+  Eigen::RowVector3d p2(0,0,0),q2(0.5,0,0),r2(0,0.5,0);
+
+  // should intersect along the vector (0,0,0) -> (sqrt(2),sqrt(2),0)
+  Eigen::RowVector3d s,t;
+  bool coplanar;
+  REQUIRE( igl::tri_tri_intersection_test_3d(p1,q1,r1,p2,q2,r2, coplanar, s, t) );
+  REQUIRE(coplanar);
+}

tests/include/igl/triangle_triangle_intersect.cpp

@@ -0,0 +1,47 @@
+#include "test_common.h"
+#include <igl/triangle_triangle_intersect.h>
+#include <igl/edge_flaps.h>
+
+TEST_CASE("triangle_triangle_intersect: shared-edge", "[igl]" )
+{
+  Eigen::MatrixXd V(4,3);
+  V<<
+    0,0,0,
+    1,0,0,
+    0,1,0,
+    -1,0,0;
+  Eigen::MatrixXi F(2,3);
+  F <<
+    0,1,2,
+    0,2,3;
+  //      2
+  //     /|\
+  //    / | \
+  //   /  |  \
+  //  3---0---1
+  const int f = 0;
+  const int c = 1;
+  const int g = 1;
+
+  Eigen::MatrixXi E,EF,EI;
+  Eigen::VectorXi EMAP;
+  igl::edge_flaps(F,E,EMAP,EF,EI);
+
+
+  bool ret;
+  ret = igl::triangle_triangle_intersect(V,F,E,EMAP,EF,f,c,V.row(F(f,c)),g);
+  REQUIRE(ret == false);
+
+  for(const double epsilon : {0.,1e-15,-1e-15})
+  {
+    //  2
+    //  |\⟍
+    //  | \ ⟍
+    //  |  \  ⟍
+    //  0---1---3
+    V.row(3) << 2.0+epsilon,0,0;
+    ret = igl::triangle_triangle_intersect(V,F,E,EMAP,EF,f,c,V.row(F(f,c)),g);
+    REQUIRE(ret == true);
+  }
+
+}

tutorial/903_FastFindSelfIntersections/main.cpp

@@ -1,58 +1,38 @@
-#include <igl/readOFF.h>
-#include <igl/combine.h>
-
+#include <igl/read_triangle_mesh.h>
+#include <igl/unique.h>
 #include <igl/opengl/glfw/Viewer.h>
-//#include <igl/opengl/glfw/imgui/ImGuiPlugin.h>
-#include <igl/opengl/glfw/imgui/ImGuiMenu.h>
-
 #include <igl/fast_find_self_intersections.h>
 
-Eigen::MatrixXd V1,V2,V;
-Eigen::MatrixXi F1,F2,F;
-
 int main(int argc, char *argv[])
 {
 
-  // Load two meshes 
-  igl::readOFF(TUTORIAL_SHARED_PATH "/planexy.off", V1, F1);
-  igl::readOFF(TUTORIAL_SHARED_PATH "/cow.off",     V2, F2);
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(argc<=1?TUTORIAL_SHARED_PATH "/cow.off":argv[1], V, F);
 
-  // Combine into one mesh (will produce self-intersections)
-  igl::combine<Eigen::MatrixXd,Eigen::MatrixXi>({V1,V2},{F1,F2}, V,F);
 
   // Plot the mesh
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V, F);
 
-  Eigen::VectorXi I;
-  Eigen::MatrixXd edges;
+  Eigen::VectorXi EI;
+  Eigen::MatrixXd EV;
+  Eigen::MatrixXi IF,EE;
 
-  if(igl::fast_find_self_intersections(V,F,I,edges))
+  if(igl::fast_find_self_intersections(V,F,false,false,IF,EV,EE,EI))
   {
-    std::cout<<"Found "<<I.sum()<<" self intersections"<<std::endl;
-
+    std::cout<<"Found "<<IF.rows()<<" self intersecting pairs"<<std::endl;
     // plot edge vertices
-    //viewer.data().add_points(edges, Eigen::RowVector3d(1,0,0));
-
-    // Plot the edges of the self intersects
-    for (unsigned i=0;i<edges.rows(); i+=2)
-    {
-      viewer.data().add_edges
-      (
-        edges.row(i),
-        edges.row(i+1),
-        Eigen::RowVector3d(1,0,0)
-      );
-    }
-    std::cout<<std::endl;
+    viewer.data().set_edges(EV,EE, Eigen::RowVector3d(1,0,0));
   }
-  viewer.data().set_data(I.cast<double>());
+  Eigen::VectorXi I;
+  igl::unique(IF,I);
+  Eigen::VectorXd D = Eigen::MatrixXd::Zero(F.rows(),1);
+  D(I).setConstant(1.0);
+  viewer.data().set_data(D,0,1,igl::COLOR_MAP_TYPE_PARULA);
+  viewer.data().set_face_based(true);
   viewer.data().double_sided=true;
 
-  igl::opengl::glfw::imgui::ImGuiMenu menu;
-  //plugin.widgets.push_back(&menu);
-  menu.callback_draw_viewer_window = [](){};
-
   // Launch the viewer
   viewer.launch();
 }

tutorial/904_FastFindIntersections/main.cpp

@@ -1,8 +1,8 @@
-#include <igl/readOFF.h>
-#include <igl/combine.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/AABB.h>
+#include <igl/unique.h>
 
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/opengl/glfw/imgui/ImGuiMenu.h>
 
 #include <igl/fast_find_intersections.h>
 
@@ -17,26 +17,24 @@ double slice_z;
 
 void update_visualization(igl::opengl::glfw::Viewer & viewer)
 {
-  Eigen::MatrixXi I;
-  Eigen::MatrixXd edges;
-
   //shifted intersection object
   Eigen::MatrixXd V2_(V2.rows(),V2.cols());
   V2_<< V2.col(0), V2.col(1), V2.col(2).array()+slice_z;
 
-  igl::fast_find_intersections(tree, V1,F1, V2_,F2, I,edges);
-  Eigen::MatrixXi edges_link=Eigen::MatrixXi::NullaryExpr(edges.rows()/2,2, [](int i,int j) { return i*2+j;});
+  Eigen::MatrixXi IF,EE;
+  Eigen::MatrixXd EV;
+  Eigen::VectorXi EI;
+  igl::fast_find_intersections(tree, V1,F1, V2_,F2,false,false,IF,EV,EE,EI);
  
   // Plot the edges of the intersects
-  viewer.data().set_edges ( edges, edges_link, Eigen::RowVector3d(1,0,0));
+  viewer.data().set_edges( EV,EE, Eigen::RowVector3d(1,0,0));
   
   // show faces which are intersected
-  Eigen::VectorXd face_data=Eigen::VectorXd::Zero(F1.rows());
-
-  for(int i=0; i<I.rows(); ++i)
-    face_data(I(i,0)) = 1.0;
-  
-  viewer.data().set_data(face_data);
+  Eigen::VectorXi I;
+  igl::unique(IF,I);
+  Eigen::VectorXd D = Eigen::MatrixXd::Zero(F1.rows(),1);
+  D(I).setConstant(1.0);
+  viewer.data().set_data(D,0,1,igl::COLOR_MAP_TYPE_PARULA);
 }
 
 
@@ -61,8 +59,10 @@ bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int mod)
 int main(int argc, char *argv[])
 {
   // Load two meshes 
-  igl::readOFF(TUTORIAL_SHARED_PATH "/cow.off",     V1, F1);
-  igl::readOFF(TUTORIAL_SHARED_PATH "/planexy.off", V2, F2);
+  igl::read_triangle_mesh(
+    argc<=1?TUTORIAL_SHARED_PATH "/cow.off"    :argv[1],V1,F1);
+  igl::read_triangle_mesh(
+    argc<=2?TUTORIAL_SHARED_PATH "/planexy.off":argv[2],V2,F2);
 
   // initialize AABB tree with first mesh (it doesn't move)
   tree.init(V1,F1);
@@ -81,13 +81,11 @@ int main(int argc, char *argv[])
   // Plot the mesh
   igl::opengl::glfw::Viewer viewer;
   viewer.data().set_mesh(V1, F1);
+  viewer.data().set_face_based(true);
   
 
   update_visualization(viewer);
 
-  igl::opengl::glfw::imgui::ImGuiMenu menu;
-  //plugin.widgets.push_back(&menu);
-  menu.callback_draw_viewer_window = [](){};
   viewer.callback_key_down = &key_down;
 
   // show the cow closer

tutorial/907_DynamicAABB/main.cpp

@@ -0,0 +1,153 @@
+#include <igl/AABB.h>
+#include <igl/box_faces.h>
+#include <igl/colon.h>
+#include <igl/find.h>
+#include <igl/get_seconds.h>
+#include <igl/per_face_normals.h>
+#include <igl/quad_edges.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/unproject_ray.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <limits>
+#include <deque>
+#include <stdio.h>
+
+
+
+int main(int argc, char *argv[])
+{
+  Eigen::MatrixXd V,N;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(
+    argc>1?argv[1]:TUTORIAL_SHARED_PATH "/decimated-knight.off",V,F);
+  // Make mesh into disconnected soup
+  V = V(Eigen::Map<Eigen::VectorXi>(F.data(),F.size()), Eigen::all).eval();
+  F = Eigen::Map<Eigen::MatrixXi>(igl::colon<int>(0,V.rows()-1).data(),V.rows()/3,3).eval();
+  // Cache normals
+  igl::per_face_normals(V,F,N);
+
+  // Build an initial tree. Store as pointer so we can update root.
+  auto * tree = new igl::AABB<Eigen::MatrixXd, 3>();
+  tree->init(V,F);
+  // Gather leaf pointers (these are stable)
+  const auto leaves = tree->gather_leaves(F.rows());
+
+  // Set up visualization and interaction
+  igl::opengl::glfw::Viewer vr;
+  Eigen::MatrixXd C = Eigen::RowVector3d(0.9,0.9,0.9).replicate(F.rows(),1);
+  vr.data().set_mesh(V,F);
+  vr.data().set_face_based(true);
+  vr.data().set_colors(C);
+
+  // draw the boxes of the tree at a given depth
+  int depth = 0;
+  Eigen::MatrixXd TV;
+  Eigen::MatrixXi TQ;
+  Eigen::VectorXi TD;
+  const auto update_edges = [&]()
+  {
+    Eigen::MatrixXi TQd = TQ(igl::find((TD.array()==depth).eval()),Eigen::all);
+    Eigen::MatrixXi TE;
+    igl::quad_edges(TQd,TE);
+    vr.data().set_edges(TV,TE,Eigen::RowVector3d(1,1,1));
+    vr.data().line_width = 2;
+  };
+  const auto update_tree_vis = [&]()
+  {
+    igl::box_faces(*tree,TV,TQ,TD);
+    update_edges();
+  };
+
+  // Update the tree for each triangle in hits which might have moved.
+  const auto update_tree = [&](const std::vector<igl::Hit> & hits)
+  {
+    for(const auto hit : hits)
+    {
+      const int fid = hit.id;
+      Eigen::AlignedBox<double,3> box;
+      box.extend(V.row(F(fid,0)).transpose());
+      box.extend(V.row(F(fid,1)).transpose());
+      box.extend(V.row(F(fid,2)).transpose());
+      // This is O(height) which is typically O(log n)
+      tree = leaves[fid]->update(box,0)->root();
+    }
+    // update the visualization. This is O(n) 🙃
+    update_tree_vis();
+  };
+
+  update_tree_vis();
+
+  std::vector<igl::Hit> hits;
+  Eigen::RowVector3d dir;
+  Eigen::RowVector3d red(1,0.2,0.2);
+  vr.callback_pre_draw = [&](decltype(vr) &)->bool
+  {
+    // While mouse is down, pull front (back) faces towards (away from) cursor.
+    if(hits.size())
+    {
+      for(const auto hit : hits)
+      {
+        for(int c : {0,1,2})
+        {
+          const int fid = hit.id;
+          V.row(F(fid,c)) += 
+            ((N.row(fid).dot(dir))>0?1:-1) *
+            0.001*dir.cast<double>().transpose();
+          C.row(fid) = red;
+        }
+      }
+      // Things have moved, so update tree ~O(#hits log n)
+      update_tree(hits);
+      // update viewer. This is O(n) 🙃
+      vr.data().set_vertices(V);
+      vr.data().set_colors(C);
+    }
+    return false;
+  };
+  vr.callback_mouse_down = [&](decltype(vr) &, int button, int modifier)->bool
+  {
+    // Shoot a ray through cursor into mesh accelerated by current AABB tree
+    const Eigen::Vector2f pos(
+      vr.current_mouse_x,vr.core().viewport(3)-vr.current_mouse_y);
+    const auto model = vr.core().view;
+    const auto proj = vr.core().proj;
+    const auto viewport = vr.core().viewport;
+    Eigen::RowVector3d src;
+    {
+      Eigen::Vector3f src_f,dir_f;
+      igl::unproject_ray(pos,model,proj,viewport,src_f,dir_f);
+      src = src_f.cast<double>().transpose();
+      dir = dir_f.cast<double>().transpose();
+    }
+    if(tree->intersect_ray(V,F,src,dir,hits))
+    {
+      vr.core().is_animating = true;
+      return true;
+    }
+    return false;
+  };
+  vr.callback_mouse_up = [&](decltype(vr) &, int button, int modifier)->bool
+  {
+    hits.clear();
+    vr.core().is_animating = false;
+    return false;
+  };
+  vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
+  {
+    switch(key)
+    {
+      case ',':
+      case '.':
+        depth = std::max(depth+(key==','?-1:1),0);
+        update_edges();
+        return true;
+      default:
+        return false;
+    }
+  };
+  printf("  ,/. to change tree-vis depth\n");
+  printf("  Click and hold to push triangles\n");
+  vr.launch();
+  // delete the tree (and all subtrees/leaves)
+  delete tree;
+}

tutorial/908_IntersectionBlockingDecimation/main.cpp

@@ -0,0 +1,168 @@
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/fast_find_self_intersections.h>
+#include <igl/unique.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/get_seconds.h>
+#include <igl/AABB.h>
+#include <igl/writePLY.h>
+#include <igl/intersection_blocking_collapse_edge_callbacks.h>
+#include <igl/qslim_optimal_collapse_edge_callbacks.h>
+#include <igl/per_vertex_point_to_plane_quadrics.h>
+#include <igl/STR.h>
+#include <igl/connect_boundary_to_infinity.h>
+#include <igl/decimate.h>
+#include <igl/max_faces_stopping_condition.h>
+#include <igl/point_simplex_squared_distance.h>
+#include <igl/edge_flaps.h>
+#include <igl/decimate_callback_types.h>
+#include <igl/find.h>
+
+int main(int argc, char *argv[])
+{
+  IGL_TICTOC_LAMBDA;
+  Eigen::MatrixXd V,V0;
+  Eigen::MatrixXi F,F0;
+  igl::read_triangle_mesh( 
+      argc<=1 ? TUTORIAL_SHARED_PATH "/octopus-low.mesh" :
+      argv[1],V,F);
+  V0 = V;F0 = F;
+  ///////////////////////////////////////////////////////////
+  /// Before collapsing starts
+  ///////////////////////////////////////////////////////////
+  tictoc();
+  igl::AABB<Eigen::MatrixXd, 3> * tree = new igl::AABB<Eigen::MatrixXd, 3>();
+  tree->init(V,F);
+  printf("tree->init: %g\n",tictoc());
+  tree->validate();
+
+  const int target_m = F.rows() * 0.1 + 1;
+  Eigen::MatrixXd dV[2];
+  Eigen::MatrixXd dC[2];
+  Eigen::MatrixXi dF[2];
+  Eigen::RowVector3d gray(0.9,0.9,0.9);
+  for(auto pass : {0,1})
+  {
+    Eigen::MatrixXd VO;
+    Eigen::MatrixXi FO;
+    igl::connect_boundary_to_infinity(V,F,VO,FO);
+    Eigen::VectorXi EMAP;
+    Eigen::MatrixXi E,EF,EI;
+    igl::edge_flaps(FO,E,EMAP,EF,EI);
+
+    igl::decimate_cost_and_placement_callback cost_and_placement;
+    igl::decimate_pre_collapse_callback  pre_collapse;
+    igl::decimate_post_collapse_callback post_collapse;
+
+    // Quadrics per vertex
+    typedef std::tuple<Eigen::MatrixXd,Eigen::RowVectorXd,double> Quadric;
+    std::vector<Quadric> quadrics;
+    igl::per_vertex_point_to_plane_quadrics(VO,FO,EMAP,EF,EI,quadrics);
+    // State variables keeping track of edge we just collapsed
+    int v1 = -1;
+    int v2 = -1;
+    // Callbacks for computing and updating metric
+    igl::qslim_optimal_collapse_edge_callbacks(
+      E,quadrics,v1,v2, cost_and_placement, pre_collapse,post_collapse);
+
+    if(pass == 1)
+    {
+      igl::intersection_blocking_collapse_edge_callbacks(
+        pre_collapse, post_collapse, // These will get copied as needed
+        tree,
+        pre_collapse, post_collapse);
+    }
+
+    int m = F.rows();
+    const int orig_m = m;
+    Eigen::MatrixXd U;
+    Eigen::MatrixXi G;
+    Eigen::VectorXi J,I;
+    tictoc();
+    const bool ret = igl::decimate(
+      VO, FO,
+      cost_and_placement,
+      igl::max_faces_stopping_condition(m,orig_m,target_m),
+      pre_collapse,
+      post_collapse,
+      E, EMAP, EF, EI,
+      U, G, J, I);
+    G = G(igl::find((J.array()<orig_m).eval()), Eigen::all).eval();
+    {
+      Eigen::VectorXi _;
+      igl::remove_unreferenced(Eigen::MatrixXd(U),Eigen::MatrixXi(G),U,G,_);
+    }
+    printf("qslim-%22s in %g secs\n",pass?"-intersection-blocking":"",tictoc());
+    igl::writePLY(STR("out-"<<pass<<".ply"),U,G);
+    dV[pass] = U;
+    dF[pass] = G;
+    {
+      Eigen::VectorXi BI;
+      {
+        Eigen::MatrixXd EV;
+        Eigen::MatrixXi IF,EE;
+        Eigen::VectorXi EI;
+        igl::fast_find_self_intersections(dV[pass],dF[pass],true,false,IF,EV,EE,EI);
+        igl::unique(IF,BI);
+      }
+      printf("  # self-intersections: %d\n",(int)BI.size());
+      dC[pass] = gray.replicate(dF[pass].rows(),1);
+      dC[pass](BI,Eigen::all) = 
+        Eigen::RowVector3d(0.95,0.15,0.15).replicate(BI.size(),1);
+    }
+  }
+  
+
+  tree->validate();
+  assert(tree == tree->root());
+  tree = tree->root();
+
+  igl::opengl::glfw::Viewer vr;
+  vr.callback_key_pressed = [&](decltype(vr) &,unsigned int key, int mod)
+  {
+    switch(key)
+    {
+      case ',':
+      case '.':
+        vr.data_list[vr.selected_data_index].is_visible = false;
+        vr.selected_data_index += (key==','?-1:1) + vr.data_list.size();
+        vr.selected_data_index %= vr.data_list.size();
+        vr.data_list[vr.selected_data_index].is_visible = true;
+        switch(vr.selected_data_index)
+        {
+          case 0:
+            printf("Original mesh\n");
+            break;
+          case 1:
+            printf("Qslim mesh\n");
+            break;
+          case 2:
+            printf("Qslim mesh with intersection blocking\n");
+            break;
+        }
+        return true;
+      default:
+        return false;
+    }
+  };
+
+  vr.data().set_mesh(V0,F0);
+  vr.data().show_lines = false;
+  vr.data().set_face_based(true);
+  vr.data().is_visible = false;
+  vr.data().set_colors(gray);
+  vr.append_mesh();
+  vr.data().set_mesh(dV[0],dF[0]);
+  vr.data().show_lines = true;
+  vr.data().set_face_based(true);
+  vr.data().is_visible = false;
+  vr.data().set_colors(dC[0]);
+  vr.append_mesh();
+  vr.data().set_mesh(dV[1],dF[1]);
+  vr.data().show_lines = true;
+  vr.data().set_face_based(true);
+  vr.data().set_colors(dC[1]);
+  vr.launch();
+
+
+}

tutorial/CMakeLists.txt

@@ -125,8 +125,10 @@ endif()
 if(LIBIGL_TUTORIALS_CHAPTER9)
     igl_add_tutorial(901_VectorFieldSmoothing igl::glfw)
     igl_add_tutorial(902_VectorParallelTransport igl::glfw)
-    igl_add_tutorial(903_FastFindSelfIntersections igl::imgui igl::glfw)
-    igl_add_tutorial(904_FastFindIntersections igl::imgui igl::glfw)
+    igl_add_tutorial(903_FastFindSelfIntersections igl::glfw)
+    igl_add_tutorial(904_FastFindIntersections igl::glfw)
     igl_add_tutorial(905_Isolines igl::imgui igl::glfw)
     igl_add_tutorial(906_TrimWithSolid igl::glfw igl_copyleft::cgal)
+    igl_add_tutorial(907_DynamicAABB igl::glfw)
+    igl_add_tutorial(908_IntersectionBlockingDecimation igl::glfw)
 endif()