libigl/include/igl/triangle_triangle_intersect.cpp

#include "triangle_triangle_intersect.h"
#include "triangle_triangle_intersect_shared_edge.h"
#include "triangle_triangle_intersect_shared_vertex.h"
#include "tri_tri_intersect.h"
#include <Eigen/Geometry>
#include <stdio.h>
#include <igl/unique_edge_map.h>
#include <igl/barycentric_coordinates.h>
#include <unordered_map>

//#define IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
// CGAL::Epeck
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#warning "🐌🐌🐌🐌🐌🐌🐌🐌 Slow debug mode for igl::triangle_triangle_intersect"
#endif

template <
  typename DerivedV,
  typename DerivedF,
  typename DerivedE,
  typename DerivedEMAP,
  typename DerivedEF,
  typename Derivedp>
IGL_INLINE bool igl::triangle_triangle_intersect(
  const Eigen::MatrixBase<DerivedV> & V,
  const Eigen::MatrixBase<DerivedF> & F,
  const Eigen::MatrixBase<DerivedE> & E,
  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
  const Eigen::MatrixBase<DerivedEF> & EF,
  const int f,
  const int c,
  const Eigen::MatrixBase<Derivedp> & p,
  const int g)
{
  static_assert(
    std::is_same<typename DerivedV::Scalar,typename Derivedp::Scalar>::value, 
    "V and p should have same Scalar type");
  assert(V.cols() == 3);
  assert(p.cols() == 3);
#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
  using Kernel = CGAL::Epeck;
  typedef CGAL::Point_3<Kernel>    Point_3;
  typedef CGAL::Segment_3<Kernel>  Segment_3;
  typedef CGAL::Triangle_3<Kernel> Triangle_3;
    bool cgal_found_intersection = false;
    Point_3 Vg[3];
    Point_3 Vf[3];
    for(int i = 0;i<3;i++)
    {
      Vg[i] = Point_3(V(F(g,i),0),V(F(g,i),1),V(F(g,i),2));
      if(i == c)
      {
        Vf[i] = Point_3(p(0),p(1),p(2));
      }else
      {
        Vf[i] = Point_3(V(F(f,i),0),V(F(f,i),1),V(F(f,i),2));
      }
    }
    Triangle_3 Tg(Vg[0],Vg[1],Vg[2]);
    Triangle_3 Tf(Vf[0],Vf[1],Vf[2]);
#endif

  // I'm leaving this debug printing stuff in for a bit until I trust this
  // better. 
  constexpr bool stinker = false;
  //const bool stinker = (f==1492 && g==1554);
  if(stinker) { printf("👀\n"); }
  bool found_intersection = false;
  // So edge opposite F(f,c) is the outer edge.
  const bool share_edge_opposite_c = [&]()
  {
    const int o = EMAP(f + c*F.rows());
    return (EF(o,0) == f && EF(o,1) == g) || (EF(o,1) == f && EF(o,0) == g);
  }();
  const int o = EMAP(f + c*F.rows());
  // Do they share the edge opposite c?
  if(share_edge_opposite_c)
  {
    if(stinker) { printf("⚠️ shares an edge\n"); }
    found_intersection = triangle_triangle_intersect_shared_edge(V,F,f,c,p,g,1e-8);
  }else
  {
    if(stinker) { printf("does not share an edge\n"); }
    // Do they share a vertex?
    int sf,sg;
    bool found_shared_vertex = false;
    for(sf = 0;sf<3;sf++)
    {
      if(sf == c){ continue;}
      for(sg = 0;sg<3;sg++)
      {
        if(F(f,sf) == F(g,sg))
        {
          found_shared_vertex = true;
          break;
        }
      }
      if(found_shared_vertex) { break;} 
    }
    if(found_shared_vertex)
    {
      if(stinker) { printf("⚠️ shares a vertex\n"); }
      found_intersection = 
        triangle_triangle_intersect_shared_vertex(V,F,f,sf,c,p,g,sg,1e-14);
    }else
    {
      bool coplanar;
      Eigen::RowVector3d i1,i2;
      found_intersection = 
        igl::tri_tri_overlap_test_3d(
                V.row(F(g,0)).template cast<double>(), 
                V.row(F(g,1)).template cast<double>(), 
                V.row(F(g,2)).template cast<double>(),
                            p.template cast<double>(),
          V.row(F(f,(c+1)%3)).template cast<double>(),
          V.row(F(f,(c+2)%3)).template cast<double>()) && 
        igl::tri_tri_intersection_test_3d(
                V.row(F(g,0)).template cast<double>(), 
                V.row(F(g,1)).template cast<double>(), 
                V.row(F(g,2)).template cast<double>(),
                            p.template cast<double>(),
          V.row(F(f,(c+1)%3)).template cast<double>(),
          V.row(F(f,(c+2)%3)).template cast<double>(),
          coplanar,
          i1,i2);
      if(stinker) { printf("tri_tri_intersection_test_3d says %s\n",found_intersection?"☠️":"✅"); }
#ifdef IGL_TRIANGLE_TRIANGLE_INTERSECT_DEBUG
      if(CGAL::do_intersect(Tg,Tf))
      {
        cgal_found_intersection = true;
        printf("  ✅ sure it's anything\n");
      }
      assert(found_intersection == cgal_found_intersection);
#endif
    }
  }
  if(stinker) { printf("%s\n",found_intersection?"☠️":"✅"); }
  return found_intersection;
}


template <
  typename DerivedV, 
  typename DerivedF, 
  typename DerivedIF,
  typename DerivedEV,
  typename DerivedEE>
void igl::triangle_triangle_intersect(
    const Eigen::MatrixBase<DerivedV> & V1,
    const Eigen::MatrixBase<DerivedF> & F1,
    const Eigen::MatrixBase<DerivedV> & V2,
    const Eigen::MatrixBase<DerivedF> & F2,
    const Eigen::MatrixBase<DerivedIF> & IF,
    Eigen::PlainObjectBase<DerivedEV> & EV,
    Eigen::PlainObjectBase<DerivedEE> & EE)
{
  using Scalar = typename DerivedEV::Scalar;
  using RowVector3S = Eigen::Matrix<Scalar,1,3>;
  // We were promised that IF is a list of non-coplanar non-degenerately
  // intersecting triangle pairs. This implies that the set of intersection is a
  // line-segment whose endpoints are defined by edge-triangle intersections.
  // 
  // Each edge-triangle intersection will be stored in a map (e,f)
  //
  std::unordered_map<std::int64_t,int> uf_to_ev;
  Eigen::VectorXi EMAP1;
  Eigen::MatrixXi uE1;
  {
    Eigen::VectorXi uEE,uEC;
    Eigen::MatrixXi E;
    igl::unique_edge_map(F1,E,uE1,EMAP1,uEE,uEC);
  }
  Eigen::VectorXi EMAP2_cpy;
  Eigen::MatrixXi uE2_cpy;
  const bool self = &F1 == &F2;
  if(!self)
  {
    Eigen::VectorXi uEE,uEC;
    Eigen::MatrixXi E;
    igl::unique_edge_map(F2,E,uE2_cpy,EMAP2_cpy,uEE,uEC);
  }
  const Eigen::VectorXi & EMAP2 = self?EMAP1:EMAP2_cpy;
  const Eigen::MatrixXi & uE2 = self?uE1:uE2_cpy;

  int num_ev = 0;
  EE.resize(IF.rows(),2);
  for(int i = 0; i<IF.rows(); i++)
  {
    // Just try all 6 edges
    Eigen::Matrix<double,6,6> B;
    Eigen::Matrix<double,6,3> X;
    Eigen::Matrix<int,6,2> uf;
    for(int p = 0;p<2;p++)
    {
      const auto consider_edges = [&B,&X,&uf]
        (const int p,
         const int f1, 
         const int f2,
         const Eigen::MatrixBase<DerivedV> & V1,
         const Eigen::MatrixBase<DerivedF> & F1,
         const Eigen::VectorXi & EMAP1,
         const Eigen::MatrixXi & uE1,
         const Eigen::MatrixBase<DerivedV> & V2,
         const Eigen::MatrixBase<DerivedF> & F2)
      {
        for(int e1 = 0;e1<3;e1++)
        {
          // intersect edge (ij) opposite vertex F(f1,e1) with triangle ABC of F(f2,:)
          const int u1 = EMAP1(f1 +(EMAP1.size()/3)*e1);
          uf.row(p*3+e1) << u1,f2;
          const int i = uE1(u1,0);
          const int j = uE1(u1,1);
          // Just copy.
          const RowVector3S Vi = V1.row(i);
          const RowVector3S Vj = V1.row(j);
          const RowVector3S VA = V2.row(F2(f2,0));
          const RowVector3S VB = V2.row(F2(f2,1));
          const RowVector3S VC = V2.row(F2(f2,2));
          // Find intersection of line (Vi,Vj) with plane of triangle (A,B,C)
          const RowVector3S n = (VB-VA).template head<3>().cross((VC-VA).template head<3>());
          const Scalar d = n.dot(VA);
          const Scalar t = (d - n.dot(Vi))/(n.dot(Vj-Vi));
          const RowVector3S x = Vi + t*(Vj-Vi);
          
          // Get barycenteric coordinates (possibly negative of X)
          RowVector3S b1;
          igl::barycentric_coordinates(x,VA,VB,VC,b1);
          B.row(p*3+e1).head<3>() = b1;
          RowVector3S b2;
          igl::barycentric_coordinates(x,
              V1.row(F1(f1,0)).template head<3>().eval(),
              V1.row(F1(f1,1)).template head<3>().eval(),
              V1.row(F1(f1,2)).template head<3>().eval(),
              b2);
          B.row(p*3+e1).tail<3>() = b2;

          X.row(p*3+e1) = x;
        }
      };


      const int f1 = IF(i,p);
      const int f2 = IF(i,(p+1)%2);
      consider_edges(p,f1,f2,
        p==0?   V1:V2,
        p==0?   F1:F2,
        p==0?EMAP1:EMAP2,
        p==0?  uE1:uE2,
        p==0?   V2:V1,
        p==0?   F2:F1);
    }

    // Find the two rows in B with the largest-smallest element
    int j1,j2;
    {
      double b_min1 = -std::numeric_limits<double>::infinity();
      double b_min2 = -std::numeric_limits<double>::infinity();
      for(int j = 0;j<6;j++)
      {
        // It's not clear that using barycentric coordinates is better than
        // point_simplex distance (though that requires inequalities).
        const double bminj = B.row(j).minCoeff();
        if(bminj > b_min1)
        {
          b_min2 = b_min1;
          j2 = j1;
          b_min1 = bminj;
          j1 = j;
        }else if(bminj > b_min2)
        {
          b_min2 = bminj;
          j2 = j;
        }
      }
    }

    const auto append_or_find = [&](
      int p, int u, int f, const RowVector3S & x,
      std::unordered_map<std::int64_t,int> & uf_to_ev)->int
    {
      const std::int64_t key = (std::int64_t)u + ((std::int64_t)f << 32) + ((std::int64_t)p << 63);
      if(uf_to_ev.find(key) == uf_to_ev.end())
      {
        if(num_ev == EV.rows())
        {
          EV.conservativeResize(2*EV.rows()+1,3);
        }
        EV.row(num_ev) = x;
        uf_to_ev[key] = num_ev;
        num_ev++;
      }
      return uf_to_ev[key];
    };

    EE.row(i) <<
      append_or_find(j1>=3,uf(j1,0),uf(j1,1),X.row(j1),uf_to_ev),
      append_or_find(j2>=3,uf(j2,0),uf(j2,1),X.row(j2),uf_to_ev);
  }
  EV.conservativeResize(num_ev,3);
}

template <
  typename DerivedV, 
  typename DerivedF, 
  typename DerivedIF,
  typename DerivedEV,
  typename DerivedEE>
void igl::triangle_triangle_intersect(
    const Eigen::MatrixBase<DerivedV> & V,
    const Eigen::MatrixBase<DerivedF> & F,
    const Eigen::MatrixBase<DerivedIF> & IF,
    Eigen::PlainObjectBase<DerivedEV> & EV,
    Eigen::PlainObjectBase<DerivedEE> & EE)
{
  // overload will take care of detecting reference equality
  return triangle_triangle_intersect(V,F,V,F,IF,EV,EE);
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::triangle_triangle_intersect<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
// generated by autoexplicit.sh
template void igl::triangle_triangle_intersect<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
template bool igl::triangle_triangle_intersect<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, -1, false> > const&, int);
template bool igl::triangle_triangle_intersect<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, int, Eigen::MatrixBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&, int);
#endif