Add and fix test for fast_find_self_intersections (#2382)

* add test case

* add coplanarity test

* another failing test

* use orient3d

* no printing

* reverting... That introduced lots of other failure cases

* subdivided knight case

* wip predicates

* promising predicates version

* tri_tri_overlap compiles, header guards, predicate find_*

* tests for predicates::find_in...

* parallel for

* note

* mv to predicates

* remove old functions

* tutorials; extracting segments is broken 904

* fix extraction bug

* missing header

* capture consts

include/igl/fast_find_intersections.cpp

@@ -1,253 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2022 Vladimir S. FONOV <[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_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 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>;
-  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++;
-    }
-
-    if(num_if >= IF.rows())
-    {
-      IF.conservativeResize(2*IF.rows()+1,2);
-    }
-    IF.row(num_if++) << f1,f2;
-  };
-
-  // 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++)
-      {
-        // Find in opposite direction on jth face
-        if(F1(g,e) == d && F1(g,(e+1)%3) == s)
-        {
-          return c;
-        }
-      }
-    }
-    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++)
-    {
-      for(sg = 0;sg<3;sg++)
-      {
-        if(F1(g,sg) == F1(f,sf))
-        {
-          return true;
-        }
-      }
-    }
-    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)
-    {
-      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
-// 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,84 +0,0 @@
-// 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/
-#pragma once
-#ifndef FAST_FIND_MESH_INTERSECT_H
-#define FAST_FIND_MESH_INTERSECT_H
-
-#include "igl_inline.h"
-#include <Eigen/Core>
-
-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  #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 <
-    typename DerivedV1,
-    typename DerivedF1,
-    typename DerivedV2,
-    typename DerivedF2,
-    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
-  /// \brief Tree built internally.
-  template <
-    typename DerivedV1,
-    typename DerivedF1,
-    typename DerivedV2,
-    typename DerivedF2,
-    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 "fast_find_intersections.cpp"
-#endif
-
-#endif

include/igl/fast_find_self_intersections.cpp

@@ -1,39 +0,0 @@
-// 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 "fast_find_intersections.h"
-
-template <
-  typename DerivedV,
-  typename DerivedF,
-  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,
-  const bool detect_only,
-  const bool first_only,
-  Eigen::PlainObjectBase<DerivedIF> & IF,
-  Eigen::PlainObjectBase<DerivedEV> & EV,
-  Eigen::PlainObjectBase<DerivedEE> & EE,
-  Eigen::PlainObjectBase<DerivedEI> & EI)
-{
-  // 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);
-}
-
-#ifdef IGL_STATIC_LIBRARY
-// Explicit template instantiation
-// 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

@@ -1,53 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2022 Vladimir S. FONOV <[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/
-#pragma once
-#ifndef FAST_FIND_SELF_INTERSECTIONS_H
-#define FAST_FIND_SELF_INTERSECTIONS_H
-
-#include "igl_inline.h"
-#include <Eigen/Core>
-
-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 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::SelfIntersectMesh
-  template <
-    typename DerivedV,
-    typename DerivedF,
-    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,
-    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 "fast_find_self_intersections.cpp"
-#endif
-
-#endif

include/igl/predicates/find_intersections.cpp

@@ -0,0 +1,293 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2022 Vladimir S. FONOV <[email protected]>
+// Copyright (C) 2024 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 "find_intersections.h"
+
+#include "../AABB.h"
+#include "../triangle_triangle_intersect_shared_edge.h"
+#include "../triangle_triangle_intersect_shared_vertex.h"
+#include "../find.h"
+#include "../list_to_matrix.h"
+#include "triangle_triangle_intersect.h"
+#include "../triangle_triangle_intersect.h"
+#include <stdio.h>
+// atomic
+#include <igl/parallel_for.h>
+#include <atomic>
+#include <mutex>
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedIF,
+  typename DerivedCP >
+IGL_INLINE bool igl::predicates::find_intersections(
+  const igl::AABB<DerivedV1,3> & tree1,
+  const Eigen::MatrixBase<DerivedV1> & V1,
+  const Eigen::MatrixBase<DerivedF1> & F1,
+  const Eigen::MatrixBase<DerivedV2> & V2,
+  const Eigen::MatrixBase<DerivedF2> & F2,
+  const bool first_only,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedCP> & CP)
+{
+  const bool detect_only = true;
+  constexpr bool stinker = false;
+  using AABBTree=igl::AABB<DerivedV1,3>;
+  using Scalar=typename DerivedV1::Scalar;
+  static_assert(
+    std::is_same<Scalar,typename DerivedV2::Scalar>::value,
+    "V1 and V2 must have the same scalar type");
+  using RowVector3S=Eigen::Matrix<Scalar,1,3>;
+
+  // 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;
+  // mutex
+  std::mutex append_mutex;
+  const auto append_intersection = 
+    [&IF,&CP,&num_if,&append_mutex]( const int f1, const int f2, const bool coplanar = false)
+  {
+    std::lock_guard<std::mutex> lock(append_mutex);
+    if(num_if >= IF.rows())
+    {
+      IF.conservativeResize(2*IF.rows()+1,2);
+      CP.conservativeResize(IF.rows());
+    }
+    CP(num_if) = coplanar;
+    IF.row(num_if) << f1,f2;
+    num_if++;
+  };
+
+  // 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++)
+      {
+        // Find in opposite direction on jth face
+        if(F1(g,e) == d && F1(g,(e+1)%3) == s)
+        {
+          return c;
+        }
+      }
+    }
+    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++)
+    {
+      for(sg = 0;sg<3;sg++)
+      {
+        if(F1(g,sg) == F1(f,sf))
+        {
+          return true;
+        }
+      }
+    }
+    return false;
+  };
+
+  RowVector3S dummy;
+
+  igl::parallel_for(F2.rows(),[&](const int 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;
+    tree1.append_intersecting_leaves(box, candidates);
+    for(const auto * candidate : candidates)
+    {
+      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);
+          }
+        }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)
+            {
+              // But wait? Couldn't these be coplanar?
+              append_intersection(f1,f2,false);
+            }
+          }
+          
+        }
+      }
+      // This logic is confusing. 
+      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 = 
+          triangle_triangle_intersect(
+            V2.row(F2(f2,0)).template head<3>().eval(),
+            V2.row(F2(f2,1)).template head<3>().eval(),
+            V2.row(F2(f2,2)).template head<3>().eval(),
+            V1.row(F1(f1,0)).template head<3>().eval(),
+            V1.row(F1(f1,1)).template head<3>().eval(),
+            V1.row(F1(f1,2)).template head<3>().eval(),
+            coplanar);
+        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,coplanar);
+        }else if(tt_found_intersection)
+        {
+          found_intersection = true;
+          append_intersection(f1,f2,coplanar);
+        }
+      }
+      if(stinker) { printf("    %s\n",found_intersection? "☠️":"❌"); }
+      if(num_if && first_only) { return; }
+    }
+    if(num_if && first_only) { return; }
+  },1000);
+  IF.conservativeResize(num_if,2);
+  CP.conservativeResize(IF.rows());
+  return IF.rows();
+}
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedIF,
+  typename DerivedCP>
+IGL_INLINE bool igl::predicates::find_intersections(
+  const Eigen::MatrixBase<DerivedV1> & V1,
+  const Eigen::MatrixBase<DerivedF1> & F1,
+  const Eigen::MatrixBase<DerivedV2> & V2,
+  const Eigen::MatrixBase<DerivedF2> & F2,
+  const bool first_only,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedCP> & CP)
+{
+  igl::AABB<DerivedV1,3> tree1;
+  tree1.init(V1,F1);
+  return find_intersections(tree1,V1,F1,V2,F2,first_only,IF,CP);
+}
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedIF,
+  typename DerivedCP,
+  typename DerivedEV,
+  typename DerivedEE,
+  typename DerivedEI>
+IGL_INLINE bool igl::predicates::find_intersections(
+  const Eigen::MatrixBase<DerivedV1> & V1,
+  const Eigen::MatrixBase<DerivedF1> & F1,
+  const Eigen::MatrixBase<DerivedV2> & V2,
+  const Eigen::MatrixBase<DerivedF2> & F2,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedCP> & CP,
+  Eigen::PlainObjectBase<DerivedEV> & EV,
+  Eigen::PlainObjectBase<DerivedEE> & EE,
+  Eigen::PlainObjectBase<DerivedEI> & EI)
+{
+  igl::AABB<DerivedV1,3> tree1;
+  tree1.init(V1,F1);
+  return find_intersections(tree1,V1,F1,V2,F2,IF,CP,EV,EE,EI);
+}
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedIF,
+  typename DerivedCP,
+  typename DerivedEV,
+  typename DerivedEE,
+  typename DerivedEI>
+IGL_INLINE bool igl::predicates::find_intersections(
+  const igl::AABB<DerivedV1,3> & tree1,
+  const Eigen::MatrixBase<DerivedV1> & V1,
+  const Eigen::MatrixBase<DerivedF1> & F1,
+  const Eigen::MatrixBase<DerivedV2> & V2,
+  const Eigen::MatrixBase<DerivedF2> & F2,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedCP> & CP,
+  Eigen::PlainObjectBase<DerivedEV> & EV,
+  Eigen::PlainObjectBase<DerivedEE> & EE,
+  Eigen::PlainObjectBase<DerivedEI> & EI)
+{
+  if(!find_intersections(tree1,V1,F1,V2,F2,false,IF,CP)) { return false; }
+  std::vector<int> EI_vec = igl::find((CP.array()==false).eval());
+  igl::list_to_matrix(EI_vec,EI);
+  const auto IF_EI = IF(EI_vec,Eigen::all).eval();
+  igl::triangle_triangle_intersect(V1,F1,V2,F2,IF_EI,EV,EE);
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::predicates::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::Array<bool, -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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Array<bool, -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> >&);
+// generated by autoexplicit.sh
+template bool igl::predicates::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::Array<bool, -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::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::Array<bool, -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> >&);
+// generated by autoexplicit.sh
+template bool igl::predicates::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::Array<bool, -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, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/predicates/find_intersections.h

@@ -0,0 +1,118 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2024 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/
+#ifndef IGL_FAST_FIND_INTERSECTIONS_H
+#define IGL_FAST_FIND_INTERSECTIONS_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl 
+{
+  template <typename DerivedV, int DIM> class AABB;
+  namespace predicates
+  {
+    /// Identify triangles where two meshes interesect 
+    /// using AABBTree and igl::predicates::triangle_triangle_intersect.
+    ///
+    /// @param[in] tree1 AABB tree for the first mesh
+    /// @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] CP #IF list of whether the intersection is coplanar
+    ///
+    /// \see copyleft::cgal::intersect_other
+    template <
+      typename DerivedV1,
+      typename DerivedF1,
+      typename DerivedV2,
+      typename DerivedF2,
+      typename DerivedIF,
+      typename DerivedCP>
+    IGL_INLINE bool find_intersections(
+      const AABB<DerivedV1,3> & tree1,
+      const Eigen::MatrixBase<DerivedV1> & V1,
+      const Eigen::MatrixBase<DerivedF1> & F1,
+      const Eigen::MatrixBase<DerivedV2> & V2,
+      const Eigen::MatrixBase<DerivedF2> & F2,
+      const bool first_only,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedCP> & CP);
+    /// \overload
+    /// \brief Tree built internally.
+    template <
+      typename DerivedV1,
+      typename DerivedF1,
+      typename DerivedV2,
+      typename DerivedF2,
+      typename DerivedIF,
+      typename DerivedCP>
+    IGL_INLINE bool find_intersections(
+      const Eigen::MatrixBase<DerivedV1> & V1,
+      const Eigen::MatrixBase<DerivedF1> & F1,
+      const Eigen::MatrixBase<DerivedV2> & V2,
+      const Eigen::MatrixBase<DerivedF2> & F2,
+      const bool first_only,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedCP> & CP);
+    /// @param[out] EV #EV by 3 list of vertex positions of intersection segments
+    /// @param[out] EE #EE by 2 list of edge indices into EV
+    /// @param[out] EI #EI by 1 list of indices into rows IF indicating source of
+    ///   intersection.
+    template <
+      typename DerivedV1,
+      typename DerivedF1,
+      typename DerivedV2,
+      typename DerivedF2,
+      typename DerivedIF,
+      typename DerivedCP,
+      typename DerivedEV,
+      typename DerivedEE,
+      typename DerivedEI>
+    IGL_INLINE bool find_intersections(
+      const Eigen::MatrixBase<DerivedV1> & V1,
+      const Eigen::MatrixBase<DerivedF1> & F1,
+      const Eigen::MatrixBase<DerivedV2> & V2,
+      const Eigen::MatrixBase<DerivedF2> & F2,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedCP> & CP,
+      Eigen::PlainObjectBase<DerivedEV> & EV,
+      Eigen::PlainObjectBase<DerivedEE> & EE,
+      Eigen::PlainObjectBase<DerivedEI> & EI);
+  template <
+    typename DerivedV1,
+    typename DerivedF1,
+    typename DerivedV2,
+    typename DerivedF2,
+    typename DerivedIF,
+    typename DerivedCP,
+    typename DerivedEV,
+    typename DerivedEE,
+    typename DerivedEI>
+  IGL_INLINE bool find_intersections(
+    const igl::AABB<DerivedV1,3> & tree1,
+    const Eigen::MatrixBase<DerivedV1> & V1,
+    const Eigen::MatrixBase<DerivedF1> & F1,
+    const Eigen::MatrixBase<DerivedV2> & V2,
+    const Eigen::MatrixBase<DerivedF2> & F2,
+    Eigen::PlainObjectBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedCP> & CP,
+    Eigen::PlainObjectBase<DerivedEV> & EV,
+    Eigen::PlainObjectBase<DerivedEE> & EE,
+    Eigen::PlainObjectBase<DerivedEI> & EI);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "find_intersections.cpp"
+#endif
+
+#endif
+

include/igl/predicates/find_self_intersections.cpp

@@ -0,0 +1,54 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2024 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 "find_intersections.h"
+#include "find_self_intersections.h"
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedIF,
+  typename DerivedCP>
+IGL_INLINE bool igl::predicates::find_self_intersections(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const bool first_only,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedCP> & CP)
+{
+  // This is really just a wrapper around fast_find_intersections which will
+  // internally detect that V,F are the second set
+  return find_intersections( V,F,V,F,first_only,IF,CP);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedIF,
+  typename DerivedCP,
+  typename DerivedEV,
+  typename DerivedEE,
+  typename DerivedEI>
+bool igl::predicates::find_self_intersections(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedCP> & CP,
+  Eigen::PlainObjectBase<DerivedEV> & EV,
+  Eigen::PlainObjectBase<DerivedEE> & EE,
+  Eigen::PlainObjectBase<DerivedEI> & EI)
+{
+  return find_intersections( V,F,V,F,IF,CP,EV,EE,EI);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::predicates::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::Array<bool, -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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Array<bool, -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> >&);
+// generated by autoexplicit.sh
+template bool igl::predicates::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::Array<bool, -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, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Array<bool, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/predicates/find_self_intersections.h

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2024 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/
+#ifndef IGL_PREDICATES_FIND_SELF_INTERSECTIONS_H
+#define IGL_PREDICATES_FIND_SELF_INTERSECTIONS_H
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace predicates
+  {
+    /// Identify triangle-triangle interesections within the same mesh using
+    /// AABBTree and igl::predicates::triangle_triangle_intersect.
+    ///
+    /// @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] CP #IF list of whether the intersection is coplanar
+    ///
+    /// \see copyleft::cgal::SelfIntersectMesh
+    template <
+      typename DerivedV,
+      typename DerivedF,
+      typename DerivedIF,
+      typename DerivedCP>
+    IGL_INLINE bool find_self_intersections(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      const bool first_only,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedCP> & CP);
+    /// @param[out] EV #EV by 3 list of vertex positions of intersection segments
+    /// @param[out] EE #EE by 2 list of edge indices into EV
+    /// @param[out] EI #EI by 1 list of indices into rows IF indicating source of
+    ///   intersection.
+    template <
+      typename DerivedV,
+      typename DerivedF,
+      typename DerivedIF,
+      typename DerivedCP,
+      typename DerivedEV,
+      typename DerivedEE,
+      typename DerivedEI>
+    IGL_INLINE bool find_self_intersections(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedCP> & CP,
+      Eigen::PlainObjectBase<DerivedEV> & EV,
+      Eigen::PlainObjectBase<DerivedEE> & EE,
+      Eigen::PlainObjectBase<DerivedEI> & EI);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "find_self_intersections.cpp"
+#endif
+
+#endif
+

include/igl/predicates/triangle_triangle_intersect.cpp

@@ -0,0 +1,318 @@
+#include "triangle_triangle_intersect.h"
+#include "predicates.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+template <typename Vector3D>
+IGL_INLINE  bool igl::predicates::triangle_triangle_intersect(
+  const Vector3D & p1,
+  const Vector3D & q1,
+  const Vector3D & r1,
+  const Vector3D & p2,
+  const Vector3D & q2,
+  const Vector3D & r2,
+  bool & coplanar)
+{
+  coplanar = false;
+  // yet another translation of tri_tri_intersect.c  [Guigue & Devillers]
+  exactinit();
+  using Vector2D = Eigen::Matrix<typename Vector3D::Scalar,2,1>;
+
+  constexpr Orientation COPLANAR = igl::predicates::Orientation::COPLANAR;
+  constexpr Orientation NEGATIVE = igl::predicates::Orientation::NEGATIVE;
+  constexpr Orientation POSITIVE = igl::predicates::Orientation::POSITIVE;
+
+  // Determine for each vertex of one triangle if it's above, below or on the
+  // plane of the other triangle.
+
+  // SUB(v1,p2,r2)
+  // SUB(v2,q2,r2)
+  // CROSS(N2,v1,v2)
+  // SUB(v1,p1,r2)
+  // dp1 = DOT(v1,N2);
+  // dp1 = (p1-r2).dot( (p2-r2).cross(q2-r2) );
+  const Orientation dp1 = orient3d(p2,q2,r2,p1);
+  const Orientation dq1 = orient3d(p2,q2,r2,q1);
+  const Orientation dr1 = orient3d(p2,q2,r2,r1);
+
+  const auto same_non_coplanar = [&NEGATIVE,&COPLANAR,&POSITIVE](
+    const Orientation a, 
+    const Orientation b, 
+    const Orientation c)
+  {
+    return (a == POSITIVE && b == POSITIVE && c == POSITIVE) ||
+           (a == NEGATIVE && b == NEGATIVE && c == NEGATIVE);
+  };
+  if(same_non_coplanar(dp1,dq1,dr1)) { return false; }
+
+
+  const Orientation dp2 = orient3d(p1,q1,r1,p2);
+  const Orientation dq2 = orient3d(p1,q1,r1,q2);
+  const Orientation dr2 = orient3d(p1,q1,r1,r2);
+
+  // Theoreticaly, this should have already been fired above
+  if(same_non_coplanar(dp2,dq2,dr2)) { return false; }
+
+  const auto tri_tri_overlap_test_2d = [&NEGATIVE,&COPLANAR,&POSITIVE](
+      const Vector2D & p1,
+      const Vector2D & q1,
+      const Vector2D & r1,
+      const Vector2D & p2,
+      const Vector2D & q2,
+      const Vector2D & r2)->bool
+  {
+    const auto ccw_tri_tri_intersection_2d = [&NEGATIVE,&COPLANAR,&POSITIVE](
+      const Vector2D & p1,
+      const Vector2D & q1,
+      const Vector2D & r1,
+      const Vector2D & p2,
+      const Vector2D & q2,
+      const Vector2D & r2)->bool
+    {
+      const auto INTERSECTION_TEST_VERTEX = [&NEGATIVE,&COPLANAR,&POSITIVE](
+        const Vector2D & P1,
+        const Vector2D & Q1,
+        const Vector2D & R1,
+        const Vector2D & P2,
+        const Vector2D & Q2,
+        const Vector2D & R2)->bool
+      {
+        if (orient2d(R2,P2,Q1) != NEGATIVE)
+          if (orient2d(R2,Q2,Q1) != POSITIVE)
+            if (orient2d(P1,P2,Q1) == POSITIVE) {
+              if (orient2d(P1,Q2,Q1) != POSITIVE) return 1; 
+              else return 0;} else {
+                if (orient2d(P1,P2,R1) != NEGATIVE)
+                  if (orient2d(Q1,R1,P2) != NEGATIVE) return 1; 
+                  else return 0;
+                else return 0;}
+            else 
+              if (orient2d(P1,Q2,Q1) != POSITIVE)
+                if (orient2d(R2,Q2,R1) != POSITIVE)
+                  if (orient2d(Q1,R1,Q2) != NEGATIVE) return 1; 
+                  else return 0;
+                else return 0;
+              else return 0;
+          else
+            if (orient2d(R2,P2,R1) != NEGATIVE) 
+              if (orient2d(Q1,R1,R2) != NEGATIVE)
+                if (orient2d(P1,P2,R1) != NEGATIVE) return 1;
+                else return 0;
+              else 
+                if (orient2d(Q1,R1,Q2) != NEGATIVE) {
+                  if (orient2d(R2,R1,Q2) != NEGATIVE) return 1; 
+                  else return 0; }
+                else return 0; 
+            else  return 0; 
+      };
+      const auto INTERSECTION_TEST_EDGE = [&NEGATIVE,&COPLANAR,&POSITIVE](
+        const Vector2D & P1,
+        const Vector2D & Q1,
+        const Vector2D & R1,
+        const Vector2D & P2,
+        const Vector2D & Q2,
+        const Vector2D & R2)->bool
+      {
+        if (orient2d(R2,P2,Q1) != NEGATIVE) {
+          if (orient2d(P1,P2,Q1) != NEGATIVE) { 
+            if (orient2d(P1,Q1,R2) != NEGATIVE) return 1; 
+            else return 0;} else { 
+              if (orient2d(Q1,R1,P2) != NEGATIVE){ 
+                if (orient2d(R1,P1,P2) != NEGATIVE) return 1; else return 0;} 
+              else return 0; } 
+        } else {
+          if (orient2d(R2,P2,R1) != NEGATIVE) {
+            if (orient2d(P1,P2,R1) != NEGATIVE) {
+              if (orient2d(P1,R1,R2) != NEGATIVE) return 1;  
+              else {
+                if (orient2d(Q1,R1,R2) != NEGATIVE) return 1; else return 0;}}
+            else  return 0; }
+          else return 0; }
+      };
+      if ( orient2d(p2,q2,p1) != NEGATIVE ) {
+        if ( orient2d(q2,r2,p1) != NEGATIVE ) {
+          if ( orient2d(r2,p2,p1) != NEGATIVE ) return 1;
+          else return INTERSECTION_TEST_EDGE(p1,q1,r1,p2,q2,r2);
+        } else {  
+          if ( orient2d(r2,p2,p1) != NEGATIVE ) 
+      return INTERSECTION_TEST_EDGE(p1,q1,r1,r2,p2,q2);
+          else return INTERSECTION_TEST_VERTEX(p1,q1,r1,p2,q2,r2);}}
+      else {
+        if ( orient2d(q2,r2,p1) != NEGATIVE ) {
+          if ( orient2d(r2,p2,p1) != NEGATIVE ) 
+      return INTERSECTION_TEST_EDGE(p1,q1,r1,q2,r2,p2);
+          else  return INTERSECTION_TEST_VERTEX(p1,q1,r1,q2,r2,p2);}
+        else return INTERSECTION_TEST_VERTEX(p1,q1,r1,r2,p2,q2);}
+      return false;
+    };
+    if ( orient2d(p1,q1,r1) == NEGATIVE )
+      if ( orient2d(p2,q2,r2) == NEGATIVE )
+        return ccw_tri_tri_intersection_2d(p1,r1,q1,p2,r2,q2);
+      else
+        return ccw_tri_tri_intersection_2d(p1,r1,q1,p2,q2,r2);
+    else
+      if ( orient2d(p2,q2,r2) == NEGATIVE )
+        return ccw_tri_tri_intersection_2d(p1,q1,r1,p2,r2,q2);
+      else
+        return ccw_tri_tri_intersection_2d(p1,q1,r1,p2,q2,r2);
+
+  };
+
+  const auto coplanar_tri_tri3d = [&tri_tri_overlap_test_2d,&NEGATIVE,&COPLANAR,&POSITIVE](
+    const Vector3D & p1,
+    const Vector3D & q1,
+    const Vector3D & r1,
+    const Vector3D & p2,
+    const Vector3D & q2,
+    const Vector3D & r2)->bool
+  {
+    Vector3D normal_1 = (q1-p1).cross(r1-p1);
+    const auto n_x = ((normal_1[0]<0)?-normal_1[0]:normal_1[0]);
+    const auto n_y = ((normal_1[1]<0)?-normal_1[1]:normal_1[1]);
+    const auto n_z = ((normal_1[2]<0)?-normal_1[2]:normal_1[2]);
+    Vector2D P1,Q1,R1,P2,Q2,R2;
+    if (( n_x > n_z ) && ( n_x >= n_y )) {
+      // Project onto plane YZ
+
+        P1[0] = q1[2]; P1[1] = q1[1];
+        Q1[0] = p1[2]; Q1[1] = p1[1];
+        R1[0] = r1[2]; R1[1] = r1[1]; 
+      
+        P2[0] = q2[2]; P2[1] = q2[1];
+        Q2[0] = p2[2]; Q2[1] = p2[1];
+        R2[0] = r2[2]; R2[1] = r2[1]; 
+
+    } else if (( n_y > n_z ) && ( n_y >= n_x )) {
+      // Project onto plane XZ
+
+      P1[0] = q1[0]; P1[1] = q1[2];
+      Q1[0] = p1[0]; Q1[1] = p1[2];
+      R1[0] = r1[0]; R1[1] = r1[2]; 
+
+      P2[0] = q2[0]; P2[1] = q2[2];
+      Q2[0] = p2[0]; Q2[1] = p2[2];
+      R2[0] = r2[0]; R2[1] = r2[2]; 
+      
+    } else {
+      // Project onto plane XY
+
+      P1[0] = p1[0]; P1[1] = p1[1]; 
+      Q1[0] = q1[0]; Q1[1] = q1[1]; 
+      R1[0] = r1[0]; R1[1] = r1[1]; 
+      
+      P2[0] = p2[0]; P2[1] = p2[1]; 
+      Q2[0] = q2[0]; Q2[1] = q2[1]; 
+      R2[0] = r2[0]; R2[1] = r2[1]; 
+    }
+
+    return tri_tri_overlap_test_2d(P1,Q1,R1,P2,Q2,R2);
+    exit(1);
+    return false;
+  };
+
+  const auto TRI_TRI_3D = [&coplanar_tri_tri3d,&coplanar,&NEGATIVE,&COPLANAR,&POSITIVE](
+    const Vector3D & p1,
+    const Vector3D & q1,
+    const Vector3D & r1,
+    const Vector3D & p2,
+    const Vector3D & q2,
+    const Vector3D & r2,
+    const Orientation dp2,
+    const Orientation dq2,
+    const Orientation dr2)->bool
+  {
+    const auto CHECK_MIN_MAX = [&NEGATIVE,&COPLANAR,&POSITIVE](
+      const Vector3D & p1,
+      const Vector3D & q1,
+      const Vector3D & r1,
+      const Vector3D & p2,
+      const Vector3D & q2,
+      const Vector3D & r2)->bool
+    {
+       if (orient3d(p2,p1,q1,q2) == POSITIVE) { return false; }
+       if (orient3d(p2,r1,p1,r2) == POSITIVE) { return false; }
+       return true;
+    };
+
+
+
+    if (dp2 == POSITIVE) { 
+       if (dq2 == POSITIVE) return CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2) ;
+       else if (dr2 == POSITIVE) return CHECK_MIN_MAX(p1,r1,q1,q2,r2,p2);
+       else return CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2);
+    } else if (dp2 == NEGATIVE) { 
+      if (dq2 == NEGATIVE) return CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2);
+      else if (dr2 == NEGATIVE) return CHECK_MIN_MAX(p1,q1,r1,q2,r2,p2);
+      else return CHECK_MIN_MAX(p1,r1,q1,p2,q2,r2);
+    } else { 
+      if (dq2 == NEGATIVE) { 
+        if (dr2 == POSITIVE || dr2 == COPLANAR)  return CHECK_MIN_MAX(p1,r1,q1,q2,r2,p2);
+        else return CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2);
+      } 
+      else if (dq2 == POSITIVE) { 
+        if (dr2 == POSITIVE) return CHECK_MIN_MAX(p1,r1,q1,p2,q2,r2);
+        else  return CHECK_MIN_MAX(p1,q1,r1,q2,r2,p2);
+      } 
+      else  { 
+        if (dr2 == POSITIVE) return CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2);
+        else if (dr2 == NEGATIVE) return CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2);
+        else
+        {
+          coplanar = coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2);
+          return coplanar;
+        }
+      }
+    }
+  };
+
+  if (dp1 == POSITIVE) {
+    if (dq1 == POSITIVE) return TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2);
+    else if (dr1 == POSITIVE) return TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2)  ;
+    else return TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2);
+  } else if (dp1 == NEGATIVE) {
+    if (dq1 == NEGATIVE) return TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2);
+    else if (dr1 == NEGATIVE) return TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2);
+    else return TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2);
+  } else {
+    if (dq1 == NEGATIVE) {
+      // why COPLANAR here? It seems so haphazard.
+      if (dr1 == POSITIVE || dr1 == COPLANAR) return TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2);
+      else return TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2);
+    }
+    else if (dq1 == POSITIVE) {
+      if (dr1 == POSITIVE) return TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2);
+      else return TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2);
+    }
+    else  {
+      if (dr1 == POSITIVE) return TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2);
+      else if (dr1 == NEGATIVE) return TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2);
+      else
+      {
+        coplanar = coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2);
+        return coplanar;
+      }
+    }
+  }
+
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// template using Eigen::Vector3d.
+template bool igl::predicates::triangle_triangle_intersect<Eigen::Vector3d>(
+    const Eigen::Vector3d & p1,
+    const Eigen::Vector3d & q1,
+    const Eigen::Vector3d & r1,
+    const Eigen::Vector3d & p2,
+    const Eigen::Vector3d & q2,
+    const Eigen::Vector3d & r2,
+    bool & coplanar);
+template bool igl::predicates::triangle_triangle_intersect<Eigen::RowVector3d>(
+    const Eigen::RowVector3d & p1,
+    const Eigen::RowVector3d & q1,
+    const Eigen::RowVector3d & r1,
+    const Eigen::RowVector3d & p2,
+    const Eigen::RowVector3d & q2,
+    const Eigen::RowVector3d & r2,
+    bool & coplanar);
+#endif

include/igl/predicates/triangle_triangle_intersect.h

@@ -0,0 +1,25 @@
+#ifndef IGL_PREDICATES_TRIANGLE_TRIANGLE_INTERSECT_H
+#define IGL_PREDICATES_TRIANGLE_TRIANGLE_INTERSECT_H
+#include "../igl_inline.h"
+
+namespace igl
+{
+  namespace predicates
+  {
+    template <typename Vector3D>
+    IGL_INLINE  bool triangle_triangle_intersect(
+      const Vector3D & a1,
+      const Vector3D & a2,
+      const Vector3D & a3,
+      const Vector3D & b1,
+      const Vector3D & b2,
+      const Vector3D & b3,
+      bool & coplanar);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "triangle_triangle_intersect.cpp"
+#endif
+
+#endif 

include/igl/tri_tri_intersect.cpp

@@ -63,6 +63,7 @@
 #define IGL_TRI_TRI_INTERSECT_CPP
 
 #include "tri_tri_intersect.h"
+#include "EPS.h"
 #include <Eigen/Geometry>
 
 // helper functions
@@ -84,49 +85,18 @@ IGL_INLINE bool ccw_tri_tri_intersection_2d(
   const Eigen::MatrixBase<DerivedP1> &p1, const Eigen::MatrixBase<DerivedQ1> &q1, const Eigen::MatrixBase<DerivedR1> &r1,
   const Eigen::MatrixBase<DerivedP2> &p2, const Eigen::MatrixBase<DerivedQ2> &q2, const Eigen::MatrixBase<DerivedR2> &r2);  
 
-
-/* some 3D macros */
-
-// #define _IGL_CROSS(dest,v1,v2)  \
-//                dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
-//                dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
-//                dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
-
-// #define _IGL_DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
-
- 
-// #define _IGL_SUB(dest,v1,v2) dest[0]=v1[0]-v2[0]; \
-//                         dest[1]=v1[1]-v2[1]; \
-//                         dest[2]=v1[2]-v2[2]; 
-
-
-// #define _IGL_SCALAR(dest,alpha,v) dest[0] = alpha * v[0]; \
-//                              dest[1] = alpha * v[1]; \
-//                              dest[2] = alpha * v[2];
-
-
-// #define _IGL_CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2) {\
-//   _IGL_SUB(v1,p2,q1)\
-//   _IGL_SUB(v2,p1,q1)\
-//   _IGL_CROSS(N1,v1,v2)\
-//   _IGL_SUB(v1,q2,q1)\
-//   if (_IGL_DOT(v1,N1) > 0.0) return false;\
-//   _IGL_SUB(v1,p2,p1)\
-//   _IGL_SUB(v2,r1,p1)\
-//   _IGL_CROSS(N1,v1,v2)\
-//   _IGL_SUB(v1,r2,p1) \
-//   if (_IGL_DOT(v1,N1) > 0.0) return false;\
-//   else return true; }
-
   template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
             typename DerivedP2,typename DerivedQ2,typename DerivedR2>
   inline bool _IGL_CHECK_MIN_MAX(
     const Eigen::MatrixBase<DerivedP1> &p1, const Eigen::MatrixBase<DerivedQ1> &q1, const Eigen::MatrixBase<DerivedR1> &r1,
     const Eigen::MatrixBase<DerivedP2> &p2, const Eigen::MatrixBase<DerivedQ2> &q2, const Eigen::MatrixBase<DerivedR2> &r2)
   {
-    auto v1=p2-q1;
-    auto v2=p1-q1;
-    auto N1=v1.cross(v2);
+    using Scalar    = typename DerivedP1::Scalar;
+    using RowVector = typename Eigen::Matrix<Scalar, 1, 3>;
+
+    RowVector v1=p2-q1;
+    RowVector v2=p1-q1;
+    RowVector N1=v1.cross(v2);
     v1=q2-q1;
 
     if (v1.dot(N1) > 0.0) return false;
@@ -152,7 +122,7 @@ IGL_INLINE bool ccw_tri_tri_intersection_2d(
     const Eigen::MatrixBase<DerivedP1> &p1, const Eigen::MatrixBase<DerivedQ1> &q1, const Eigen::MatrixBase<DerivedR1> &r1,
     const Eigen::MatrixBase<DerivedP2> &p2, const Eigen::MatrixBase<DerivedQ2> &q2, const Eigen::MatrixBase<DerivedR2> &r2,
     DP2 dp2, DQ2 dq2,DR2 dr2,
-    const Eigen::MatrixBase<DerivedP2> &N1)
+    const Eigen::MatrixBase<DerivedN1> &N1)
   {
     if (dp2 > 0.0) { 
       if (dq2 > 0.0) return _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2);
@@ -179,31 +149,7 @@ IGL_INLINE bool ccw_tri_tri_intersection_2d(
   }
 
 
-// #define _IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2) { \
-//   if (dp2 > 0.0) { \
-//      if (dq2 > 0.0) _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2) \
-//      else if (dr2 > 0.0) _IGL_CHECK_MIN_MAX(p1,r1,q1,q2,r2,p2)\
-//      else _IGL_CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2) }\
-//   else if (dp2 < 0.0) { \
-//     if (dq2 < 0.0) _IGL_CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2)\
-//     else if (dr2 < 0.0) _IGL_CHECK_MIN_MAX(p1,q1,r1,q2,r2,p2)\
-//     else _IGL_CHECK_MIN_MAX(p1,r1,q1,p2,q2,r2)\
-//   } else { \
-//     if (dq2 < 0.0) { \
-//       if (dr2 >= 0.0)  _IGL_CHECK_MIN_MAX(p1,r1,q1,q2,r2,p2)\
-//       else _IGL_CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2)\
-//     } \
-//     else if (dq2 > 0.0) { \
-//       if (dr2 > 0.0) _IGL_CHECK_MIN_MAX(p1,r1,q1,p2,q2,r2)\
-//       else  _IGL_CHECK_MIN_MAX(p1,q1,r1,q2,r2,p2)\
-//     } \
-//     else  { \
-//       if (dr2 > 0.0) _IGL_CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2)\
-//       else if (dr2 < 0.0) _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2)\
-//       else return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);\
-//      }}}
   
-//}
 } //igl
 
 } // internal
@@ -224,19 +170,19 @@ IGL_INLINE bool igl::tri_tri_overlap_test_3d(
   const Eigen::MatrixBase<DerivedR2> &  r2)
 {
   using Scalar    = typename DerivedP1::Scalar;
-  //using RowVector = typename Eigen::Matrix<Scalar, 1, 3>;
+  using RowVector = typename Eigen::Matrix<Scalar, 1, 3>;
  
   Scalar dp1, dq1, dr1, dp2, dq2, dr2;
-  //RowVector v1, v2;
-  //RowVector N1, N2; 
+  RowVector v1, v2;
+  RowVector N1, N2; 
   
   /* Compute distance signs  of p1, q1 and r1 to the plane of
      triangle(p2,q2,r2) */
 
 
-  auto v1=p2-r2;
-  auto v2=q2-r2;
-  auto N2=v1.cross(v2);
+  v1=p2-r2;
+  v2=q2-r2;
+  N2=v1.cross(v2);
 
   v1=p1-r2;
   dp1 = v1.dot(N2);
@@ -252,7 +198,7 @@ IGL_INLINE bool igl::tri_tri_overlap_test_3d(
 
   v1=q1-p1;
   v2=r1-p1;
-  auto N1=v1.cross(v2);
+  N1=v1.cross(v2);
 
   v1=p2-r1;
   dp2 = v1.dot(N1);
@@ -265,28 +211,27 @@ IGL_INLINE bool igl::tri_tri_overlap_test_3d(
 
   /* Permutation in a canonical form of T1's vertices */
 
-
   if (dp1 > 0.0) {
-    if (dq1 > 0.0) return _IGL_TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2,N1);
-    else if (dr1 > 0.0) return _IGL_TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2,N1);
-    else return _IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2,N1);
+    if (dq1 > 0.0) return internal::_IGL_TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2,N1);
+    else if (dr1 > 0.0) return internal::_IGL_TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2,N1);
+    else return internal::_IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2,N1);
   } else if (dp1 < 0.0) {
-    if (dq1 < 0.0) return _IGL_TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2,N1);
-    else if (dr1 < 0.0) return _IGL_TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2,N1);
-    else return _IGL_TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2,N1);
+    if (dq1 < 0.0) return internal::_IGL_TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2,N1);
+    else if (dr1 < 0.0) return internal::_IGL_TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2,N1);
+    else return internal::_IGL_TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2,N1);
   } else {
     if (dq1 < 0.0) {
-      if (dr1 >= 0.0) return _IGL_TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2,N1);
-      else return _IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2,N1);
+      if (dr1 >= 0.0) return internal::_IGL_TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2,N1);
+      else return internal::_IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2,N1);
     }
     else if (dq1 > 0.0) {
-      if (dr1 > 0.0) return _IGL_TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2,N1);
-      else return _IGL_TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2,N1);
+      if (dr1 > 0.0) return internal::_IGL_TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2,N1);
+      else return internal::_IGL_TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2,N1);
     }
     else  {
-      if (dr1 > 0.0) return _IGL_TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2,N1);
-      else if (dr1 < 0.0) return _IGL_TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2,N1);
-      else return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
+      if (dr1 > 0.0) return internal::_IGL_TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2,N1);
+      else if (dr1 < 0.0) return internal::_IGL_TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2,N1);
+      else return internal::coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
     }
   }
 };
@@ -667,8 +612,24 @@ IGL_INLINE bool igl::tri_tri_intersection_test_3d(
   dq2 = v1.dot(N1);
   v1=r2-r1;
   dr2 = v1.dot(N1);
+
   
   if (((dp2 * dq2) > 0.0) && ((dp2 * dr2) > 0.0)) return false;
+  // Alec: it's hard to believe this will ever be perfectly robust, but checking
+  // 1e-22 against zero seems like a recipe for bad logic.
+  // Switching all these 0.0s to epsilons makes other tests fail. My claim is
+  // that the series of logic below is a bad way of determining coplanarity, so
+  // instead just check for it right away.
+  const Scalar eps = igl::EPS<Scalar>();
+  using std::abs;
+  if(
+      abs(dp1) < eps && abs(dq1) < eps && abs(dr1) < eps &&
+      abs(dp2) < eps && abs(dq2) < eps && abs(dr2) < eps)
+  { 
+    coplanar = true;
+    return internal::coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
+  };
+
 
   // Permutation in a canonical form of T1's vertices
   if (dp1 > 0.0) {
@@ -693,7 +654,7 @@ IGL_INLINE bool igl::tri_tri_intersection_test_3d(
       if (dr1 > 0.0) return internal::_IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2,coplanar,source,target,N1,N2);
       else if (dr1 < 0.0) return internal::_IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2,coplanar,source,target,N1,N2);
       else {
-        // triangles are co-planar
+        // triangles are co-planar (should have been caught above).
 
         coplanar = true;
         return internal::coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
@@ -894,4 +855,8 @@ template bool igl::tri_tri_intersection_test_3d<Eigen::Block<Eigen::Matrix<doubl
 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> >&);
+
+template bool igl::tri_tri_overlap_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::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&);
+template bool igl::tri_tri_overlap_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::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&);
+template bool igl::tri_tri_overlap_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::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&);
 #endif

include/igl/triangle_triangle_intersect.cpp

@@ -4,6 +4,9 @@
 #include "tri_tri_intersect.h"
 #include <Eigen/Geometry>
 #include <stdio.h>
+#include <igl/unique_edge_map.h>
+#include <igl/barycentric_coordinates.h>
+#include <unordered_map>
 
 //#define IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
 #ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
@@ -104,7 +107,15 @@ IGL_INLINE bool igl::triangle_triangle_intersect(
     {
       bool coplanar;
       Eigen::RowVector3d i1,i2;
-      found_intersection = igl::tri_tri_intersection_test_3d(
+      found_intersection = 
+        igl::tri_tri_overlap_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>()) && 
+        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>(),
@@ -128,8 +139,190 @@ IGL_INLINE bool igl::triangle_triangle_intersect(
   return found_intersection;
 }
 
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedIF,
+  typename DerivedEV,
+  typename DerivedEE>
+void igl::triangle_triangle_intersect(
+    const Eigen::MatrixBase<DerivedV> & V1,
+    const Eigen::MatrixBase<DerivedF> & F1,
+    const Eigen::MatrixBase<DerivedV> & V2,
+    const Eigen::MatrixBase<DerivedF> & F2,
+    const Eigen::MatrixBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedEV> & EV,
+    Eigen::PlainObjectBase<DerivedEE> & EE)
+{
+  using Scalar = typename DerivedEV::Scalar;
+  using RowVector3S = Eigen::Matrix<Scalar,1,3>;
+  // We were promised that IF is a list of non-coplanar non-degenerately
+  // intersecting triangle pairs. This implies that the set of intersection is a
+  // line-segment whose endpoints are defined by edge-triangle intersections.
+  // 
+  // Each edge-triangle intersection will be stored in a map (e,f)
+  //
+  std::unordered_map<std::int64_t,int> uf_to_ev;
+  Eigen::VectorXi EMAP1;
+  Eigen::MatrixXi uE1;
+  {
+    Eigen::VectorXi uEE,uEC;
+    Eigen::MatrixXi E;
+    igl::unique_edge_map(F1,E,uE1,EMAP1,uEE,uEC);
+  }
+  Eigen::VectorXi EMAP2_cpy;
+  Eigen::MatrixXi uE2_cpy;
+  const bool self = &F1 == &F2;
+  if(!self)
+  {
+    Eigen::VectorXi uEE,uEC;
+    Eigen::MatrixXi E;
+    igl::unique_edge_map(F2,E,uE2_cpy,EMAP2_cpy,uEE,uEC);
+  }
+  const Eigen::VectorXi & EMAP2 = self?EMAP1:EMAP2_cpy;
+  const Eigen::MatrixXi & uE2 = self?uE1:uE2_cpy;
+
+  int num_ev = 0;
+  EE.resize(IF.rows(),2);
+  for(int i = 0; i<IF.rows(); i++)
+  {
+    // Just try all 6 edges
+    Eigen::Matrix<double,6,6> B;
+    Eigen::Matrix<double,6,3> X;
+    Eigen::Matrix<int,6,2> uf;
+    for(int p = 0;p<2;p++)
+    {
+      const auto consider_edges = [&B,&X,&uf]
+        (const int p,
+         const int f1, 
+         const int f2,
+         const Eigen::MatrixBase<DerivedV> & V1,
+         const Eigen::MatrixBase<DerivedF> & F1,
+         const Eigen::VectorXi & EMAP1,
+         const Eigen::MatrixXi & uE1,
+         const Eigen::MatrixBase<DerivedV> & V2,
+         const Eigen::MatrixBase<DerivedF> & F2)
+      {
+        for(int e1 = 0;e1<3;e1++)
+        {
+          // intersect edge (ij) opposite vertex F(f1,e1) with triangle ABC of F(f2,:)
+          const int u1 = EMAP1(f1 +(EMAP1.size()/3)*e1);
+          uf.row(p*3+e1) << u1,f2;
+          const int i = uE1(u1,0);
+          const int j = uE1(u1,1);
+          // Just copy.
+          const RowVector3S Vi = V1.row(i);
+          const RowVector3S Vj = V1.row(j);
+          const RowVector3S VA = V2.row(F2(f2,0));
+          const RowVector3S VB = V2.row(F2(f2,1));
+          const RowVector3S VC = V2.row(F2(f2,2));
+          // Find intersection of line (Vi,Vj) with plane of triangle (A,B,C)
+          const RowVector3S n = (VB-VA).template head<3>().cross((VC-VA).template head<3>());
+          const Scalar d = n.dot(VA);
+          const Scalar t = (d - n.dot(Vi))/(n.dot(Vj-Vi));
+          const RowVector3S x = Vi + t*(Vj-Vi);
+          
+          // Get barycenteric coordinates (possibly negative of X)
+          RowVector3S b1;
+          igl::barycentric_coordinates(x,VA,VB,VC,b1);
+          B.row(p*3+e1).head<3>() = b1;
+          RowVector3S b2;
+          igl::barycentric_coordinates(x,
+              V1.row(F1(f1,0)).template head<3>().eval(),
+              V1.row(F1(f1,1)).template head<3>().eval(),
+              V1.row(F1(f1,2)).template head<3>().eval(),
+              b2);
+          B.row(p*3+e1).tail<3>() = b2;
+
+          X.row(p*3+e1) = x;
+        }
+      };
+
+
+      const int f1 = IF(i,p);
+      const int f2 = IF(i,(p+1)%2);
+      consider_edges(p,f1,f2,
+        p==0?   V1:V2,
+        p==0?   F1:F2,
+        p==0?EMAP1:EMAP2,
+        p==0?  uE1:uE2,
+        p==0?   V2:V1,
+        p==0?   F2:F1);
+    }
+
+    // Find the two rows in B with the largest-smallest element
+    int j1,j2;
+    {
+      double b_min1 = -std::numeric_limits<double>::infinity();
+      double b_min2 = -std::numeric_limits<double>::infinity();
+      for(int j = 0;j<6;j++)
+      {
+        // It's not clear that using barycentric coordinates is better than
+        // point_simplex distance (though that requires inequalities).
+        const double bminj = B.row(j).minCoeff();
+        if(bminj > b_min1)
+        {
+          b_min2 = b_min1;
+          j2 = j1;
+          b_min1 = bminj;
+          j1 = j;
+        }else if(bminj > b_min2)
+        {
+          b_min2 = bminj;
+          j2 = j;
+        }
+      }
+    }
+
+    const auto append_or_find = [&](
+      int p, int u, int f, const RowVector3S & x,
+      std::unordered_map<std::int64_t,int> & uf_to_ev)->int
+    {
+      const std::int64_t key = (std::int64_t)u + ((std::int64_t)f << 32) + ((std::int64_t)p << 63);
+      if(uf_to_ev.find(key) == uf_to_ev.end())
+      {
+        if(num_ev == EV.rows())
+        {
+          EV.conservativeResize(2*EV.rows()+1,3);
+        }
+        EV.row(num_ev) = x;
+        uf_to_ev[key] = num_ev;
+        num_ev++;
+      }
+      return uf_to_ev[key];
+    };
+
+    EE.row(i) <<
+      append_or_find(j1>=3,uf(j1,0),uf(j1,1),X.row(j1),uf_to_ev),
+      append_or_find(j2>=3,uf(j2,0),uf(j2,1),X.row(j2),uf_to_ev);
+  }
+  EV.conservativeResize(num_ev,3);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedIF,
+  typename DerivedEV,
+  typename DerivedEE>
+void igl::triangle_triangle_intersect(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedIF> & IF,
+    Eigen::PlainObjectBase<DerivedEV> & EV,
+    Eigen::PlainObjectBase<DerivedEE> & EE)
+{
+  // overload will take care of detecting reference equality
+  return triangle_triangle_intersect(V,F,V,F,IF,EV,EE);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void 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<double, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void 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<double, -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&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 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

@@ -41,7 +41,8 @@ namespace igl
   /// 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
+  /// \see edge_flaps, tri_tri_intersect,
+  /// predicates::triangle_triangle_intersect
   template <
     typename DerivedV,
     typename DerivedF,
@@ -59,6 +60,46 @@ namespace igl
     const int c,
     const Eigen::MatrixBase<Derivedp> & p,
     const int g);
+  /// Assuming that we've computed a list of non-coplanar intersecting triangle
+  /// pairs between (V1,F1) and (V2,F2), compute the 
+  ///
+  /// @param[in] V1 #V1 by 3 list of vertex positions
+  /// @param[in] F1 #F1 by 3 list of triangle indices into rows of V1
+  /// @param[in] V2 #V2 by 3 list of vertex positions
+  /// @param[in] F2 #F2 by 3 list of triangle indices into rows of V2
+  /// @param[in] IF #IF by 2 list of intersecting triangle pairs so that IF(i,:)
+  ///   [f1 f2] means that the f1th triangle in F1 intersects the f2th triangle in
+  ///   F2 and that they're non-coplanar
+  /// @param[out] EV #EV by 3 list of vertex positions of intersection segments
+  /// @param[out] EE #EE by 2 list of edge indices into EV
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedIF,
+    typename DerivedEV,
+    typename DerivedEE>
+  IGL_INLINE void triangle_triangle_intersect(
+      const Eigen::MatrixBase<DerivedV> & V1,
+      const Eigen::MatrixBase<DerivedF> & F1,
+      const Eigen::MatrixBase<DerivedV> & V2,
+      const Eigen::MatrixBase<DerivedF> & F2,
+      const Eigen::MatrixBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedEV> & EV,
+      Eigen::PlainObjectBase<DerivedEE> & EE);
+  /// \overload
+  /// \brief Wrapper for triangle_triangle_intersect with self-intersections
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedIF,
+    typename DerivedEV,
+    typename DerivedEE>
+  IGL_INLINE void triangle_triangle_intersect(
+      const Eigen::MatrixBase<DerivedV> & V,
+      const Eigen::MatrixBase<DerivedF> & F,
+      const Eigen::MatrixBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedEV> & EV,
+      Eigen::PlainObjectBase<DerivedEE> & EE);
 }
 
 #ifndef IGL_STATIC_LIBRARY

tests/include/igl/fast_find_intersections.cpp

@@ -1,53 +0,0 @@
-#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("fast_find_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::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( (IF.col(1).array()==0).all() );
-  
-  // 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

@@ -1,129 +0,0 @@
-#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/predicates/find_intersections.cpp

@@ -0,0 +1,82 @@
+#include "test_common.h"
+#include <igl/predicates/find_intersections.h>
+#include <igl/combine.h>
+#include <igl/triangle_triangle_intersect.h>
+#include <igl/unique.h>
+#include <igl/sortrows.h>
+#include <igl/matlab_format.h>
+#include <iostream>
+
+
+
+TEST_CASE("find_intersections: cube-triangle", "[igl/predicates]")
+{
+  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 IF;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  igl::predicates::find_intersections(V,F,Vp,Fp,false,IF,CP);
+  {
+    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( (IF.col(1).array()==0).all() );
+}
+
+TEST_CASE("find_intersections: extract", "[igl/predicates]")
+{
+  Eigen::MatrixXd V1(4,3);
+  V1 << 0,0,0,
+       1,0,0,
+       1,1,0,
+       0,1,0;
+  Eigen::MatrixXi F1(2,3);
+  F1 << 0,1,2,
+       0,2,3;
+  Eigen::MatrixXd V2(3,3);
+  V2 << 1,0,1,
+       0,1,1,
+       0.5,0.5,-1;
+  Eigen::MatrixXi F2(1,3);
+  F2 << 0,1,2;
+
+  Eigen::MatrixXi IF,EE; Eigen::MatrixXd EV;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  bool found = igl::predicates::find_intersections(V1,F1,V2,F2,false,IF,CP);
+  REQUIRE( (CP==false).all() );
+  REQUIRE( found );
+  REQUIRE( IF.rows() == 2);
+  Eigen::MatrixXi IF_gt(2,2);
+  IF_gt << 
+    0,0,
+    1,0;
+  test_common::assert_near_rows(IF,IF_gt,0);
+
+  igl::triangle_triangle_intersect(V1,F1,V2,F2,IF,EV,EE);
+  Eigen::MatrixXd EV_gt(3,3);
+  EV_gt<<
+    0.5 ,0.5,0,
+    0.25,0.75,0,
+    0.75,0.25,0;
+  test_common::assert_near_rows(EV,EV_gt,0);
+  REQUIRE( EE.rows() == 2);
+}

tests/include/igl/predicates/find_self_intersections.cpp

@@ -0,0 +1,228 @@
+#include "test_common.h"
+#include <igl/predicates/find_self_intersections.h>
+#include <igl/upsample.h>
+#include <igl/triangle_triangle_intersect.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("find_self_intersections: cube", "[igl/predicates]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+
+  igl::read_triangle_mesh(test_common::data_path("cube.obj"), V, F);
+
+  Eigen::MatrixXi IF;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  REQUIRE( !igl::predicates::find_self_intersections(V,F,true,IF,CP) );
+}
+
+TEST_CASE("find_self_intersections: cube-triangle", "[igl/predicates]")
+{
+  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;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  REQUIRE( igl::predicates::find_self_intersections(V_,F_,false,IF,CP) );
+  {
+    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() );
+  }
+}
+
+
+TEST_CASE("find_self_intersections: rose", "[igl/predicates]")
+{
+  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;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  REQUIRE( igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  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("find_self_intersections: shared-edge", "[igl/predicates]")
+{
+  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;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  REQUIRE(!igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==0 );
+  REQUIRE( CP.rows()==0 );
+  V.row(3) << 2.0,0,0;
+  REQUIRE( igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==1 );
+  REQUIRE( CP.rows()==1 );
+  REQUIRE( CP(0) );
+}
+
+TEST_CASE("find_self_intersections: coplanar", "[igl/predicates]")
+{
+  Eigen::MatrixXd V(6,3);
+  V<< 
+  0.30947001278400399,0.80785250663757346,0.47595100104808802,
+  0.299046009778976,0.79801350831985496,0.47843150794506101,
+  0.32418551295995701,0.79794725775718722,0.48203899711370451,
+  0.30425801128148999,0.80293300747871421,0.47719125449657451,
+  0.31897351145744302,0.79302775859832797,0.48327925056219101,
+  0.31376150995492902,0.78810825943946872,0.4845195040106775;
+  Eigen::MatrixXi F(2,3);
+  F<< 
+    0,4,2,
+    1,5,3;
+
+  Eigen::MatrixXi IF;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  REQUIRE( !igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==0 );
+  REQUIRE( CP.rows()==0 );
+}
+
+TEST_CASE("find_self_intersections: non-intersecting", "[igl/predicates]")
+{
+  Eigen::MatrixXd V(6,3);
+  V<< 
+    0.39234799146652199,0.38443601131439198,0.44744500517845198,
+    0.38924700021743752,0.385637506842613,0.45762349665164948,
+    0.38700349628925301,0.38789276033639897,0.45634675025939975,
+    0.39079749584197976,0.38503675907850249,0.45253425091505073,
+    0.38769650459289529,0.38623825460672351,0.46271274238824822,
+    0.39279749989509577,0.37299175560474374,0.45553924888372399;
+  Eigen::MatrixXi F(2,3);
+  F<< 
+    0,3,2,
+    1,5,4;
+
+  Eigen::MatrixXi IF;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  REQUIRE( !igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==0 );
+  REQUIRE( CP.rows()==0 );
+}
+
+TEST_CASE("find_self_intersections: upsampled-knight", "[igl/predicates]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::read_triangle_mesh(test_common::data_path("decimated-knight.obj"),V,F);
+
+  Eigen::MatrixXi IF;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+
+  REQUIRE( !igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==0 );
+  REQUIRE( CP.rows()==0 );
+
+  igl::upsample(Eigen::MatrixXd(V),Eigen::MatrixXi(F),V,F);
+
+
+  REQUIRE( !igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==0 );
+  REQUIRE( CP.rows()==0 );
+
+  igl::upsample(Eigen::MatrixXd(V),Eigen::MatrixXi(F),V,F);
+
+  REQUIRE( !igl::predicates::find_self_intersections(V,F,false,IF,CP));
+  REQUIRE( IF.rows()==0 );
+  REQUIRE( CP.rows()==0 );
+
+}
+
+TEST_CASE("find_self_intersections: extract", "[igl/predicates]")
+{
+  Eigen::MatrixXd V(4+3,3);
+  V << 0,0,0,
+       1,0,0,
+       1,1,0,
+       0,1,0,
+       1,0,1,
+       0,1,1,
+       0.5,0.5,-1;
+  Eigen::MatrixXi F(3,3);
+  F << 0,1,2,
+       0,2,3,
+    4,5,6;
+  Eigen::MatrixXi IF,EE; Eigen::MatrixXd EV;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  bool found = igl::predicates::find_self_intersections(V,F,false,IF,CP);
+
+  igl::triangle_triangle_intersect(V,F,IF,EV,EE);
+  REQUIRE( (CP==false).all() );
+  REQUIRE( found );
+  REQUIRE( IF.rows() == 2);
+  Eigen::MatrixXi IF_gt(2,2);
+  IF_gt << 
+    0,2,
+    1,2;
+  test_common::assert_near_rows(IF,IF_gt,0);
+
+  Eigen::MatrixXd EV_gt(3,3);
+  EV_gt<<
+    0.5 ,0.5,0,
+    0.25,0.75,0,
+    0.75,0.25,0;
+  test_common::assert_near_rows(EV,EV_gt,0);
+  REQUIRE( EE.rows() == 2);
+  
+}

tests/include/igl/predicates/triangle_triangle_intersect.cpp

@@ -0,0 +1,90 @@
+#include "test_common.h"
+#include <igl/predicates/triangle_triangle_intersect.h>
+
+TEST_CASE("triangle_triangle_intersect intersect", "[igl/predicates]")
+{
+  // 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);
+
+    bool coplanar;
+    REQUIRE( igl::predicates::triangle_triangle_intersect(p1,q1,r1, p2,q2,r2,coplanar));
+  }
+}
+
+TEST_CASE("triangle_triangle_intersect not intersect", "[igl/predicates]")
+{
+  // 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);
+
+  bool coplanar;
+  REQUIRE( !igl::predicates::triangle_triangle_intersect(p1,q1,r1,p2,q2,r2,coplanar) );
+}
+
+
+TEST_CASE("triangle_triangle_intersect coplanar", "[igl/predicates]")
+{
+  // 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)
+  bool coplanar;
+  REQUIRE( igl::predicates::triangle_triangle_intersect(p1,q1,r1,p2,q2,r2,coplanar) );
+  REQUIRE( coplanar );
+}
+
+TEST_CASE("triangle_triangle_intersect: non-intersecting", "[igl]")
+{
+  Eigen::Matrix<double,Eigen::Dynamic,3> V(6,3);
+  V<< 
+    0.39234799146652199,0.38443601131439198,0.44744500517845198,
+    0.38924700021743752,0.385637506842613,0.45762349665164948,
+    0.38700349628925301,0.38789276033639897,0.45634675025939975,
+    0.39079749584197976,0.38503675907850249,0.45253425091505073,
+    0.38769650459289529,0.38623825460672351,0.46271274238824822,
+    0.39279749989509577,0.37299175560474374,0.45553924888372399;
+  Eigen::MatrixXi F(2,3);
+  F<< 
+    0,3,2,
+    1,5,4;
+
+  bool coplanar;
+  REQUIRE(!igl::predicates::triangle_triangle_intersect(
+        V.row(F(0,0)).eval(),
+        V.row(F(0,1)).eval(),
+        V.row(F(0,2)).eval(),
+        V.row(F(1,0)).eval(),
+        V.row(F(1,1)).eval(),
+        V.row(F(1,2)).eval(),
+        coplanar));
+}
+
+
+TEST_CASE("triangle_triangle_intersect: coplanar", "[igl]")
+{
+  Eigen::Matrix<double,Eigen::Dynamic,3> V(6,3);
+  V<< 
+  0.30947001278400399,0.80785250663757346,0.47595100104808802,
+  0.299046009778976,0.79801350831985496,0.47843150794506101,
+  0.32418551295995701,0.79794725775718722,0.48203899711370451,
+  0.30425801128148999,0.80293300747871421,0.47719125449657451,
+  0.31897351145744302,0.79302775859832797,0.48327925056219101,
+  0.31376150995492902,0.78810825943946872,0.4845195040106775;
+  Eigen::MatrixXi F(2,3);
+  F<< 
+    0,4,2,
+    1,5,3;
+
+  bool coplanar;
+  REQUIRE(!igl::predicates::triangle_triangle_intersect(
+        V.row(F(0,0)).eval(),
+        V.row(F(0,1)).eval(),
+        V.row(F(0,2)).eval(),
+        V.row(F(1,0)).eval(),
+        V.row(F(1,1)).eval(),
+        V.row(F(1,2)).eval(),coplanar));
+}

tests/test_common.h

@@ -6,6 +6,7 @@
 #include <igl/readDMAT.h>
 
 #include <igl/find.h>
+#include <igl/sortrows.h>
 
 #include <Eigen/Core>
 #include <catch2/catch.hpp>
@@ -204,6 +205,7 @@ namespace test_common
     const Eigen::MatrixBase<DerivedB> & B,
     const EpsType & eps)
   {
+    if(eps == 0){ return assert_eq(A,B); }
     // Sizes should match
     REQUIRE(A.rows() == B.rows());
     REQUIRE(A.cols() == B.cols());
@@ -232,4 +234,19 @@ namespace test_common
     }
   }
 
+  // assert_near after sortrows on each input
+  template <typename DerivedA, typename DerivedB, typename EpsType>
+  void assert_near_rows(
+    const Eigen::MatrixBase<DerivedA> & A,
+    const Eigen::MatrixBase<DerivedB> & B,
+    const EpsType & eps)
+  {
+    Eigen::VectorXi _;
+    DerivedA sA;
+    DerivedB sB;
+    igl::sortrows(A,false,sA,_);
+    igl::sortrows(B,false,sB,_);
+    return assert_near(sA,sB,eps);
+  }
+
 }

tutorial/903_FastFindSelfIntersections/main.cpp

@@ -1,7 +1,7 @@
 #include <igl/read_triangle_mesh.h>
 #include <igl/unique.h>
 #include <igl/opengl/glfw/Viewer.h>
-#include <igl/fast_find_self_intersections.h>
+#include <igl/predicates/find_self_intersections.h>
 
 int main(int argc, char *argv[])
 {
@@ -18,12 +18,13 @@ int main(int argc, char *argv[])
   Eigen::VectorXi EI;
   Eigen::MatrixXd EV;
   Eigen::MatrixXi IF,EE;
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
 
-  if(igl::fast_find_self_intersections(V,F,false,false,IF,EV,EE,EI))
+  if(igl::predicates::find_self_intersections(V,F,IF,CP,EV,EE,EI))
   {
     std::cout<<"Found "<<IF.rows()<<" self intersecting pairs"<<std::endl;
     // plot edge vertices
-    viewer.data().set_edges(EV,EE, Eigen::RowVector3d(1,0,0));
+    viewer.data().set_edges(EV,EE, Eigen::RowVector3d(1,1,1));
   }
   Eigen::VectorXi I;
   igl::unique(IF,I);

tutorial/904_FastFindIntersections/main.cpp

@@ -4,10 +4,11 @@
 
 #include <igl/opengl/glfw/Viewer.h>
 
-#include <igl/fast_find_intersections.h>
+#include <igl/predicates/find_intersections.h>
 
 Eigen::MatrixXd V1,V2;
 Eigen::MatrixXi F1,F2;
+int mid_1, mid_2;
 
 igl::AABB<Eigen::MatrixXd,3> tree;
 
@@ -20,21 +21,29 @@ void update_visualization(igl::opengl::glfw::Viewer & viewer)
   //shifted intersection object
   Eigen::MatrixXd V2_(V2.rows(),V2.cols());
   V2_<< V2.col(0), V2.col(1), V2.col(2).array()+slice_z;
+  viewer.data_list[mid_2].set_vertices(V2_);
 
   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);
+  Eigen::Array<bool,Eigen::Dynamic,1> CP;
+  igl::predicates::find_intersections(tree, V1,F1, V2_,F2,IF,CP,EV,EE,EI);
  
   // Plot the edges of the intersects
-  viewer.data().set_edges( EV,EE, Eigen::RowVector3d(1,0,0));
+  viewer.data_list[mid_1].set_edges( EV,EE, Eigen::RowVector3d(1,1,1));
   
   // show faces which are intersected
   Eigen::VectorXi I;
-  igl::unique(IF,I);
+  igl::unique(IF.col(0).eval(),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);
+  viewer.data_list[mid_1].set_data(D,0,1,igl::COLOR_MAP_TYPE_PARULA);
+
+  igl::unique(IF.col(1).eval(),I);
+  D = Eigen::MatrixXd::Zero(F2.rows(),1);
+  D(I).setConstant(1.0);
+  viewer.data_list[mid_2].set_data(D,0,1,igl::COLOR_MAP_TYPE_JET);
+
 }
 
 
@@ -71,17 +80,28 @@ int main(int argc, char *argv[])
   std::cout<<"'.'/','  push back/pull forward slicing plane."<<std::endl;
   std::cout<<std::endl;
 
-  min_z=V1.col(2).minCoeff();
-  max_z=V1.col(2).maxCoeff();
-  slice_z=(max_z+min_z)/2;
+  if(argc<=2)
+  {
+    min_z=V1.col(2).minCoeff();
+    max_z=V1.col(2).maxCoeff();
+    slice_z=(max_z+min_z)/2;
+  }else
+  {
+    slice_z = 0;
+  }
 
   //center slicing object
   V2.col(2).array() -= V2.col(2).array().mean();
 
   // Plot the mesh
   igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V2, F2);
+  viewer.data().set_face_based(true);
+  mid_2 = viewer.selected_data_index;
+  viewer.append_mesh();
   viewer.data().set_mesh(V1, F1);
   viewer.data().set_face_based(true);
+  mid_1 = viewer.selected_data_index;
   
 
   update_visualization(viewer);

tutorial/908_IntersectionBlockingDecimation/main.cpp

@@ -1,6 +1,6 @@
 #include <igl/opengl/glfw/Viewer.h>
 #include <igl/read_triangle_mesh.h>
-#include <igl/fast_find_self_intersections.h>
+#include <igl/predicates/find_self_intersections.h>
 #include <igl/unique.h>
 #include <igl/remove_unreferenced.h>
 #include <igl/get_seconds.h>
@@ -99,10 +99,9 @@ int main(int argc, char *argv[])
     {
       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);
+        Eigen::MatrixXi IF;
+        Eigen::Array<bool,Eigen::Dynamic,1> CP;
+        igl::predicates::find_self_intersections(dV[pass],dF[pass],false,IF,CP);
         igl::unique(IF,BI);
       }
       printf("  # self-intersections: %d\n",(int)BI.size());

tutorial/CMakeLists.txt

@@ -122,10 +122,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::glfw)
-    igl_add_tutorial(904_FastFindIntersections igl::glfw)
+    igl_add_tutorial(903_FastFindSelfIntersections igl::glfw igl::predicates )
+    igl_add_tutorial(904_FastFindIntersections igl::glfw igl::predicates)
     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)
+    igl_add_tutorial(908_IntersectionBlockingDecimation igl::glfw igl::predicates)
 endif()