Fast mesh-to-mesh intersection and mesh self intersection without CGAL (#2109)

* Added code for fast triangle-triangle intersection checking and function for fast detection of mesh self-intersections and mesh-to-mesh intersections withoug CGAL

* Fixed auto parameters in lambda helper function

* Replaced cbegin/cend with begin/end

* Fixed auto parameters in lambda helper function

* Added tests for igl::tri_tri_intersection_test_3d

* Added more tests, converted macros in Guigue2003_tri_tri_intersect.cpp to proper c++

* Renamed files and function names to follow IGL guidelines,
added reference to the original license for tri_tri_intersect

include/igl/fast_find_intersections.cpp

@@ -0,0 +1,125 @@
+// 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"
+
+
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedI,
+  typename DerivedE>
+IGL_INLINE void igl::fast_find_intersections(
+  const Eigen::MatrixBase<DerivedV1>& V1,
+  const Eigen::MatrixBase<DerivedF1>& F1,
+  const Eigen::MatrixBase<DerivedV2>& V2,
+  const Eigen::MatrixBase<DerivedF2>& F2,
+        Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
+        Eigen::PlainObjectBase<DerivedE>& edges )
+{
+  using AABBTree=igl::AABB<DerivedV1,3>;
+  AABBTree tree;
+  tree.init(V1,F1);
+
+  fast_find_intersections(tree, V1, F1, V2, F2, intersect_pairs, edges);
+}
+
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedI,
+  typename DerivedE>
+IGL_INLINE void igl::fast_find_intersections(
+  const AABB<DerivedV1,3>           & tree,
+  const Eigen::MatrixBase<DerivedV1>& V1,
+  const Eigen::MatrixBase<DerivedF1>& F1,
+  const Eigen::MatrixBase<DerivedV2>& V2,
+  const Eigen::MatrixBase<DerivedF2>& F2,
+        Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
+        Eigen::PlainObjectBase<DerivedE>& edges )
+{
+    using AABBTree=igl::AABB<DerivedV1,3>;
+    using Scalar=typename DerivedV1::Scalar;
+    using BBOX=Eigen::AlignedBox<Scalar,3>;
+    using VERTEX=Eigen::Matrix<typename DerivedE::Scalar,1,3,Eigen::RowMajor>;
+
+    std::vector<VERTEX> _edges;
+    std::vector<int>    _intersect_pairs;
+
+    for(int i=0; i<F2.rows(); ++i)
+    {
+      BBOX tri_box;
+
+      for(int j=0;j<3;++j)
+        tri_box.extend( V2.row( F2(i,j) ).transpose() );
+      
+      // find leaf nodes containing intersecting tri_box
+      // need to specify exact type instead of auto to allow recursion
+      std::function<void(const AABBTree &,int)> check_intersect = 
+        [&](const AABBTree &t,int d) -> void
+      {
+        if(t.m_primitive != -1) //check for the actual intersection //t.is_leaf()
+        {
+          bool coplanar=false;
+          VERTEX edge1,edge2;
+
+          if(igl::tri_tri_intersection_test_3d(
+            V2.row(F2(i,0)),            V2.row(F2(i,1)),            V2.row(F2(i,2)),
+            V1.row(F1(t.m_primitive,0)),V1.row(F1(t.m_primitive,1)),V1.row(F1(t.m_primitive,2)),
+            coplanar,
+            edge1,edge2))
+          {
+            if(!coplanar)
+            {
+              _intersect_pairs.push_back(t.m_primitive);
+              _intersect_pairs.push_back(i);
+
+              _edges.push_back(edge1);
+              _edges.push_back(edge2);
+            }
+          }
+        } else {
+          if(t.m_box.intersects( tri_box )) 
+          { // need to check all branches 
+            check_intersect(*t.m_left, d+1);
+            check_intersect(*t.m_right,d+1);
+          }
+        }
+      };
+
+      // run actual search
+      check_intersect(tree, 0);
+    }
+    edges.resize(_edges.size(), 3);
+
+    for(int i=0; i!=_edges.size(); ++i)
+    {
+      edges.row(i) = _edges[i];
+    }
+
+    intersect_pairs.resize(_intersect_pairs.size()/2,2);
+    for(int i=0; i!=_intersect_pairs.size(); ++i)
+    {
+      intersect_pairs(i/2, i%2) = _intersect_pairs[i];
+    }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::fast_find_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::fast_find_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+
+#endif

include/igl/fast_find_intersections.h

@@ -0,0 +1,78 @@
+// 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_MESH_INTERSECT_H
+#define FAST_FIND_MESH_INTERSECT_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include "AABB.h"
+
+namespace igl {
+
+
+  // identify triangles where two meshes interesect 
+  // using AABBTree and tri_tri_intersection_test_3d
+  // Inputs:
+  //   V1  #V by 3 list representing vertices on the first mesh
+  //   F1  #F by 3 list representing triangles on the first mesh
+  //   V2  #V by 3 list representing vertices on the second mesh
+  //   F2  #F by 3 list representing triangles on the second mesh
+  // Output:
+  //   intersect_pairs  correspondance list of intersecting triangles
+  //                    column 0 - mesh 1, column 1 - mesh2  
+  //   edges      list of pairs of intersection edges
+  template <
+    typename DerivedV1,
+    typename DerivedF1,
+    typename DerivedV2,
+    typename DerivedF2,
+    typename DerivedI,
+    typename DerivedE>
+  IGL_INLINE void fast_find_intersections(
+    const Eigen::MatrixBase<DerivedV1>& V1,
+    const Eigen::MatrixBase<DerivedF1>& F1,
+    const Eigen::MatrixBase<DerivedV2>& V2,
+    const Eigen::MatrixBase<DerivedF2>& F2,
+          Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
+          Eigen::PlainObjectBase<DerivedE>& edges );
+
+  // identify triangles where two meshes interesect 
+  // using AABBTree and tri_tri_intersection_test_3d
+  // Inputs:
+  //   tree - AABB tree bult from the first mesh
+  //   V1  #V by 3 list representing vertices on the first mesh
+  //   F1  #F by 3 list representing triangles on the first mesh
+  //   V2  #V by 3 list representing vertices on the second mesh
+  //   F2  #F by 3 list representing triangles on the second mesh
+  // Output:
+  //   intersect_pairs  correspondance list of intersecting triangles
+  //                    column 0 - mesh 1, column 1 - mesh2  
+  //   edges      list of pairs of intersection edges
+  template <
+    typename DerivedV1,
+    typename DerivedF1,
+    typename DerivedV2,
+    typename DerivedF2,
+    typename DerivedI,
+    typename DerivedE>
+  IGL_INLINE void fast_find_intersections(
+    const AABB<DerivedV1,3>           & tree,
+    const Eigen::MatrixBase<DerivedV1>& V1,
+    const Eigen::MatrixBase<DerivedF1>& F1,
+    const Eigen::MatrixBase<DerivedV2>& V2,
+    const Eigen::MatrixBase<DerivedF2>& F2,
+          Eigen::PlainObjectBase<DerivedI>& intersect_pairs,
+          Eigen::PlainObjectBase<DerivedE>& edges );
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "fast_find_intersections.cpp"
+#endif
+
+#endif

include/igl/fast_find_self_intersections.cpp

@@ -0,0 +1,202 @@
+// 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_self_intersections.h"
+#include "AABB.h"
+#include "tri_tri_intersect.h"
+
+namespace igl{
+namespace internal {
+
+    // helper function, to check if two faces have shared vertices
+    template <typename Derived>
+    bool adjacent_faces(const Eigen::MatrixBase<Derived>& A,
+                        const Eigen::MatrixBase<Derived>& B)
+    {
+        for(int i=0;i<3;++i)
+          for(int j=0;j<3;j++)
+          {
+            if(A(i)==B(j))
+              return true;
+          }
+        return false;
+    }
+
+}
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI>
+IGL_INLINE bool igl::fast_find_self_intersections(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+        Eigen::PlainObjectBase<DerivedI>& intersect)
+{
+    using Scalar=typename DerivedV::Scalar;
+    using BBOX=Eigen::AlignedBox<Scalar,3>;
+    using AABBTree=igl::AABB<DerivedV,3>;
+    AABBTree tree;
+
+    tree.init(V,F);
+    bool _intersects=false;
+
+    intersect.resize(F.rows(),1);
+    intersect.setConstant(0);
+
+    for(int i=0; i<F.rows(); ++i)
+    {
+      if( intersect(i) ) continue;
+
+      BBOX tri_box;
+
+      for(int j=0;j<3;++j)
+        tri_box.extend( V.row( F(i,j) ).transpose() );
+      
+      // find leaf nodes containing intersecting tri_box
+      // need to declare full type to enable recursion
+      std::function<bool(const AABBTree &,int)> check_intersect = 
+        [&](const AABBTree &t,int d) -> bool
+      {
+        if(t.m_primitive != -1) //check for the actual intersection (is_leaf)
+        {
+          if(t.m_primitive==i) //itself
+              return false;
+          if(igl::internal::adjacent_faces(F.row(i), F.row(t.m_primitive)) )
+              return false;
+
+          bool coplanar=false;
+          Eigen::Matrix<Scalar,1,3,Eigen::RowMajor> edge1,edge2;
+
+          if(igl::tri_tri_intersection_test_3d(
+            V.row(F(i,0)),V.row(F(i,1)),V.row(F(i,2)),
+            V.row(F(t.m_primitive,0)),V.row(F(t.m_primitive,1)),V.row(F(t.m_primitive,2)),
+            coplanar,
+            edge1,edge2))
+          {
+            if(!coplanar)
+            {
+              intersect(i)=1;
+              intersect(t.m_primitive)=1;
+              return true;
+            }
+          }
+        } else {
+          if(t.m_box.intersects(tri_box)) {
+            // need to check both subtrees
+            bool r1=check_intersect(*t.m_left ,d+1);
+            bool r2=check_intersect(*t.m_right,d+1);
+            return r1||r2;
+          }
+        }
+        return false;
+      };
+
+      bool r=check_intersect(tree,0);
+      _intersects = _intersects || r;
+    }
+    return _intersects;
+}
+
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedI,
+  typename DerivedE>
+IGL_INLINE bool igl::fast_find_self_intersections(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const Eigen::MatrixBase<DerivedF>& F,
+        Eigen::PlainObjectBase<DerivedI>& intersect,
+        Eigen::PlainObjectBase<DerivedE>& edges )
+{
+    using Scalar=typename DerivedV::Scalar;
+    using BBOX=Eigen::AlignedBox<Scalar,3>;
+    using AABBTree=igl::AABB<DerivedV,3>;
+    using EDGE=Eigen::Matrix<typename DerivedE::Scalar,1,3,Eigen::RowMajor>;
+
+    std::vector<EDGE> _edges;
+    AABBTree tree;
+
+    tree.init(V,F);
+    bool _intersects=false;
+
+    intersect.resize(F.rows(),1);
+    intersect.setConstant(0);
+
+    for(int i=0; i<F.rows(); ++i)
+    {
+      int _inter=-1;
+      if( intersect(i) ) continue;
+
+      BBOX tri_box;
+
+      for(int j=0;j<3;++j)
+        tri_box.extend( V.row( F(i,j) ).transpose() );
+      
+      // find leaf nodes containing intersecting tri_box
+      std::function<bool(const AABBTree &,int)> check_intersect = 
+        [&](const AABBTree &t,int d) -> bool
+      {
+        if(t.m_primitive != -1) //check for the actual intersection //t.is_leaf()
+        {
+          if(t.m_primitive==i) //itself
+              return false;
+          if(igl::internal::adjacent_faces(F.row(i), F.row(t.m_primitive)) )
+              return false;
+
+          bool coplanar=false;
+          EDGE edge1,edge2;
+
+          if(igl::tri_tri_intersection_test_3d(
+            V.row(F(i,0)),            V.row(F(i,1)),            V.row(F(i,2)),
+            V.row(F(t.m_primitive,0)),V.row(F(t.m_primitive,1)),V.row(F(t.m_primitive,2)),
+            coplanar,
+            edge1,edge2))
+          {
+            if(!coplanar)
+            {
+              intersect(i)=1;
+              intersect(t.m_primitive)=1;
+              _edges.push_back(edge1);
+              _edges.push_back(edge2);
+              return true;
+            }
+          }
+        } else {
+          if(t.m_box.intersects( tri_box ))
+          {
+            bool r1=check_intersect(*t.m_left, d+1);
+            bool r2=check_intersect(*t.m_right,d+1);
+            return r1||r2;
+          }
+        }
+        return false;
+      };
+
+      // run actual search
+      bool r=check_intersect(tree, 0);
+      _intersects = _intersects || r;
+    }
+
+    edges.resize(_edges.size(),3);
+    int i=0;
+
+    for(auto e=std::begin(_edges);e!=std::end(_edges);++e,++i)
+    {
+      edges.row(i) = *e;
+    }
+    return _intersects;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template bool igl::fast_find_self_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template bool igl::fast_find_self_intersections<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/fast_find_self_intersections.h

@@ -0,0 +1,62 @@
+// 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 mesh intersects itself
+  // using AABBTree and tri_tri_intersection_test_3d
+  // Inputs:
+  //   V  #V by 3 list representing vertices
+  //   F  #F by 3 list representing triangles.
+  // Output:
+  //   intersect  #F by 1 indicator that triangle intersects anothe triangle
+  // Returns:
+  //   self-interection found ?
+    template <
+      typename DerivedV,
+      typename DerivedF,
+      typename DerivedI>
+    IGL_INLINE bool fast_find_self_intersections(
+      const Eigen::MatrixBase<DerivedV>& V,
+      const Eigen::MatrixBase<DerivedF>& F,
+      Eigen::PlainObjectBase<DerivedI>& intersect);
+
+  // identify triangles where mesh intersects itself
+  // using AABBTree and tri_tri_intersection_test_3d
+  // Inputs:
+  //   V  #V by 3 list representing vertices
+  //   F  #F by 3 list representing triangles.
+  // Output:
+  //   intersect  #F by 1 indicator that triangle intersects anothe triangle
+  //   edges      list of pairs of intersection edges
+  // Returns:
+  //   self-interection found ?
+    template <
+      typename DerivedV,
+      typename DerivedF,
+      typename DerivedI,
+      typename DerivedE>
+    IGL_INLINE bool fast_find_self_intersections(
+      const Eigen::MatrixBase<DerivedV>& V,
+      const Eigen::MatrixBase<DerivedF>& F,
+      Eigen::PlainObjectBase<DerivedI>& intersect,
+      Eigen::PlainObjectBase<DerivedE>& edges );
+
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "fast_find_self_intersections.cpp"
+#endif
+
+#endif

include/igl/tri_tri_intersect.cpp

@@ -0,0 +1,894 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2021 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/.
+
+/*
+*  
+*  C++ version based on the routines published in    
+*  "Fast and Robust Triangle-Triangle Overlap Test      
+*   Using Orientation Predicates"  P. Guigue - O. Devillers
+*  
+*  Works with Eigen data structures instead of plain C arrays
+*  returns bool values
+* 
+*  Code is rewritten to get rid of the macros and use C++ lambda and 
+*  inline functions instead
+*  
+*  Original notice: 
+*
+*  Triangle-Triangle Overlap Test Routines        
+*  July, 2002                                                          
+*  Updated December 2003
+*
+*  Updated by Vladimir S. FONOV 
+*  March, 2023                                             
+*                                                                       
+*  This file contains C implementation of algorithms for                
+*  performing two and three-dimensional triangle-triangle intersection test 
+*  The algorithms and underlying theory are described in                    
+*                                                                           
+* "Fast and Robust Triangle-Triangle Overlap Test 
+*  Using Orientation Predicates"  P. Guigue - O. Devillers
+*                                                 
+*  Journal of Graphics Tools, 8(1), 2003                                    
+*                                                                           
+*  Several geometric predicates are defined.  Their parameters are all      
+*  points.  Each point is an array of two or three double precision         
+*  floating point numbers. The geometric predicates implemented in          
+*  this file are:                                                            
+*                                                                           
+*    int tri_tri_overlap_test_3d(p1,q1,r1,p2,q2,r2)                         
+*    int tri_tri_overlap_test_2d(p1,q1,r1,p2,q2,r2)                         
+*                                                                           
+*    int tri_tri_intersection_test_3d(p1,q1,r1,p2,q2,r2,
+*                                     coplanar,source,target)               
+*                                                                           
+*       is a version that computes the segment of intersection when            
+*       the triangles overlap (and are not coplanar)                        
+*                                                                           
+*    each function returns 1 if the triangles (including their              
+*    boundary) intersect, otherwise 0                                       
+*                                                                           
+*                                                                           
+*  Other information are available from the Web page                        
+*  http://www.acm.org/jgt/papers/GuigueDevillers03/                         
+*                                                                           
+*/
+
+#ifndef IGL_TRI_TRI_INTERSECT_CPP
+#define IGL_TRI_TRI_INTERSECT_CPP
+
+#include "tri_tri_intersect.h"
+#include <Eigen/Geometry>
+
+// helper functions
+namespace igl {
+
+namespace internal {
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+typename DerivedN1>
+IGL_INLINE bool coplanar_tri_tri3d(
+  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,
+  const Eigen::MatrixBase<DerivedN1> &normal_1);
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2>
+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);
+    v1=q2-q1;
+
+    if (v1.dot(N1) > 0.0) return false;
+
+    v1=p2-p1;
+    v2=r1-p1;
+    N1=v1.cross(v2);
+    v1=r2-p1;
+
+    if (v1.dot(N1) > 0.0) return false;
+    else return true;
+  }
+
+
+
+/* Permutation in a canonical form of T2's vertices */
+
+  template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+            typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+            typename DP2,typename DQ2,typename DR2,
+            typename DerivedN1>
+  inline bool _IGL_TRI_TRI_3D(
+    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)
+  {
+    if (dp2 > 0.0) { 
+      if (dq2 > 0.0) return _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2);
+      else if (dr2 > 0.0) return _IGL_CHECK_MIN_MAX(p1,r1,q1,q2,r2,p2);
+      else return _IGL_CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2); }
+    else if (dp2 < 0.0) { 
+      if (dq2 < 0.0) return _IGL_CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2);
+      else if (dr2 < 0.0) return _IGL_CHECK_MIN_MAX(p1,q1,r1,q2,r2,p2);
+      else return _IGL_CHECK_MIN_MAX(p1,r1,q1,p2,q2,r2);
+    } else { 
+      if (dq2 < 0.0) { 
+        if (dr2 >= 0.0)  return  _IGL_CHECK_MIN_MAX(p1,r1,q1,q2,r2,p2);
+        else return _IGL_CHECK_MIN_MAX(p1,q1,r1,p2,q2,r2);
+      }
+      else if (dq2 > 0.0) { 
+        if (dr2 > 0.0) return _IGL_CHECK_MIN_MAX(p1,r1,q1,p2,q2,r2);
+        else  return _IGL_CHECK_MIN_MAX(p1,q1,r1,q2,r2,p2);
+      } 
+      else  { 
+        if (dr2 > 0.0) return _IGL_CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2);
+        else if (dr2 < 0.0) return _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2);
+        else return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
+      }}
+  }
+
+
+// #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
+/*
+*
+*  Three-dimensional Triangle-Triangle Overlap Test
+*
+*/
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2> 
+IGL_INLINE bool igl::tri_tri_overlap_test_3d(
+  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)
+{
+  using Scalar    = typename DerivedP1::Scalar;
+  //using RowVector = typename Eigen::Matrix<Scalar, 1, 3>;
+ 
+  Scalar dp1, dq1, dr1, dp2, dq2, dr2;
+  //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=p1-r2;
+  dp1 = v1.dot(N2);
+  v1=q1-r2;
+  dq1 = v1.dot(N2);
+  v1=r1-r2;
+  dr1 = v1.dot(N2);
+  
+  if (((dp1 * dq1) > 0.0) && ((dp1 * dr1) > 0.0))  return false; 
+
+  /* Compute distance signs  of p2, q2 and r2 to the plane of
+     triangle(p1,q1,r1) */
+
+  v1=q1-p1;
+  v2=r1-p1;
+  auto N1=v1.cross(v2);
+
+  v1=p2-r1;
+  dp2 = v1.dot(N1);
+  v1=q2-r1;
+  dq2 = v1.dot(N1);
+  v1=r2-r1;
+  dr2 = v1.dot(N1);
+  
+  if (((dp2 * dq2) > 0.0) && ((dp2 * dr2) > 0.0)) return false;
+
+  /* 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);
+  } 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);
+  } 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);
+    }
+    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);
+    }
+    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);
+    }
+  }
+};
+
+
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+typename DerivedN1>
+IGL_INLINE bool igl::internal::coplanar_tri_tri3d(
+  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,
+  const Eigen::MatrixBase<DerivedN1> &normal_1)
+{
+
+  using Scalar= typename DerivedP1::Scalar;
+  using RowVector2D = typename Eigen::Matrix<Scalar,1,2>;
+ 
+  RowVector2D P1,Q1,R1;
+  RowVector2D P2,Q2,R2;
+
+  Scalar n_x, n_y, n_z;
+
+  n_x = ((normal_1[0]<0.0)?-normal_1[0]:normal_1[0]);
+  n_y = ((normal_1[1]<0.0)?-normal_1[1]:normal_1[1]);
+  n_z = ((normal_1[2]<0.0)?-normal_1[2]:normal_1[2]);
+
+
+  /* Projection of the triangles in 3D onto 2D such that the area of
+     the projection is maximized. */
+
+
+  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);
+    
+};
+
+
+namespace igl 
+{
+  namespace internal {
+/*
+*                                                                
+*  Three-dimensional Triangle-Triangle Intersection              
+*
+*/
+
+/*
+   This macro is called when the triangles surely intersect
+   It constructs the segment of intersection of the two triangles
+   if they are not coplanar.
+*/
+
+// NOTE: a faster, but possibly less precise, method of computing
+// point B is described here: https://github.com/erich666/jgt-code/issues/5
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+typename DerivedS,typename DerivedT,
+typename DerivedN1,typename DerivedN2>
+bool _IGL_CONSTRUCT_INTERSECTION(
+  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,
+        Eigen::MatrixBase<DerivedS>  &source, Eigen::MatrixBase<DerivedT> &target,
+  const Eigen::MatrixBase<DerivedN1> &N1,const Eigen::MatrixBase<DerivedN2> &N2)
+{
+  using Scalar = typename DerivedP1::Scalar;
+  using RowVector3D = typename Eigen::Matrix<Scalar,1,3>;
+
+  RowVector3D v,v1,v2,N;
+
+  v1=q1-p1;
+  v2=r2-p1;
+  N=v1.cross(v2);
+  v=p2-p1;
+  if (v.dot(N) > 0.0) {
+    v1=r1-p1;
+    N=v1.cross(v2);
+    if (v.dot(N) <= 0.0) {
+      v2=q2-p1;
+      N=v1.cross(v2);
+      if (v.dot(N) > 0.0) {
+        v1=p1-p2;
+        v2=p1-r1;
+        Scalar alpha = v1.dot(N2) / v2.dot(N2);
+        v1=v2*alpha;
+        source=p1-v1;
+        v1=p2-p1;
+        v2=p2-r2;
+        alpha = v1.dot(N1) / v2.dot(N1);
+        v1=v2*alpha;
+        target=p2-v1;
+        return true;
+      } else { 
+        v1=p2-p1;
+        v2=p2-q2;
+        Scalar alpha = v1.dot(N1) / v2.dot(N1);
+        v1=v2*alpha;
+        source=p2-v1;
+        v1=p2-p1;
+        v2=p2-r2;
+        alpha = v1.dot(N1) / v2.dot(N1);
+        v1=v2*alpha;
+        target=p2-v1;
+        return true;
+      } 
+    } else {
+      return false;
+    } 
+  } else { 
+    v2=q2-p1;
+    N=v1.cross(v2);
+    if (v.dot(N) < 0.0) {
+      return false;
+    } else {
+      v1=r1-p1;
+      N=v1.cross(v2);
+      if (v.dot(N) >= 0.0) { 
+        v1=p1-p2;
+        v2=p1-r1;
+        Scalar alpha = v1.dot(N2) / v2.dot(N2);
+        v1=v2*alpha;
+        source=p1-v1;
+        v1=p1-p2;
+        v2=p1-q1;
+        alpha = v1.dot(N2) / v2.dot(N2);
+        v1=v2*alpha;
+        target=p1-v1 ;
+        return true; 
+      } else { 
+        v1=p2-p1 ;
+        v2=p2-q2 ;
+        Scalar alpha = v1.dot(N1) / v2.dot(N1);
+        v1=v2*alpha;
+        source=p2-v1;
+        v1=p1-p2;
+        v2=p1-q1;
+        alpha = v1.dot(N2) / v2.dot(N2);
+        v1=v2*alpha;
+        target=p1-v1 ;
+        return true;
+    }}}
+}
+
+
+
+// #define _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2) { \
+//   _IGL_SUB(v1,q1,p1) \
+//   _IGL_SUB(v2,r2,p1) \
+//   _IGL_CROSS(N,v1,v2) \
+//   _IGL_SUB(v,p2,p1) \
+//   if (_IGL_DOT(v,N) > 0.0) {\
+//     _IGL_SUB(v1,r1,p1) \
+//     _IGL_CROSS(N,v1,v2) \
+//     if (_IGL_DOT(v,N) <= 0.0) { \
+//       _IGL_SUB(v2,q2,p1) \
+//       _IGL_CROSS(N,v1,v2) \
+//       if (_IGL_DOT(v,N) > 0.0) { \
+//   _IGL_SUB(v1,p1,p2) \
+//   _IGL_SUB(v2,p1,r1) \
+//   alpha = _IGL_DOT(v1,N2) / _IGL_DOT(v2,N2); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(source,p1,v1) \
+//   _IGL_SUB(v1,p2,p1) \
+//   _IGL_SUB(v2,p2,r2) \
+//   alpha = _IGL_DOT(v1,N1) / _IGL_DOT(v2,N1); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(target,p2,v1) \
+//   return true; \
+//       } else { \
+//   _IGL_SUB(v1,p2,p1) \
+//   _IGL_SUB(v2,p2,q2) \
+//   alpha = _IGL_DOT(v1,N1) / _IGL_DOT(v2,N1); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(source,p2,v1) \
+//   _IGL_SUB(v1,p2,p1) \
+//   _IGL_SUB(v2,p2,r2) \
+//   alpha = _IGL_DOT(v1,N1) / _IGL_DOT(v2,N1); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(target,p2,v1) \
+//   return true; \
+//       } \
+//     } else { \
+//       return false; \
+//     } \
+//   } else { \
+//     _IGL_SUB(v2,q2,p1) \
+//     _IGL_CROSS(N,v1,v2) \
+//     if (_IGL_DOT(v,N) < 0.0) { \
+//       return false; \
+//     } else { \
+//       _IGL_SUB(v1,r1,p1) \
+//       _IGL_CROSS(N,v1,v2) \
+//       if (_IGL_DOT(v,N) >= 0.0) { \
+//   _IGL_SUB(v1,p1,p2) \
+//   _IGL_SUB(v2,p1,r1) \
+//   alpha = _IGL_DOT(v1,N2) / _IGL_DOT(v2,N2); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(source,p1,v1) \
+//   _IGL_SUB(v1,p1,p2) \
+//   _IGL_SUB(v2,p1,q1) \
+//   alpha = _IGL_DOT(v1,N2) / _IGL_DOT(v2,N2); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(target,p1,v1) \
+//   return true; \
+//       } else { \
+//   _IGL_SUB(v1,p2,p1) \
+//   _IGL_SUB(v2,p2,q2) \
+//   alpha = _IGL_DOT(v1,N1) / _IGL_DOT(v2,N1); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(source,p2,v1) \
+//   _IGL_SUB(v1,p1,p2) \
+//   _IGL_SUB(v2,p1,q1) \
+//   alpha = _IGL_DOT(v1,N2) / _IGL_DOT(v2,N2); \
+//   _IGL_SCALAR(v1,alpha,v2) \
+//   _IGL_SUB(target,p1,v1) \
+//   return true; \
+//       }}}}
+
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+            typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+            typename DP2,typename DQ2,typename DR2,
+            typename DerivedS,typename DerivedT,
+            typename DerivedN1,typename DerivedN2
+            >
+  inline bool _IGL_TRI_TRI_INTER_3D(
+    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,
+    bool & coplanar,
+    Eigen::MatrixBase<DerivedS> &source, Eigen::MatrixBase<DerivedT> &target,
+    const Eigen::MatrixBase<DerivedN1> &N1,const Eigen::MatrixBase<DerivedN2> &N2
+    )
+  {
+    if (dp2 > 0.0) { 
+     if (dq2 > 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,r2,p2,q2,source,target,N1,N2);
+     else if (dr2 > 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,q2,r2,p2,source,target,N1,N2);
+          else return _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2,source,target,N1,N2); }
+    else if (dp2 < 0.0) { 
+      if (dq2 < 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,r2,p2,q2,source,target,N1,N2);
+      else if (dr2 < 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,q2,r2,p2,source,target,N1,N2);
+          else return _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,p2,q2,r2,source,target,N1,N2);
+    } else { 
+      if (dq2 < 0.0) { 
+        if (dr2 >= 0.0)  return _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,q2,r2,p2,source,target,N1,N2);
+        else return _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2,source,target,N1,N2);
+    } 
+    else if (dq2 > 0.0) {
+      if (dr2 > 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,p2,q2,r2,source,target,N1,N2);
+      else  return _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,q2,r2,p2,source,target,N1,N2);
+    } 
+    else  {
+      if (dr2 > 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,r2,p2,q2,source,target,N1,N2);
+      else if (dr2 < 0.0) return _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,r2,p2,q2,source,target,N1,N2);
+      else { 
+        coplanar = true;
+        return igl::internal::coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
+     }
+  }}
+
+  }
+
+// #define _IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2) { \
+//   if (dp2 > 0.0) { \
+//      if (dq2 > 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,r2,p2,q2) \
+//      else if (dr2 > 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,q2,r2,p2)\
+//      else _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2) }\
+//   else if (dp2 < 0.0) { \
+//     if (dq2 < 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,r2,p2,q2)\
+//     else if (dr2 < 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,q2,r2,p2)\
+//     else _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,p2,q2,r2)\
+//   } else { \
+//     if (dq2 < 0.0) { \
+//       if (dr2 >= 0.0)  _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,q2,r2,p2)\
+//       else _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2)\
+//     } \
+//     else if (dq2 > 0.0) { \
+//       if (dr2 > 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,p2,q2,r2)\
+//       else  _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,q2,r2,p2)\
+//     } \
+//     else  { \
+//       if (dr2 > 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,r2,p2,q2)\
+//       else if (dr2 < 0.0) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,r2,p2,q2)\
+//       else { \
+//         coplanar = true; \
+//   return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);\
+//      } \
+//   }} }
+  
+  } //internal 
+} //igl
+
+/*
+   The following version computes the segment of intersection of the
+   two triangles if it exists. 
+   coplanar returns whether the triangles are coplanar
+   source and target are the endpoints of the line segment of intersection 
+*/
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+typename DerivedS,typename DerivedT>
+IGL_INLINE bool igl::tri_tri_intersection_test_3d(
+    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,
+    bool & coplanar,
+    Eigen::MatrixBase<DerivedS> & source, Eigen::MatrixBase<DerivedT> & target )
+{
+  using Scalar = typename DerivedP1::Scalar;
+  using RowVector3D = typename Eigen::Matrix<Scalar, 1, 3>;
+
+  Scalar dp1, dq1, dr1, dp2, dq2, dr2;
+  RowVector3D v1, v2, v;
+  RowVector3D N1, N2, N;
+  Scalar alpha;
+  // Compute distance signs  of p1, q1 and r1 
+  // to the plane of triangle(p2,q2,r2)
+
+  v1=p2-r2;
+  v2=q2-r2;
+  N2=v1.cross(v2);
+
+  v1=p1-r2;
+  dp1 = v1.dot(N2);
+  v1=q1-r2;
+  dq1 = v1.dot(N2);
+  v1=r1-r2;
+  dr1 = v1.dot(N2);
+  
+  coplanar = false;
+
+  if (((dp1 * dq1) > 0.0) && ((dp1 * dr1) > 0.0))  return false; 
+
+  // Compute distance signs  of p2, q2 and r2 
+  // to the plane of triangle(p1,q1,r1)
+
+  
+  v1=q1-p1;
+  v2=r1-p1;
+  N1=v1.cross(v2);
+
+  v1=p2-r1;
+  dp2 = v1.dot(N1);
+  v1=q2-r1;
+  dq2 = v1.dot(N1);
+  v1=r2-r1;
+  dr2 = v1.dot(N1);
+  
+  if (((dp2 * dq2) > 0.0) && ((dp2 * dr2) > 0.0)) return false;
+
+  // Permutation in a canonical form of T1's vertices
+  if (dp1 > 0.0) {
+    if (dq1 > 0.0) return internal::_IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2,coplanar,source,target,N1,N2);
+    else if (dr1 > 0.0) return internal::_IGL_TRI_TRI_INTER_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2,coplanar,source,target,N1,N2);
+  
+    else return internal::_IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2,coplanar,source,target,N1,N2);
+  } else if (dp1 < 0.0) {
+    if (dq1 < 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(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2,coplanar,source,target,N1,N2);
+    else return internal::_IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2,coplanar,source,target,N1,N2);
+  } else {
+    if (dq1 < 0.0) {
+      if (dr1 >= 0.0) return internal::_IGL_TRI_TRI_INTER_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2,coplanar,source,target,N1,N2);
+      else return internal::_IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2,coplanar,source,target,N1,N2);
+    }
+    else if (dq1 > 0.0) {
+      if (dr1 > 0.0) return internal::_IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2,coplanar,source,target,N1,N2);
+      else return internal::_IGL_TRI_TRI_INTER_3D(q1,r1,p1,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,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
+
+        coplanar = true;
+        return internal::coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
+      }
+    }
+  }
+};
+
+
+
+namespace igl {
+  namespace internal {
+/*
+*
+*  Two dimensional Triangle-Triangle Overlap Test    
+*
+*/
+
+
+/* some 2D macros */
+
+//#define _IGL_ORIENT_2D(a, b, c)  ((a[0]-c[0])*(b[1]-c[1])-(a[1]-c[1])*(b[0]-c[0]))
+template <typename DerivedA,typename DerivedB,typename DerivedC>
+  inline typename Eigen::MatrixBase<DerivedA>::Scalar _IGL_ORIENT_2D(
+    const Eigen::MatrixBase<DerivedA> &  a, 
+    const Eigen::MatrixBase<DerivedB> &  b,
+    const Eigen::MatrixBase<DerivedC> &  c)
+    {
+      return  ((a[0]-c[0])*(b[1]-c[1])-(a[1]-c[1])*(b[0]-c[0]));
+    }
+
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+          typename DerivedP2,typename DerivedQ2,typename DerivedR2>
+bool _IGL_INTERSECTION_TEST_VERTEX(
+  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
+)
+{
+  if (_IGL_ORIENT_2D(R2,P2,Q1) >= 0.0)
+    if (_IGL_ORIENT_2D(R2,Q2,Q1) <= 0.0)
+      if (_IGL_ORIENT_2D(P1,P2,Q1) > 0.0) {
+        if (_IGL_ORIENT_2D(P1,Q2,Q1) <= 0.0) return true;
+        else return false;} 
+      else {
+        if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0)
+          if (_IGL_ORIENT_2D(Q1,R1,P2) >= 0.0) return true; 
+          else return false;
+        else return false;
+      }
+      else 
+        if (_IGL_ORIENT_2D(P1,Q2,Q1) <= 0.0)
+          if (_IGL_ORIENT_2D(R2,Q2,R1) <= 0.0)
+            if (_IGL_ORIENT_2D(Q1,R1,Q2) >= 0.0) return true; 
+            else return false;
+          else return false;
+        else return false;
+      else
+        if (_IGL_ORIENT_2D(R2,P2,R1) >= 0.0) 
+          if (_IGL_ORIENT_2D(Q1,R1,R2) >= 0.0)
+            if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0) return true;
+            else return false;
+          else 
+            if (_IGL_ORIENT_2D(Q1,R1,Q2) >= 0.0) {
+              if (_IGL_ORIENT_2D(R2,R1,Q2) >= 0.0) return true; 
+              else return false; 
+            }
+        else return false; 
+  else  return false; 
+}
+
+// #define INTERSECTION_TEST_VERTEX(P1, Q1, R1, P2, Q2, R2) {\
+//   if (_IGL_ORIENT_2D(R2,P2,Q1) >= 0.0)\
+//     if (_IGL_ORIENT_2D(R2,Q2,Q1) <= 0.0)\
+//       if (_IGL_ORIENT_2D(P1,P2,Q1) > 0.0) {\
+//   if (_IGL_ORIENT_2D(P1,Q2,Q1) <= 0.0) return true; \
+//   else return false;} else {\
+//   if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0)\
+//     if (_IGL_ORIENT_2D(Q1,R1,P2) >= 0.0) return true; \
+//     else return false;\
+//   else return false;}\
+//     else \
+//       if (_IGL_ORIENT_2D(P1,Q2,Q1) <= 0.0)\
+//   if (_IGL_ORIENT_2D(R2,Q2,R1) <= 0.0)\
+//     if (_IGL_ORIENT_2D(Q1,R1,Q2) >= 0.0) return true; \
+//     else return false;\
+//   else return false;\
+//       else return false;\
+//   else\
+//     if (_IGL_ORIENT_2D(R2,P2,R1) >= 0.0) \
+//       if (_IGL_ORIENT_2D(Q1,R1,R2) >= 0.0)\
+//   if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0) return true;\
+//   else return false;\
+//       else \
+//   if (_IGL_ORIENT_2D(Q1,R1,Q2) >= 0.0) {\
+//     if (_IGL_ORIENT_2D(R2,R1,Q2) >= 0.0) return true; \
+//     else return false; }\
+//   else return false; \
+//     else  return false; \
+//  };
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+          typename DerivedP2,typename DerivedQ2,typename DerivedR2>
+bool _IGL_INTERSECTION_TEST_EDGE(
+  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
+)
+{
+  if (_IGL_ORIENT_2D(R2,P2,Q1) >= 0.0) {
+    if (_IGL_ORIENT_2D(P1,P2,Q1) >= 0.0) {
+        if (_IGL_ORIENT_2D(P1,Q1,R2) >= 0.0) return true;
+        else return false;} else { 
+      if (_IGL_ORIENT_2D(Q1,R1,P2) >= 0.0){ 
+  if (_IGL_ORIENT_2D(R1,P1,P2) >= 0.0) return true; else return false;} 
+      else return false; } 
+  } else {
+    if (_IGL_ORIENT_2D(R2,P2,R1) >= 0.0) {
+      if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0) {
+  if (_IGL_ORIENT_2D(P1,R1,R2) >= 0.0) return true;  
+  else {
+    if (_IGL_ORIENT_2D(Q1,R1,R2) >= 0.0) return true; else return false;}}
+      else  return false; }
+    else return false; }
+}
+
+
+// #define _IGL_INTERSECTION_TEST_EDGE(P1, Q1, R1, P2, Q2, R2) { \
+//   if (_IGL_ORIENT_2D(R2,P2,Q1) >= 0.0) {\
+//     if (_IGL_ORIENT_2D(P1,P2,Q1) >= 0.0) { \
+//         if (_IGL_ORIENT_2D(P1,Q1,R2) >= 0.0) return true; \
+//         else return false;} else { \
+//       if (_IGL_ORIENT_2D(Q1,R1,P2) >= 0.0){ \
+//   if (_IGL_ORIENT_2D(R1,P1,P2) >= 0.0) return true; else return false;} \
+//       else return false; } \
+//   } else {\
+//     if (_IGL_ORIENT_2D(R2,P2,R1) >= 0.0) {\
+//       if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0) {\
+//   if (_IGL_ORIENT_2D(P1,R1,R2) >= 0.0) return true;  \
+//   else {\
+//     if (_IGL_ORIENT_2D(Q1,R1,R2) >= 0.0) return true; else return false;}}\
+//       else  return false; }\
+//     else return false; }}
+
+
+
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2>
+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) 
+{
+  if ( _IGL_ORIENT_2D(p2,q2,p1) >= 0.0 ) {
+    if ( _IGL_ORIENT_2D(q2,r2,p1) >= 0.0 ) {
+      if ( _IGL_ORIENT_2D(r2,p2,p1) >= 0.0 ) return true;
+      else return _IGL_INTERSECTION_TEST_EDGE(p1,q1,r1,p2,q2,r2);
+    } else {  
+      if ( _IGL_ORIENT_2D(r2,p2,p1) >= 0.0 ) 
+      return _IGL_INTERSECTION_TEST_EDGE(p1,q1,r1,r2,p2,q2);
+      else return _IGL_INTERSECTION_TEST_VERTEX(p1,q1,r1,p2,q2,r2);}}
+  else {
+    if ( _IGL_ORIENT_2D(q2,r2,p1) >= 0.0 ) {
+      if ( _IGL_ORIENT_2D(r2,p2,p1) >= 0.0 ) 
+        return _IGL_INTERSECTION_TEST_EDGE(p1,q1,r1,q2,r2,p2);
+      else  return _IGL_INTERSECTION_TEST_VERTEX(p1,q1,r1,q2,r2,p2);}
+    else return _IGL_INTERSECTION_TEST_VERTEX(p1,q1,r1,r2,p2,q2);}
+};
+
+  }//internal
+} //igl
+
+
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2>
+IGL_INLINE bool igl::tri_tri_overlap_test_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) 
+{
+  if ( igl::internal::_IGL_ORIENT_2D(p1,q1,r1) < 0.0)
+    if ( igl::internal::_IGL_ORIENT_2D(p2,q2,r2) < 0.0)
+      return igl::internal::ccw_tri_tri_intersection_2d(p1,r1,q1,p2,r2,q2);
+    else
+      return igl::internal::ccw_tri_tri_intersection_2d(p1,r1,q1,p2,q2,r2);
+  else
+    if ( igl::internal::_IGL_ORIENT_2D(p2,q2,r2) < 0.0 )
+      return igl::internal::ccw_tri_tri_intersection_2d(p1,q1,r1,p2,r2,q2);
+    else
+      return igl::internal::ccw_tri_tri_intersection_2d(p1,q1,r1,p2,q2,r2);
+};
+
+#endif //IGL_TRI_TRI_INTERSECT_CPP
+
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::tri_tri_intersection_test_3d<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1> const, 1, -1, false> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template bool igl::tri_tri_intersection_test_3d<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, bool&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+#endif

include/igl/tri_tri_intersect.h

@@ -0,0 +1,161 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2021 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/
+
+/*
+*  
+*  C++ version based on the routines published in    
+*  "Fast and Robust Triangle-Triangle Overlap Test      
+*   Using Orientation Predicates"  P. Guigue - O. Devillers
+*  
+*  Works with Eigen data structures instead of plain C arrays
+*  returns bool values
+* 
+*  Code is rewritten to get rid of the macros and use C++ lambda and 
+*  inline functions instead
+*  
+*  Original notice: 
+*
+*  Triangle-Triangle Overlap Test Routines        
+*  July, 2002                                                          
+*  Updated December 2003
+*
+*  Updated by Vladimir S. FONOV 
+*  March, 2023                                             
+*                                                                       
+*  This file contains C implementation of algorithms for                
+*  performing two and three-dimensional triangle-triangle intersection test 
+*  The algorithms and underlying theory are described in                    
+*                                                                           
+* "Fast and Robust Triangle-Triangle Overlap Test 
+*  Using Orientation Predicates"  P. Guigue - O. Devillers
+*                                                 
+*  Journal of Graphics Tools, 8(1), 2003                                    
+*                                                                           
+*  Several geometric predicates are defined.  Their parameters are all      
+*  points.  Each point is an array of two or three double precision         
+*  floating point numbers. The geometric predicates implemented in          
+*  this file are:                                                            
+*                                                                           
+*    int tri_tri_overlap_test_3d(p1,q1,r1,p2,q2,r2)                         
+*    int tri_tri_overlap_test_2d(p1,q1,r1,p2,q2,r2)                         
+*                                                                           
+*    int tri_tri_intersection_test_3d(p1,q1,r1,p2,q2,r2,
+*                                     coplanar,source,target)               
+*                                                                           
+*       is a version that computes the segment of intersection when            
+*       the triangles overlap (and are not coplanar)                        
+*                                                                           
+*    each function returns 1 if the triangles (including their              
+*    boundary) intersect, otherwise 0                                       
+*                                                                           
+*                                                                           
+*  Other information are available from the Web page                        
+*  http://www.acm.org/jgt/papers/GuigueDevillers03/                         
+*                                                                           
+*/
+
+/*
+Original MIT License (from https://github.com/erich666/jgt-code)
+Copyright (c) <year> <copyright holders>
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of this software
+and associated documentation files (the "Software"), to deal in the Software without
+restriction, including without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the
+Software is furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all copies or
+substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING
+BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
+DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+*/
+
+#pragma once
+#ifndef IGL_TRI_TRI_INTERSECT_H
+#define IGL_TRI_TRI_INTERSECT_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl {
+
+
+// Three-dimensional Triangle-Triangle overlap test
+//   if triangles are co-planar
+//
+// Input:
+//   p1,q1,r1  - vertices of the 1st triangle (3D)
+//   p2,q2,r2  - vertices of the 2nd triangle (3D)
+//
+// Output:
+// 
+//   Return true if two triangles overlap
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2> 
+IGL_INLINE bool tri_tri_overlap_test_3d(
+  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);
+
+
+// Three-dimensional Triangle-Triangle Intersection Test
+// additionaly computes the segment of intersection of the two triangles if it exists. 
+// coplanar returns whether the triangles are coplanar, 
+// source and target are the endpoints of the line segment of intersection 
+//
+// Input:
+//   p1,q1,r1  - vertices of the 1st triangle (3D)
+//   p2,q2,r2  - vertices of the 2nd triangle (3D)
+//
+// Output:
+//   coplanar - flag if two triangles are coplanar
+//   source - 1st point of intersection (if exists)
+//   target - 2nd point in intersection (if exists)
+//
+//   Return true if two triangles intersect
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2,
+typename DerivedS,typename DerivedT>
+IGL_INLINE bool tri_tri_intersection_test_3d(
+    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,
+    bool & coplanar, 
+    Eigen::MatrixBase<DerivedS> & source, 
+    Eigen::MatrixBase<DerivedT> & target );
+
+
+
+// Two dimensional Triangle-Triangle Overlap Test
+// Input:
+//   p1,q1,r1  - vertices of the 1st triangle (2D)
+//   p2,q2,r2  - vertices of the 2nd triangle (2D)
+//
+// Output:
+//   Return true if two triangles overlap
+template <typename DerivedP1,typename DerivedQ1,typename DerivedR1,
+typename DerivedP2,typename DerivedQ2,typename DerivedR2>
+IGL_INLINE bool tri_tri_overlap_test_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);
+
+
+};
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "tri_tri_intersect.cpp"
+#endif
+
+#endif // IGL_TRI_TRI_INTERSECT_H

tests/include/igl/fast_find_intersections.cpp

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

tutorial/723_FastFindSelfIntersections/main.cpp

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

tutorial/724_FastFindIntersections/main.cpp

@@ -0,0 +1,98 @@
+#include <igl/readOFF.h>
+#include <igl/combine.h>
+
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/opengl/glfw/imgui/ImGuiMenu.h>
+
+#include <igl/fast_find_intersections.h>
+
+Eigen::MatrixXd V1,V2;
+Eigen::MatrixXi F1,F2;
+
+igl::AABB<Eigen::MatrixXd,3> tree;
+
+double min_z,max_z;
+double slice_z;
+
+
+void update_visualization(igl::opengl::glfw::Viewer & viewer)
+{
+  Eigen::MatrixXi I;
+  Eigen::MatrixXd edges;
+
+  //shifted intersection object
+  Eigen::MatrixXd V2_(V2.rows(),V2.cols());
+  V2_<< V2.col(0), V2.col(1), V2.col(2).array()+slice_z;
+
+  igl::fast_find_intersections(tree, V1,F1, V2_,F2, I,edges);
+  Eigen::MatrixXi edges_link=Eigen::MatrixXi::NullaryExpr(edges.rows()/2,2, [](int i,int j) { return i*2+j;});
+ 
+  // Plot the edges of the intersects
+  viewer.data().set_edges ( edges, edges_link, Eigen::RowVector3d(1,0,0));
+  
+  // show faces which are intersected
+  Eigen::VectorXd face_data=Eigen::VectorXd::Zero(F1.rows());
+
+  for(int i=0; i<I.rows(); ++i)
+    face_data(I(i,0)) = 1.0;
+  
+  viewer.data().set_data(face_data);
+}
+
+
+bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int mod)
+{
+  switch(key)
+  {
+    default:
+      return false;
+    case '.':
+      slice_z = std::min(slice_z+0.01,max_z);
+      break;
+    case ',':
+      slice_z = std::max(slice_z-0.01,min_z);
+      break;
+  }
+  update_visualization(viewer);
+  return true;
+}
+
+
+int main(int argc, char *argv[])
+{
+  // Load two meshes 
+  igl::readOFF(TUTORIAL_SHARED_PATH "/cow.off",     V1, F1);
+  igl::readOFF(TUTORIAL_SHARED_PATH "/planexy.off", V2, F2);
+
+  // initialize AABB tree with first mesh (it doesn't move)
+  tree.init(V1,F1);
+
+  std::cout<<"Usage:"<<std::endl;
+  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;
+
+  //center slicing object
+  V2.col(2).array() -= V2.col(2).array().mean();
+
+  // Plot the mesh
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V1, F1);
+  
+
+  update_visualization(viewer);
+
+  igl::opengl::glfw::imgui::ImGuiMenu menu;
+  //plugin.widgets.push_back(&menu);
+  menu.callback_draw_viewer_window = [](){};
+  viewer.callback_key_down = &key_down;
+
+  // show the cow closer
+  viewer.core().camera_zoom = 2.0;
+
+  // Launch the viewer
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -116,4 +116,6 @@ if(LIBIGL_TUTORIALS_CHAPTER7)
     igl_add_tutorial(720_BlueNoise igl::glfw)
     igl_add_tutorial(721_VectorFieldSmoothing igl::glfw)
     igl_add_tutorial(722_VectorParallelTransport igl::glfw)
+    igl_add_tutorial(723_FastFindSelfIntersections igl::imgui igl::glfw)
+    igl_add_tutorial(724_FastFindIntersections igl::imgui igl::glfw)
 endif()