refactor to return meshed intersections

Former-commit-id: ecf4227e1df7db5ce705dfa514e2b008f9023fda

include/igl/cgal/intersect_other.cpp

@@ -8,123 +8,261 @@
 #include "intersect_other.h"
 #include "CGAL_includes.hpp"
 #include "mesh_to_cgal_triangle_list.h"
+#include "remesh_intersections.h"
 
 #ifndef IGL_FIRST_HIT_EXCEPTION
 #define IGL_FIRST_HIT_EXCEPTION 10
 #endif
 
-IGL_INLINE void igl::cgal::intersect_other(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  const Eigen::MatrixXd & U,
-  const Eigen::MatrixXi & G,
-  const bool first_only,
-  Eigen::MatrixXi & IF)
+// Un-exposed helper functions
+namespace igl
 {
-  using namespace std;
-  using namespace Eigen;
-
-
-  typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
-  // 3D Primitives
-  typedef CGAL::Point_3<Kernel>    Point_3;
-  typedef CGAL::Segment_3<Kernel>  Segment_3; 
-  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
-  typedef CGAL::Plane_3<Kernel>    Plane_3;
-  typedef CGAL::Tetrahedron_3<Kernel> Tetrahedron_3; 
-  // 2D Primitives
-  typedef CGAL::Point_2<Kernel>    Point_2;
-  typedef CGAL::Segment_2<Kernel>  Segment_2; 
-  typedef CGAL::Triangle_2<Kernel> Triangle_2; 
-  // 2D Constrained Delaunay Triangulation types
-  typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
-  typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
-  typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
-  typedef CGAL::Exact_intersections_tag Itag;
-  typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
-    CDT_2;
-  typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
-  // Axis-align boxes for all-pairs self-intersection detection
-  typedef std::vector<Triangle_3> Triangles;
-  typedef typename Triangles::iterator TrianglesIterator;
-  typedef typename Triangles::const_iterator TrianglesConstIterator;
-  typedef 
-    CGAL::Box_intersection_d::Box_with_handle_d<double,3,TrianglesIterator> 
-    Box;
-
-  Triangles TF,TG;
-  // Compute and process self intersections
-  mesh_to_cgal_triangle_list(V,F,TF);
-  mesh_to_cgal_triangle_list(U,G,TG);
-  // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
-  // Create the corresponding vector of bounding boxes
-  std::vector<Box> F_boxes,G_boxes;
-  const auto box_up = [](Triangles & T, std::vector<Box> & boxes)
+  namespace cgal
   {
-    boxes.reserve(T.size());
-    for ( 
-      TrianglesIterator tit = T.begin(); 
-      tit != T.end(); 
-      ++tit)
+    template <typename DerivedF>
+    static IGL_INLINE void push_result(
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      const int f,
+      const CGAL::Object & result,
+      std::map<
+        typename DerivedF::Index,
+        std::pair<typename DerivedF::Index,
+          std::vector<CGAL::Object> > > & offending,
+      std::map<
+        std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+        std::vector<typename DerivedF::Index> > & edge2faces)
     {
-      boxes.push_back(Box(tit->bbox(), tit));
-    }
-  };
-  box_up(TF,F_boxes);
-  box_up(TG,G_boxes);
-  std::list<int> lIF;
-  const auto cb = [&](const Box &a, const Box &b)
-  {
-    using namespace std;
-    // index in F and T
-    int fa = a.handle()-TF.begin();
-    int fb = b.handle()-TG.begin();
-    const Triangle_3 & A = *a.handle();
-    const Triangle_3 & B = *b.handle();
-    if(CGAL::do_intersect(A,B))
-    {
-      // There was an intersection
-      lIF.push_back(fa);
-      lIF.push_back(fb);
-      if(first_only)
+      typedef typename DerivedF::Index Index;
+      typedef std::pair<Index,Index> EMK;
+      if(offending.count(f) == 0)
       {
-        throw IGL_FIRST_HIT_EXCEPTION;
+        // first time marking, initialize with new id and empty list
+        Index id = offending.size();
+        offending[f] = {id,{}};
+        for(Index e = 0; e<3;e++)
+        {
+          // append face to edge's list
+          Index i = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+1)%3) : F(f,(e+2)%3);
+          Index j = F(f,(e+1)%3) < F(f,(e+2)%3) ? F(f,(e+2)%3) : F(f,(e+1)%3);
+          edge2faces[EMK(i,j)].push_back(f);
+        }
       }
+      offending[f].second.push_back(result);
     }
-  };
-  try{
-    CGAL::box_intersection_d(
-      F_boxes.begin(), F_boxes.end(),
-      G_boxes.begin(), G_boxes.end(),
-      cb);
-  }catch(int e)
-  {
-    // Rethrow if not FIRST_HIT_EXCEPTION
-    if(e != IGL_FIRST_HIT_EXCEPTION)
+    template <
+      typename Kernel,
+      typename DerivedVA,
+      typename DerivedFA,
+      typename DerivedVB,
+      typename DerivedFB,
+      typename DerivedIF,
+      typename DerivedVVA,
+      typename DerivedFFA,
+      typename DerivedJA,
+      typename DerivedIMA,
+      typename DerivedVVB,
+      typename DerivedFFB,
+      typename DerivedJB,
+      typename DerivedIMB>
+    static IGL_INLINE bool intersect_other_helper(
+      const Eigen::PlainObjectBase<DerivedVA> & VA,
+      const Eigen::PlainObjectBase<DerivedFA> & FA,
+      const Eigen::PlainObjectBase<DerivedVB> & VB,
+      const Eigen::PlainObjectBase<DerivedFB> & FB,
+      const RemeshSelfIntersectionsParam & params,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedVVA> & VVA,
+      Eigen::PlainObjectBase<DerivedFFA> & FFA,
+      Eigen::PlainObjectBase<DerivedJA>  & JA,
+      Eigen::PlainObjectBase<DerivedIMA> & IMA,
+      Eigen::PlainObjectBase<DerivedVVB> & VVB,
+      Eigen::PlainObjectBase<DerivedFFB> & FFB,
+      Eigen::PlainObjectBase<DerivedJB>  & JB,
+      Eigen::PlainObjectBase<DerivedIMB> & IMB)
     {
-      throw e;
+
+      using namespace std;
+      using namespace Eigen;
+
+      typedef typename DerivedFA::Index Index;
+      // 3D Primitives
+      typedef CGAL::Point_3<Kernel>    Point_3;
+      typedef CGAL::Segment_3<Kernel>  Segment_3; 
+      typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+      typedef CGAL::Plane_3<Kernel>    Plane_3;
+      typedef CGAL::Tetrahedron_3<Kernel> Tetrahedron_3; 
+      // 2D Primitives
+      typedef CGAL::Point_2<Kernel>    Point_2;
+      typedef CGAL::Segment_2<Kernel>  Segment_2; 
+      typedef CGAL::Triangle_2<Kernel> Triangle_2; 
+      // 2D Constrained Delaunay Triangulation types
+      typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
+      typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTAB_2;
+      typedef CGAL::Triangulation_data_structure_2<TVB_2,CTAB_2> TDS_2;
+      typedef CGAL::Exact_intersections_tag Itag;
+      typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
+        CDT_2;
+      typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
+      // Axis-align boxes for all-pairs self-intersection detection
+      typedef std::vector<Triangle_3> Triangles;
+      typedef typename Triangles::iterator TrianglesIterator;
+      typedef typename Triangles::const_iterator TrianglesConstIterator;
+      typedef 
+        CGAL::Box_intersection_d::Box_with_handle_d<double,3,TrianglesIterator> 
+        Box;
+      typedef 
+        std::map<Index,std::pair<Index,std::vector<CGAL::Object> > > 
+        OffendingMap;
+      typedef std::map<std::pair<Index,Index>,std::vector<Index> >  EdgeMap;
+      typedef std::pair<Index,Index> EMK;
+
+      Triangles TA,TB;
+      // Compute and process self intersections
+      mesh_to_cgal_triangle_list(VA,FA,TA);
+      mesh_to_cgal_triangle_list(VB,FB,TB);
+      // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
+      // Create the corresponding vector of bounding boxes
+      std::vector<Box> A_boxes,B_boxes;
+      const auto box_up = [](Triangles & T, std::vector<Box> & boxes)
+      {
+        boxes.reserve(T.size());
+        for ( 
+          TrianglesIterator tit = T.begin(); 
+          tit != T.end(); 
+          ++tit)
+        {
+          boxes.push_back(Box(tit->bbox(), tit));
+        }
+      };
+      box_up(TA,A_boxes);
+      box_up(TB,B_boxes);
+      OffendingMap offendingA,offendingB;
+      EdgeMap edge2facesA,edge2facesB;
+
+      std::list<int> lIF;
+      const auto cb = [&](const Box &a, const Box &b)
+      {
+        using namespace std;
+        // index in F and T
+        int fa = a.handle()-TA.begin();
+        int fb = b.handle()-TB.begin();
+        const Triangle_3 & A = *a.handle();
+        const Triangle_3 & B = *b.handle();
+        if(CGAL::do_intersect(A,B))
+        {
+          // There was an intersection
+          lIF.push_back(fa);
+          lIF.push_back(fb);
+          if(params.first_only)
+          {
+            throw IGL_FIRST_HIT_EXCEPTION;
+          }
+          if(!params.detect_only)
+          {
+            CGAL::Object result = CGAL::intersection(A,B);
+
+            push_result(FA,fa,result,offendingA,edge2facesA);
+            push_result(FB,fb,result,offendingB,edge2facesB);
+          }
+        }
+      };
+      try{
+        CGAL::box_intersection_d(
+          A_boxes.begin(), A_boxes.end(),
+          B_boxes.begin(), B_boxes.end(),
+          cb);
+      }catch(int e)
+      {
+        // Rethrow if not FIRST_HIT_EXCEPTION
+        if(e != IGL_FIRST_HIT_EXCEPTION)
+        {
+          throw e;
+        }
+        // Otherwise just fall through
+      }
+
+      // Convert lIF to Eigen matrix
+      assert(lIF.size()%2 == 0);
+      IF.resize(lIF.size()/2,2);
+      {
+        int i=0;
+        for(
+          list<int>::const_iterator ifit = lIF.begin();
+          ifit!=lIF.end();
+          )
+        {
+          IF(i,0) = (*ifit);
+          ifit++; 
+          IF(i,1) = (*ifit);
+          ifit++;
+          i++;
+        }
+      }
+      if(!params.detect_only)
+      {
+        remesh_intersections(VA,FA,TA,offendingA,edge2facesA,VVA,FFA,JA,IMA);
+        remesh_intersections(VB,FB,TB,offendingB,edge2facesB,VVB,FFB,JB,IMB);
+      }
+
+      return IF.rows() > 0;
     }
-    // Otherwise just fall through
   }
+}
 
-  // Convert lIF to Eigen matrix
-  assert(lIF.size()%2 == 0);
-  IF.resize(lIF.size()/2,2);
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedIF,
+  typename DerivedVVA,
+  typename DerivedFFA,
+  typename DerivedJA,
+  typename DerivedIMA,
+  typename DerivedVVB,
+  typename DerivedFFB,
+  typename DerivedJB,
+  typename DerivedIMB>
+IGL_INLINE bool igl::cgal::intersect_other(
+  const Eigen::PlainObjectBase<DerivedVA> & VA,
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<DerivedVB> & VB,
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  const RemeshSelfIntersectionsParam & params,
+  Eigen::PlainObjectBase<DerivedIF> & IF,
+  Eigen::PlainObjectBase<DerivedVVA> & VVA,
+  Eigen::PlainObjectBase<DerivedFFA> & FFA,
+  Eigen::PlainObjectBase<DerivedJA>  & JA,
+  Eigen::PlainObjectBase<DerivedIMA> & IMA,
+  Eigen::PlainObjectBase<DerivedVVB> & VVB,
+  Eigen::PlainObjectBase<DerivedFFB> & FFB,
+  Eigen::PlainObjectBase<DerivedJB>  & JB,
+  Eigen::PlainObjectBase<DerivedIMB> & IMB)
+{
+  if(params.detect_only)
   {
-    int i=0;
-    for(
-      list<int>::const_iterator ifit = lIF.begin();
-      ifit!=lIF.end();
-      )
-    {
-      IF(i,0) = (*ifit);
-      ifit++; 
-      IF(i,1) = (*ifit);
-      ifit++;
-      i++;
-    }
+    return intersect_other_helper<CGAL::Epick>
+      (VA,FA,VB,FB,params,IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+  }else
+  {
+    return intersect_other_helper<CGAL::Epeck>
+      (VA,FA,VB,FB,params,IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
   }
+}
 
+IGL_INLINE bool igl::cgal::intersect_other(
+  const Eigen::MatrixXd & VA,
+  const Eigen::MatrixXi & FA,
+  const Eigen::MatrixXd & VB,
+  const Eigen::MatrixXi & FB,
+  const bool first_only,
+  Eigen::MatrixXi & IF)
+{
+  Eigen::MatrixXd VVA,VVB;
+  Eigen::MatrixXi FFA,FFB;
+  Eigen::VectorXi JA,JB,IMA,IMB;
+  return intersect_other(
+    VA,FA,VB,FB,{true,first_only},IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/cgal/intersect_other.h

@@ -1,13 +1,14 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2015 Alec Jacobson <alecjacobson@gmail.com>
 // 
 // This Source Code Form is subject to the terms of the Mozilla Public License 
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_CGAL_INTERSECT_OTHER_H
 #define IGL_CGAL_INTERSECT_OTHER_H
-#include <igl/igl_inline.h>
+#include "../igl_inline.h"
+#include "RemeshSelfIntersectionsParam.h"
 
 #include <Eigen/Dense>
 
@@ -22,27 +23,76 @@ namespace igl
 {
   namespace cgal
   {
-    // INTERSECT Given a triangle mesh (V,F) and another mesh (U,G) find all pairs
-    // of intersecting faces. Note that self-intersections are ignored.
+    // INTERSECT_OTHER Given a triangle mesh (VA,FA) and another mesh (VB,FB)
+    // find all pairs of intersecting faces. Note that self-intersections are
+    // ignored.
     // 
-    // [VV,FF,IF] = selfintersect(V,F,'ParameterName',ParameterValue, ...)
+    // Inputs:
+    //   VA  #V by 3 list of vertex positions
+    //   FA  #F by 3 list of triangle indices into VA
+    //   VB  #V by 3 list of vertex positions
+    //   FB  #F by 3 list of triangle indices into VB
+    //   params   whether to detect only and then whether to only find first
+    //     intersection
+    // Outputs:
+    //   IF  #intersecting face pairs by 2 list of intersecting face pairs,
+    //     indexing FA and FB
+    //   VVA  #VVA by 3 list of vertex positions
+    //   FFA  #FFA by 3 list of triangle indices into VVA
+    //   JA  #FFA list of indices into FA denoting birth triangle
+    //   IMA  #VVA list of indices into VVA of unique vertices.
+    //   VVB  #VVB by 3 list of vertex positions
+    //   FFB  #FFB by 3 list of triangle indices into VVB
+    //   JB  #FFB list of indices into FB denoting birth triangle
+    //   IMB  #VVB list of indices into VVB of unique vertices.
+    template <
+      typename DerivedVA,
+      typename DerivedFA,
+      typename DerivedVB,
+      typename DerivedFB,
+      typename DerivedIF,
+      typename DerivedVVA,
+      typename DerivedFFA,
+      typename DerivedJA,
+      typename DerivedIMA,
+      typename DerivedVVB,
+      typename DerivedFFB,
+      typename DerivedJB,
+      typename DerivedIMB>
+    IGL_INLINE bool intersect_other(
+      const Eigen::PlainObjectBase<DerivedVA> & VA,
+      const Eigen::PlainObjectBase<DerivedFA> & FA,
+      const Eigen::PlainObjectBase<DerivedVB> & VB,
+      const Eigen::PlainObjectBase<DerivedFB> & FB,
+      const RemeshSelfIntersectionsParam & params,
+      Eigen::PlainObjectBase<DerivedIF> & IF,
+      Eigen::PlainObjectBase<DerivedVVA> & VVA,
+      Eigen::PlainObjectBase<DerivedFFA> & FFA,
+      Eigen::PlainObjectBase<DerivedJA>  & JA,
+      Eigen::PlainObjectBase<DerivedIMA> & IMA,
+      Eigen::PlainObjectBase<DerivedVVB> & VVB,
+      Eigen::PlainObjectBase<DerivedFFB> & FFB,
+      Eigen::PlainObjectBase<DerivedJB>  & JB,
+      Eigen::PlainObjectBase<DerivedIMB> & IMB);
+    // Legacy wrapper for detect only using common types.
     //
     // Inputs:
-    //   V  #V by 3 list of vertex positions
-    //   F  #F by 3 list of triangle indices into V
-    //   U  #U by 3 list of vertex positions
-    //   G  #G by 3 list of triangle indices into U
+    //   VA  #V by 3 list of vertex positions
+    //   FA  #F by 3 list of triangle indices into VA
+    //   VB  #V by 3 list of vertex positions
+    //   FB  #F by 3 list of triangle indices into VB
     //   first_only  whether to only detect the first intersection.
     // Outputs:
     //   IF  #intersecting face pairs by 2 list of intersecting face pairs,
-    //     indexing F and G
+    //     indexing FA and FB
+    // Returns true if any intersections were found
     //
-    // See also: selfintersect
-    IGL_INLINE void intersect_other(
-      const Eigen::MatrixXd & V,
-      const Eigen::MatrixXi & F,
-      const Eigen::MatrixXd & U,
-      const Eigen::MatrixXi & G,
+    // See also: remesh_self_intersections
+    IGL_INLINE bool intersect_other(
+      const Eigen::MatrixXd & VA,
+      const Eigen::MatrixXi & FA,
+      const Eigen::MatrixXd & VB,
+      const Eigen::MatrixXi & FB,
       const bool first_only,
       Eigen::MatrixXi & IF);
   }