Merge branch 'tri_tri_intersect' of github.com:vfonov/libigl into vfonov-tri_tri_intersect

include/igl/Guigue2003_tri_tri_intersect.cpp

@@ -0,0 +1,633 @@
+// 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
+*  
+*  Original notice: 
+*
+*  Triangle-Triangle Overlap Test Routines        
+*  July, 2002                                                          
+*  Updated December 2003                                                
+*                                                                       
+*  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 "Guigue2003_tri_tri_intersect.h"
+
+// helper functions
+namespace igl {
+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.0f) 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.0f) return false;\
+  else return true; }
+
+
+
+/* Permutation in a canonical form of T2's vertices */
+
+#define _IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2) { \
+  if (dp2 > 0.0f) { \
+     if (dq2 > 0.0f) _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2) \
+     else if (dr2 > 0.0f) _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.0f) { \
+    if (dq2 < 0.0f) _IGL_CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2)\
+    else if (dr2 < 0.0f) _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.0f) { \
+      if (dr2 >= 0.0f)  _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.0f) { \
+      if (dr2 > 0.0f) _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.0f) _IGL_CHECK_MIN_MAX(p1,q1,r1,r2,p2,q2)\
+      else if (dr2 < 0.0f) _IGL_CHECK_MIN_MAX(p1,r1,q1,r2,p2,q2)\
+      else return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);\
+     }}}
+  
+
+
+/*
+*
+*  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) */
+
+
+  _IGL_SUB(v1,p2,r2)
+  _IGL_SUB(v2,q2,r2)
+  _IGL_CROSS(N2,v1,v2)
+
+  _IGL_SUB(v1,p1,r2)
+  dp1 = _IGL_DOT(v1,N2);
+  _IGL_SUB(v1,q1,r2)
+  dq1 = _IGL_DOT(v1,N2);
+  _IGL_SUB(v1,r1,r2)
+  dr1 = _IGL_DOT(v1,N2);
+  
+  if (((dp1 * dq1) > 0.0f) && ((dp1 * dr1) > 0.0f))  return false; 
+
+  /* Compute distance signs  of p2, q2 and r2 to the plane of
+     triangle(p1,q1,r1) */
+
+  _IGL_SUB(v1,q1,p1)
+  _IGL_SUB(v2,r1,p1)
+  _IGL_CROSS(N1,v1,v2)
+
+  _IGL_SUB(v1,p2,r1)
+  dp2 = _IGL_DOT(v1,N1);
+  _IGL_SUB(v1,q2,r1)
+  dq2 = _IGL_DOT(v1,N1);
+  _IGL_SUB(v1,r2,r1)
+  dr2 = _IGL_DOT(v1,N1);
+  
+  if (((dp2 * dq2) > 0.0f) && ((dp2 * dr2) > 0.0f)) return false;
+
+  /* Permutation in a canonical form of T1's vertices */
+
+
+  if (dp1 > 0.0f) {
+    if (dq1 > 0.0f) _IGL_TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2)
+    else if (dr1 > 0.0f) _IGL_TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2)  
+    else _IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2)
+  } else if (dp1 < 0.0f) {
+    if (dq1 < 0.0f) _IGL_TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2)
+    else if (dr1 < 0.0f) _IGL_TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2)
+    else _IGL_TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2)
+  } else {
+    if (dq1 < 0.0f) {
+      if (dr1 >= 0.0f) _IGL_TRI_TRI_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2)
+      else _IGL_TRI_TRI_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2)
+    }
+    else if (dq1 > 0.0f) {
+      if (dr1 > 0.0f) _IGL_TRI_TRI_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2)
+      else _IGL_TRI_TRI_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2)
+    }
+    else  {
+      if (dr1 > 0.0f) _IGL_TRI_TRI_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2)
+      else if (dr1 < 0.0f) _IGL_TRI_TRI_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2)
+      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::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);
+    
+};
+
+
+
+/*
+*                                                                
+*  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
+
+#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.0f) {\
+    _IGL_SUB(v1,r1,p1) \
+    _IGL_CROSS(N,v1,v2) \
+    if (_IGL_DOT(v,N) <= 0.0f) { \
+      _IGL_SUB(v2,q2,p1) \
+      _IGL_CROSS(N,v1,v2) \
+      if (_IGL_DOT(v,N) > 0.0f) { \
+  _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.0f) { \
+      return false; \
+    } else { \
+      _IGL_SUB(v1,r1,p1) \
+      _IGL_CROSS(N,v1,v2) \
+      if (_IGL_DOT(v,N) >= 0.0f) { \
+  _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; \
+      }}}} 
+
+                
+
+#define _IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2) { \
+  if (dp2 > 0.0f) { \
+     if (dq2 > 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,r2,p2,q2) \
+     else if (dr2 > 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,q2,r2,p2)\
+     else _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2) }\
+  else if (dp2 < 0.0f) { \
+    if (dq2 < 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,r2,p2,q2)\
+    else if (dr2 < 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,q2,r2,p2)\
+    else _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,p2,q2,r2)\
+  } else { \
+    if (dq2 < 0.0f) { \
+      if (dr2 >= 0.0f)  _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,q2,r2,p2)\
+      else _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,p2,q2,r2)\
+    } \
+    else if (dq2 > 0.0f) { \
+      if (dr2 > 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,p2,q2,r2)\
+      else  _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,q2,r2,p2)\
+    } \
+    else  { \
+      if (dr2 > 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,q1,r1,r2,p2,q2)\
+      else if (dr2 < 0.0f) _IGL_CONSTRUCT_INTERSECTION(p1,r1,q1,r2,p2,q2)\
+      else { \
+        coplanar = true; \
+  return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);\
+     } \
+  }} }
+  
+
+/*
+   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)
+
+
+  _IGL_SUB(v1,p2,r2)
+  _IGL_SUB(v2,q2,r2)
+  _IGL_CROSS(N2,v1,v2)
+
+  _IGL_SUB(v1,p1,r2)
+  dp1 = _IGL_DOT(v1,N2);
+  _IGL_SUB(v1,q1,r2)
+  dq1 = _IGL_DOT(v1,N2);
+  _IGL_SUB(v1,r1,r2)
+  dr1 = _IGL_DOT(v1,N2);
+  
+  coplanar = false;
+
+  if (((dp1 * dq1) > 0.0f) && ((dp1 * dr1) > 0.0f))  return false; 
+
+  // Compute distance signs  of p2, q2 and r2 
+  // to the plane of triangle(p1,q1,r1)
+
+  
+  _IGL_SUB(v1,q1,p1)
+  _IGL_SUB(v2,r1,p1)
+  _IGL_CROSS(N1,v1,v2)
+
+  _IGL_SUB(v1,p2,r1)
+  dp2 = _IGL_DOT(v1,N1);
+  _IGL_SUB(v1,q2,r1)
+  dq2 = _IGL_DOT(v1,N1);
+  _IGL_SUB(v1,r2,r1)
+  dr2 = _IGL_DOT(v1,N1);
+  
+  if (((dp2 * dq2) > 0.0f) && ((dp2 * dr2) > 0.0f)) return false;
+
+  // Permutation in a canonical form of T1's vertices
+
+
+  if (dp1 > 0.0f) {
+    if (dq1 > 0.0f) _IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2)
+    else if (dr1 > 0.0f) _IGL_TRI_TRI_INTER_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2)
+  
+    else _IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2)
+  } else if (dp1 < 0.0f) {
+    if (dq1 < 0.0f) _IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2)
+    else if (dr1 < 0.0f) _IGL_TRI_TRI_INTER_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2)
+    else _IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2)
+  } else {
+    if (dq1 < 0.0f) {
+      if (dr1 >= 0.0f) _IGL_TRI_TRI_INTER_3D(q1,r1,p1,p2,r2,q2,dp2,dr2,dq2)
+      else _IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,q2,r2,dp2,dq2,dr2)
+    }
+    else if (dq1 > 0.0f) {
+      if (dr1 > 0.0f) _IGL_TRI_TRI_INTER_3D(p1,q1,r1,p2,r2,q2,dp2,dr2,dq2)
+      else _IGL_TRI_TRI_INTER_3D(q1,r1,p1,p2,q2,r2,dp2,dq2,dr2)
+    }
+    else  {
+      if (dr1 > 0.0f) _IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,q2,r2,dp2,dq2,dr2)
+      else if (dr1 < 0.0f) _IGL_TRI_TRI_INTER_3D(r1,p1,q1,p2,r2,q2,dp2,dr2,dq2)
+      else {
+  // triangles are co-planar
+
+  coplanar = true;
+  return coplanar_tri_tri3d(p1,q1,r1,p2,q2,r2,N1);
+      }
+    }
+  }
+};
+
+
+
+
+
+/*
+*
+*  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]))
+
+
+#define INTERSECTION_TEST_VERTEX(P1, Q1, R1, P2, Q2, R2) {\
+  if (_IGL_ORIENT_2D(R2,P2,Q1) >= 0.0f)\
+    if (_IGL_ORIENT_2D(R2,Q2,Q1) <= 0.0f)\
+      if (_IGL_ORIENT_2D(P1,P2,Q1) > 0.0f) {\
+  if (_IGL_ORIENT_2D(P1,Q2,Q1) <= 0.0f) return true; \
+  else return false;} else {\
+  if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0f)\
+    if (_IGL_ORIENT_2D(Q1,R1,P2) >= 0.0f) return true; \
+    else return false;\
+  else return false;}\
+    else \
+      if (_IGL_ORIENT_2D(P1,Q2,Q1) <= 0.0f)\
+  if (_IGL_ORIENT_2D(R2,Q2,R1) <= 0.0f)\
+    if (_IGL_ORIENT_2D(Q1,R1,Q2) >= 0.0f) return true; \
+    else return false;\
+  else return false;\
+      else return false;\
+  else\
+    if (_IGL_ORIENT_2D(R2,P2,R1) >= 0.0f) \
+      if (_IGL_ORIENT_2D(Q1,R1,R2) >= 0.0f)\
+  if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0f) return true;\
+  else return false;\
+      else \
+  if (_IGL_ORIENT_2D(Q1,R1,Q2) >= 0.0f) {\
+    if (_IGL_ORIENT_2D(R2,R1,Q2) >= 0.0f) 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.0f) {\
+    if (_IGL_ORIENT_2D(P1,P2,Q1) >= 0.0f) { \
+        if (_IGL_ORIENT_2D(P1,Q1,R2) >= 0.0f) return true; \
+        else return false;} else { \
+      if (_IGL_ORIENT_2D(Q1,R1,P2) >= 0.0f){ \
+  if (_IGL_ORIENT_2D(R1,P1,P2) >= 0.0f) return true; else return false;} \
+      else return false; } \
+  } else {\
+    if (_IGL_ORIENT_2D(R2,P2,R1) >= 0.0f) {\
+      if (_IGL_ORIENT_2D(P1,P2,R1) >= 0.0f) {\
+  if (_IGL_ORIENT_2D(P1,R1,R2) >= 0.0f) return true;  \
+  else {\
+    if (_IGL_ORIENT_2D(Q1,R1,R2) >= 0.0f) 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 igl::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.0f ) {
+    if ( _IGL_ORIENT_2D(q2,r2,p1) >= 0.0f ) {
+      if ( _IGL_ORIENT_2D(r2,p2,p1) >= 0.0f ) return true;
+      else _IGL_INTERSECTION_TEST_EDGE(p1,q1,r1,p2,q2,r2)
+    } else {  
+      if ( _IGL_ORIENT_2D(r2,p2,p1) >= 0.0f ) 
+  _IGL_INTERSECTION_TEST_EDGE(p1,q1,r1,r2,p2,q2)
+      else INTERSECTION_TEST_VERTEX(p1,q1,r1,p2,q2,r2)}}
+  else {
+    if ( _IGL_ORIENT_2D(q2,r2,p1) >= 0.0f ) {
+      if ( _IGL_ORIENT_2D(r2,p2,p1) >= 0.0f ) 
+  _IGL_INTERSECTION_TEST_EDGE(p1,q1,r1,q2,r2,p2)
+      else  INTERSECTION_TEST_VERTEX(p1,q1,r1,q2,r2,p2)}
+    else INTERSECTION_TEST_VERTEX(p1,q1,r1,r2,p2,q2)}
+};
+
+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_ORIENT_2D(p1,q1,r1) < 0.0f )
+    if ( _IGL_ORIENT_2D(p2,q2,r2) < 0.0f )
+      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 ( _IGL_ORIENT_2D(p2,q2,r2) < 0.0f )
+      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);
+
+};
+
+//cleanup
+#undef _IGL_ORIENT_2D
+#undef _IGL_INTERSECTION_TEST_VERTEX
+#undef _IGL_TRI_TRI_INTER_3D
+#undef _IGL_SUB
+#undef _IGL_DOT
+#undef _IGL_CROSS
+#undef _IGL_SCALAR
+#undef _IGL_CHECK_MIN_MAX
+#undef _IGL_INTERSECTION_TEST_EDGE
+
+
+#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> >&);
+#endif

include/igl/Guigue2003_tri_tri_intersect.h

@@ -0,0 +1,86 @@
+// 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/
+#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 "Guigue2003_tri_tri_intersect.cpp"
+#endif
+
+#endif // IGL_TRI_TRI_INTERSECT_H

include/igl/fast_find_mesh_intersect.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_mesh_intersect.h"
+#include "AABB.h"
+#include "Guigue2003_tri_tri_intersect.h"
+
+
+
+template <
+  typename DerivedV1,
+  typename DerivedF1,
+  typename DerivedV2,
+  typename DerivedF2,
+  typename DerivedI,
+  typename DerivedE>
+IGL_INLINE void igl::fast_find_mesh_intersect(
+  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_mesh_intersect(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_mesh_intersect(
+  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_mesh_intersect<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_mesh_intersect<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_mesh_intersect.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_mesh_intersect(
+    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_mesh_intersect(
+    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_mesh_intersect.cpp"
+#endif
+
+#endif

include/igl/fast_find_mesh_selfintersect.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_mesh_selfintersect.h"
+#include "AABB.h"
+#include "Guigue2003_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_mesh_selfintersect(
+  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_mesh_selfintersect(
+  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_mesh_selfintersect<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_mesh_selfintersect<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_mesh_selfintersect.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_MESH_SELFINTERSECT_H
+#define FAST_FIND_MESH_SELFINTERSECT_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_mesh_selfintersect(
+      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_mesh_selfintersect(
+      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_mesh_selfintersect.cpp"
+#endif
+
+#endif

tests/include/igl/mesh_intersections.cpp

@@ -0,0 +1,86 @@
+#include <test_common.h>
+#include <igl/fast_find_mesh_selfintersect.h>
+#include <igl/fast_find_mesh_intersect.h>
+#include <igl/combine.h>
+
+TEST_CASE("fast_find_mesh_selfintersect: 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_mesh_selfintersect(V,F,I) );
+
+  REQUIRE ( I.sum()==0);
+}
+
+TEST_CASE("fast_find_mesh_selfintersect: 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_mesh_selfintersect(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_mesh_intersect:", "[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_mesh_intersect(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_ShowSelfIntersections/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_mesh_selfintersect.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_mesh_selfintersect(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_ShowIntersections/main.cpp

@@ -0,0 +1,99 @@
+#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_mesh_intersect.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_mesh_intersect(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_ShowSelfIntersections igl::imgui igl::glfw)
+    igl_add_tutorial(724_ShowIntersections igl::imgui igl::glfw)
 endif()