trim with solid mesh

Former-commit-id: e7273737423f05c4f75f89e9628abe8cbb8b53be

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -109,7 +109,7 @@ namespace igl
           // Make a short name for the edge map
           typedef std::map<EMK,EMV> EdgeMap;
           // Maps edges of offending faces to all incident offending faces
-          EdgeMap edge2faces;
+          //EdgeMap edge2faces;
           std::vector<std::pair<const Box, const Box> > candidate_box_pairs;
 
         public:
@@ -330,7 +330,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   T(),
   lIF(),
   offending(),
-  edge2faces(),
+  //edge2faces(),
   params(params)
 {
   using namespace std;
@@ -427,7 +427,7 @@ inline igl::copyleft::cgal::SelfIntersectMesh<
   }
 
   remesh_intersections(
-    V,F,T,offending,edge2faces,params.stitch_all,VV,FF,J,IM);
+    V,F,T,offending,params.stitch_all,VV,FF,J,IM);
 
 #ifdef IGL_SELFINTERSECTMESH_DEBUG
   log_time("remesh_intersection");
@@ -477,7 +477,7 @@ inline void igl::copyleft::cgal::SelfIntersectMesh<
       // 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);
+      //edge2faces[EMK(i,j)].push_back(f);
     }
   }
 }

include/igl/copyleft/cgal/intersect_other.cpp

@@ -9,6 +9,8 @@
 #include "CGAL_includes.hpp"
 #include "mesh_to_cgal_triangle_list.h"
 #include "remesh_intersections.h"
+#include "../../slice_mask.h"
+#include "../../remove_unreferenced.h"
 
 #ifndef IGL_FIRST_HIT_EXCEPTION
 #define IGL_FIRST_HIT_EXCEPTION 10
@@ -30,10 +32,10 @@ namespace igl
         std::map<
           typename DerivedF::Index,
           std::vector<std::pair<typename DerivedF::Index, CGAL::Object> > > &
-          offending,
-        std::map<
-          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-          std::vector<typename DerivedF::Index> > & edge2faces)
+          offending)
+        //std::map<
+        //  std::pair<typename DerivedF::Index,typename DerivedF::Index>,
+        //  std::vector<typename DerivedF::Index> > & edge2faces)
       {
         typedef typename DerivedF::Index Index;
         typedef std::pair<Index,Index> EMK;
@@ -46,7 +48,7 @@ namespace igl
             // 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);
+            //edge2faces[EMK(i,j)].push_back(f);
           }
         }
         offending[f].push_back({f_other,result});
@@ -58,14 +60,10 @@ namespace igl
         typename DerivedVB,
         typename DerivedFB,
         typename DerivedIF,
-        typename DerivedVVA,
-        typename DerivedFFA,
-        typename DerivedJA,
-        typename DerivedIMA,
-        typename DerivedVVB,
-        typename DerivedFFB,
-        typename DerivedJB,
-        typename DerivedIMB>
+        typename DerivedVVAB,
+        typename DerivedFFAB,
+        typename DerivedJAB,
+        typename DerivedIMAB>
       static IGL_INLINE bool intersect_other_helper(
         const Eigen::PlainObjectBase<DerivedVA> & VA,
         const Eigen::PlainObjectBase<DerivedFA> & FA,
@@ -73,14 +71,10 @@ namespace igl
         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)
+        Eigen::PlainObjectBase<DerivedVVAB> & VVAB,
+        Eigen::PlainObjectBase<DerivedFFAB> & FFAB,
+        Eigen::PlainObjectBase<DerivedJAB>  & JAB,
+        Eigen::PlainObjectBase<DerivedIMAB> & IMAB)
       {
 
         using namespace std;
@@ -139,7 +133,7 @@ namespace igl
         box_up(TA,A_boxes);
         box_up(TB,B_boxes);
         OffendingMap offendingA,offendingB;
-        EdgeMap edge2facesA,edge2facesB;
+        //EdgeMap edge2facesA,edge2facesB;
 
         std::list<int> lIF;
         const auto cb = [&](const Box &a, const Box &b) -> void
@@ -163,8 +157,8 @@ namespace igl
             {
               CGAL::Object result = CGAL::intersection(A,B);
 
-              push_result(FA,fa,fb,result,offendingA,edge2facesA);
-              push_result(FB,fb,fa,result,offendingB,edge2facesB);
+              push_result(FA,fa,fb,result,offendingA);
+              push_result(FB,fb,fa,result,offendingB);
             }
           }
         };
@@ -202,8 +196,38 @@ namespace igl
         }
         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);
+          // Obsolete, now remesh_intersections expects a single mesh
+          // remesh_intersections(VA,FA,TA,offendingA,VVA,FFA,JA,IMA);
+          // remesh_intersections(VB,FB,TB,offendingB,VVB,FFB,JB,IMB);
+          // Combine mesh and offending maps
+          DerivedVA VAB(VA.rows()+VB.rows(),3);
+          VAB<<VA,VB;
+          DerivedFA FAB(FA.rows()+FB.rows(),3);
+          FAB<<FA,(FB.array()+VA.rows());
+          Triangles TAB;
+          TAB.reserve(TA.size()+TB.size());
+          TAB.insert(TAB.end(),TA.begin(),TA.end());
+          TAB.insert(TAB.end(),TB.begin(),TB.end());
+          OffendingMap offending;
+          //offending.reserve(offendingA.size() + offendingB.size());
+          for (const auto itr : offendingA)
+          {
+            // Remap offenders in FB to FAB
+            auto offenders = itr.second;
+            for(auto & offender : offenders)
+            {
+              offender.first += FA.rows();
+            }
+            offending[itr.first] = offenders;
+          }
+          for (const auto itr : offendingB)
+          {
+            // Store offenders for FB according to place in FAB
+            offending[FA.rows() + itr.first] = itr.second;
+          }
+
+          remesh_intersections(
+            VAB,FAB,TAB,offending,params.stitch_all,VVAB,FFAB,JAB,IMAB);
         }
 
         return IF.rows() > 0;
@@ -218,38 +242,30 @@ template <
   typename DerivedVB,
   typename DerivedFB,
   typename DerivedIF,
-  typename DerivedVVA,
-  typename DerivedFFA,
-  typename DerivedJA,
-  typename DerivedIMA,
-  typename DerivedVVB,
-  typename DerivedFFB,
-  typename DerivedJB,
-  typename DerivedIMB>
+  typename DerivedVVAB,
+  typename DerivedFFAB,
+  typename DerivedJAB,
+  typename DerivedIMAB>
 IGL_INLINE bool igl::copyleft::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)
+    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<DerivedVVAB> & VVAB,
+    Eigen::PlainObjectBase<DerivedFFAB> & FFAB,
+    Eigen::PlainObjectBase<DerivedJAB>  & JAB,
+    Eigen::PlainObjectBase<DerivedIMAB> & IMAB)
 {
   if(params.detect_only)
   {
     return intersect_other_helper<CGAL::Epick>
-      (VA,FA,VB,FB,params,IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+      (VA,FA,VB,FB,params,IF,VVAB,FFAB,JAB,IMAB);
   }else
   {
     return intersect_other_helper<CGAL::Epeck>
-      (VA,FA,VB,FB,params,IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+      (VA,FA,VB,FB,params,IF,VVAB,FFAB,JAB,IMAB);
   }
 }
 
@@ -261,11 +277,11 @@ IGL_INLINE bool igl::copyleft::cgal::intersect_other(
   const bool first_only,
   Eigen::MatrixXi & IF)
 {
-  Eigen::MatrixXd VVA,VVB;
-  Eigen::MatrixXi FFA,FFB;
-  Eigen::VectorXi JA,JB,IMA,IMB;
+  Eigen::MatrixXd VVAB;
+  Eigen::MatrixXi FFAB;
+  Eigen::VectorXi JAB,IMAB;
   return intersect_other(
-    VA,FA,VB,FB,{true,first_only},IF,VVA,FFA,JA,IMA,VVB,FFB,JB,IMB);
+    VA,FA,VB,FB,{true,first_only},IF,VVAB,FFAB,JAB,IMAB);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/copyleft/cgal/intersect_other.h

@@ -39,28 +39,21 @@ namespace igl
       // 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.
+      //   VVAB  #VVAB by 3 list of vertex positions
+      //   FFAB  #FFAB by 3 list of triangle indices into VVA
+      //   JAB  #FFAB list of indices into [FA;FB] denoting birth triangle
+      //   IMAB  #VVAB list of indices stitching duplicates (resulting from
+      //     mesh intersections) together
       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>
+        typename DerivedVVAB,
+        typename DerivedFFAB,
+        typename DerivedJAB,
+        typename DerivedIMAB>
       IGL_INLINE bool intersect_other(
         const Eigen::PlainObjectBase<DerivedVA> & VA,
         const Eigen::PlainObjectBase<DerivedFA> & FA,
@@ -68,14 +61,10 @@ namespace igl
         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);
+        Eigen::PlainObjectBase<DerivedVVAB> & VVAB,
+        Eigen::PlainObjectBase<DerivedFFAB> & FFAB,
+        Eigen::PlainObjectBase<DerivedJAB>  & JAB,
+        Eigen::PlainObjectBase<DerivedIMAB> & IMAB);
       // Legacy wrapper for detect only using common types.
       //
       // Inputs:

include/igl/copyleft/cgal/point_mesh_squared_distance.cpp

@@ -8,14 +8,21 @@
 #include "point_mesh_squared_distance.h"
 #include "mesh_to_cgal_triangle_list.h"
 
-template <typename Kernel>
+template <
+  typename Kernel,
+  typename DerivedP,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedC>
 IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
-  const Eigen::MatrixXd & P,
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
-  Eigen::VectorXd & sqrD,
-  Eigen::VectorXi & I,
-  Eigen::MatrixXd & C)
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+        Eigen::PlainObjectBase<DerivedI> & I,
+        Eigen::PlainObjectBase<DerivedC> & C)
 {
   using namespace std;
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
@@ -29,10 +36,10 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
   return point_mesh_squared_distance(P,tree,T,sqrD,I,C);
 }
 
-template <typename Kernel>
+template <typename Kernel, typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance_precompute(
-  const Eigen::MatrixXd & V,
-  const Eigen::MatrixXi & F,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
   CGAL::AABB_tree<
     CGAL::AABB_traits<Kernel, 
       CGAL::AABB_triangle_primitive<Kernel, 
@@ -77,9 +84,14 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance_precompute(
   tree.closest_point_and_primitive(Point_3(0,0,0));
 }
 
-template <typename Kernel>
+template <
+  typename Kernel,
+  typename DerivedP,
+  typename DerivedsqrD,
+  typename DerivedI,
+  typename DerivedC>
 IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
-  const Eigen::MatrixXd & P,
+  const Eigen::PlainObjectBase<DerivedP> & P,
   const CGAL::AABB_tree<
     CGAL::AABB_traits<Kernel, 
       CGAL::AABB_triangle_primitive<Kernel, 
@@ -88,9 +100,9 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
     >
   > & tree,
   const std::vector<CGAL::Triangle_3<Kernel> > & T,
-  Eigen::VectorXd & sqrD,
-  Eigen::VectorXi & I,
-  Eigen::MatrixXd & C)
+  Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C)
 {
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
   typedef typename std::vector<Triangle_3>::iterator Iterator;
@@ -110,18 +122,14 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance(
     // Find closest point and primitive id
     Point_and_primitive_id pp = tree.closest_point_and_primitive(query);
     Point_3 closest_point = pp.first;
-    C(p,0) = CGAL::to_double(closest_point[0]);
-    C(p,1) = CGAL::to_double(closest_point[1]);
-    C(p,2) = CGAL::to_double(closest_point[2]);
-    sqrD(p) = CGAL::to_double((closest_point-query).squared_length());
+    assign_scalar(closest_point[0],C(p,0));
+    assign_scalar(closest_point[1],C(p,1));
+    assign_scalar(closest_point[2],C(p,2));
+    assign_scalar((closest_point-query).squared_length(),sqrD(p));
     I(p) = pp.second - T.begin();
   }
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
-template void igl::copyleft::cgal::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
-template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
-template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template void igl::copyleft::cgal::point_mesh_squared_distance<CGAL::Simple_cartesian<double> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, 1, 0, -1, 1>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
 #endif

include/igl/copyleft/cgal/point_mesh_squared_distance.h

@@ -33,24 +33,35 @@ namespace igl
       //
       // Known bugs: This only computes distances to triangles. So unreferenced
       // vertices and degenerate triangles (segments) are ignored.
-      template <typename Kernel>
+      template <
+        typename Kernel,
+        typename DerivedP,
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedsqrD,
+        typename DerivedI,
+        typename DerivedC>
       IGL_INLINE void point_mesh_squared_distance(
-        const Eigen::MatrixXd & P,
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & F,
-        Eigen::VectorXd & sqrD,
-        Eigen::VectorXi & I,
-        Eigen::MatrixXd & C);
+        const Eigen::PlainObjectBase<DerivedP> & P,
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
+              Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+              Eigen::PlainObjectBase<DerivedI> & I,
+              Eigen::PlainObjectBase<DerivedC> & C);
       // Probably can do this in a way that we don't pass around `tree` and `T`
       //
       // Outputs:
       //   tree  CGAL's AABB tree
       //   T  list of CGAL triangles in order of F (for determining which was found
       //     in computation)
-      template <typename Kernel>
+      template <
+        typename Kernel,
+        typename DerivedV,
+        typename DerivedF
+        >
       IGL_INLINE void point_mesh_squared_distance_precompute(
-        const Eigen::MatrixXd & V,
-        const Eigen::MatrixXi & F,
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedF> & F,
         CGAL::AABB_tree<
           CGAL::AABB_traits<Kernel, 
             CGAL::AABB_triangle_primitive<Kernel, 
@@ -63,9 +74,14 @@ namespace igl
       //  see above
       // Outputs:
       //  see above
-      template <typename Kernel>
+      template <
+        typename Kernel,
+        typename DerivedP,
+        typename DerivedsqrD,
+        typename DerivedI,
+        typename DerivedC>
       IGL_INLINE void point_mesh_squared_distance(
-        const Eigen::MatrixXd & P,
+        const Eigen::PlainObjectBase<DerivedP> & P,
         const CGAL::AABB_tree<
           CGAL::AABB_traits<Kernel, 
             CGAL::AABB_triangle_primitive<Kernel, 
@@ -74,9 +90,9 @@ namespace igl
           >
         > & tree,
         const std::vector<CGAL::Triangle_3<Kernel> > & T,
-        Eigen::VectorXd & sqrD,
-        Eigen::VectorXi & I,
-        Eigen::MatrixXd & C);
+        Eigen::PlainObjectBase<DerivedsqrD> & sqrD,
+        Eigen::PlainObjectBase<DerivedI> & I,
+        Eigen::PlainObjectBase<DerivedC> & C);
     }
   }
 }

include/igl/copyleft/cgal/point_solid_signed_squared_distance.cpp

@@ -0,0 +1,50 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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/.
+#include "point_solid_signed_squared_distance.h"
+#include "points_inside_component.h"
+#include "point_mesh_squared_distance.h"
+#include "../../list_to_matrix.h"
+#include "../../slice_mask.h"
+#include <vector>
+#include <Eigen/Core>
+
+template <
+  typename DerivedQ,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedD>
+IGL_INLINE void igl::copyleft::cgal::point_solid_signed_squared_distance(
+  const Eigen::PlainObjectBase<DerivedQ> & Q,
+  const Eigen::PlainObjectBase<DerivedVB> & VB,
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  Eigen::PlainObjectBase<DerivedD> & D)
+{
+  // compute unsigned distances
+  Eigen::VectorXi I;
+  DerivedVB C;
+  point_mesh_squared_distance<CGAL::Epeck>(Q,VB,FB,D,I,C);
+  // Collect queries that have non-zero distance
+  Eigen::Array<bool,Eigen::Dynamic,1> NZ = D.array()!=0;
+  // Compute sign for non-zero distance queries
+  DerivedQ QNZ;
+  slice_mask(Q,NZ,1,QNZ);
+  Eigen::Array<bool,Eigen::Dynamic,1> DNZ;
+  igl::copyleft::cgal::points_inside_component(VB,FB,QNZ,DNZ);
+  // Apply sign to distances
+  DerivedD S = DerivedD::Zero(Q.rows(),1);
+  {
+    int k = 0;
+    for(int q = 0;q<Q.rows();q++)
+    {
+      if(NZ(q))
+      {
+        D(q) *= DNZ(k++) ? -1. : 1.;
+      }
+    }
+  }
+}

include/igl/copyleft/cgal/point_solid_signed_squared_distance.h

@@ -0,0 +1,48 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_COPYLEFT_CGAL_POINT_SOLID_SIGNED_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_POINT_SOLID_SIGNED_SQUARED_DISTANCE_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // POINT_SOLID_SIGNED_SQUARED_DISTANCE Given a set of points (Q) and the
+      // boundary mesh (VB,FB) of a solid (as defined in [Zhou et al. 2016],
+      // determine the signed squared distance for each point q in Q so that d(q,B) is
+      // negative if inside and positive if outside.
+      //
+      // Inputs:
+      //   Q  #Q by 3 list of query point positions
+      //   VB  #VB by 3 list of mesh vertex positions of B
+      //   FB  #FB by 3 list of mesh triangle indices into VB
+      // Outputs:
+      //   D
+      template <
+        typename DerivedQ,
+        typename DerivedVB,
+        typename DerivedFB,
+        typename DerivedD>
+      IGL_INLINE void point_solid_signed_squared_distance(
+        const Eigen::PlainObjectBase<DerivedQ> & Q,
+        const Eigen::PlainObjectBase<DerivedVB> & VB,
+        const Eigen::PlainObjectBase<DerivedFB> & FB,
+        Eigen::PlainObjectBase<DerivedD> & D);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_solid_signed_squared_distance.cpp"
+#endif
+
+#endif

include/igl/copyleft/cgal/remesh_intersections.h

@@ -30,7 +30,6 @@ namespace igl
       //   offending #offending map taking face indices into F to pairs of order
       //     of first finding and list of intersection objects from all
       //     intersections
-      //   edge2faces  #edges <= #offending*3 to incident offending faces 
       //   stitch_all  if true, merge all vertices with thte same coordiante.
       // Outputs:
       //   VV  #VV by 3 list of vertex positions, if stitch_all = false then
@@ -58,9 +57,6 @@ namespace igl
           typename DerivedF::Index,
             std::vector<
             std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
-        const std::map<
-          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-          std::vector<typename DerivedF::Index> > & edge2faces,
         bool stitch_all,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
@@ -83,9 +79,6 @@ namespace igl
           typename DerivedF::Index,
             std::vector<
             std::pair<typename DerivedF::Index, CGAL::Object> > > & offending,
-        const std::map<
-          std::pair<typename DerivedF::Index,typename DerivedF::Index>,
-          std::vector<typename DerivedF::Index> > & edge2faces,
         Eigen::PlainObjectBase<DerivedVV> & VV,
         Eigen::PlainObjectBase<DerivedFF> & FF,
         Eigen::PlainObjectBase<DerivedJ> & J,

include/igl/copyleft/cgal/trim_with_solid.cpp

@@ -0,0 +1,106 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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/.
+#include "trim_with_solid.h"
+#include "assign_scalar.h"
+#include "intersect_other.h"
+#include "point_solid_signed_squared_distance.h"
+
+#include "../../extract_manifold_patches.h"
+#include "../../list_to_matrix.h"
+#include "../../remove_unreferenced.h"
+#include "../../slice_mask.h"
+
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#include <vector>
+
+template <
+  typename DerivedVA,
+  typename DerivedFA,
+  typename DerivedVB,
+  typename DerivedFB,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedD,
+  typename DerivedJ>
+IGL_INLINE void igl::copyleft::cgal::trim_with_solid(
+  const Eigen::PlainObjectBase<DerivedVA> & VA,
+  const Eigen::PlainObjectBase<DerivedFA> & FA,
+  const Eigen::PlainObjectBase<DerivedVB> & VB,
+  const Eigen::PlainObjectBase<DerivedFB> & FB,
+  Eigen::PlainObjectBase<DerivedV> & Vd,
+  Eigen::PlainObjectBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedD> & D,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  // resolve intersections using exact representation
+  typedef Eigen::Matrix<CGAL::Epeck::FT,Eigen::Dynamic,3> MatrixX3E;
+  typedef Eigen::Matrix<CGAL::Epeck::FT,Eigen::Dynamic,1> VectorXE;
+  typedef Eigen::Matrix<CGAL::Epeck::FT,1,3> RowVector3E;
+  MatrixX3E V;
+  Eigen::MatrixXi _1;
+  Eigen::VectorXi _2;
+  // Intersect A and B meshes and stitch together new faces
+  igl::copyleft::cgal::intersect_other(
+    VA,FA,VB,FB,{false,false,true},_1,V,F,J,_2);
+  // Partition result into manifold patches
+  Eigen::VectorXi P;
+  const size_t num_patches = igl::extract_manifold_patches(F,P);
+  // only keep faces from A
+  Eigen::Matrix<bool,Eigen::Dynamic,1> A = J.array()< FA.rows();
+  igl::slice_mask(Eigen::MatrixXi(F),A,1,F);
+  igl::slice_mask(Eigen::VectorXi(P),A,1,P);
+  igl::slice_mask(Eigen::VectorXi(J),A,1,J);
+  // Agregate representative query points for each patch
+  std::vector<bool> flag(num_patches);
+  std::vector<std::vector<CGAL::Epeck::FT> > vQ;
+  Eigen::VectorXi P2Q(num_patches);
+  for(int f = 0;f<P.rows();f++)
+  {
+    const auto p = P(f);
+    // if not yet processed this patch
+    if(!flag[p])
+    {
+      P2Q(p) = vQ.size();
+      std::vector<CGAL::Epeck::FT> q = {
+        (V(F(f,0),0)+ V(F(f,1),0)+ V(F(f,2),0))/3.,
+        (V(F(f,0),1)+ V(F(f,1),1)+ V(F(f,2),1))/3.,
+        (V(F(f,0),2)+ V(F(f,1),2)+ V(F(f,2),2))/3.};
+      vQ.emplace_back(q);
+      flag[p] = true;
+    }
+  }
+  MatrixX3E Q;
+  igl::list_to_matrix(vQ,Q);
+  VectorXE SP;
+  point_solid_signed_squared_distance(Q,VB,FB,SP);
+  Eigen::Matrix<bool,Eigen::Dynamic,1> DP = SP.array()>0;
+  // distribute flag to all faces
+  D.resize(F.rows());
+  for(int f = 0;f<F.rows();f++)
+  {
+    D(f) = DP(P2Q(P(f)));
+  }
+  Eigen::VectorXi _;
+  igl::remove_unreferenced(MatrixX3E(V),DerivedF(F),V,F,_);
+  const auto & assign = [](
+    const MatrixX3E & V, 
+    Eigen::PlainObjectBase<DerivedV> & Vd)
+  {
+    Vd.resize(V.rows(),3);
+    for(int v = 0;v<V.rows();v++)
+    {
+      for(int d = 0;d<3;d++) 
+      {
+        igl::copyleft::cgal::assign_scalar(V(v,d),Vd(v,d));
+      }
+    }
+  };
+  assign(V,Vd);
+}
+

include/igl/copyleft/cgal/trim_with_solid.h

@@ -0,0 +1,63 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 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_COPYLEFT_CGAL_TRIM_WITH_SOLID_H
+#define IGL_COPYLEFT_CGAL_TRIM_WITH_SOLID_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // TRIM_WITH_SOLID Given an arbitrary mesh (VA,FA) and the boundary mesh
+      // (VB,FB) of a solid (as defined in [Zhou et al. 2016]), Resolve intersections
+      // between A and B subdividing faces of A so that intersections with B exists
+      // only along edges and vertices (and coplanar faces). Then determine whether
+      // each of these faces is inside or outside of B. This can be used to extract
+      // the part of A inside or outside of B.
+      //
+      // Inputs:
+      //   VA  #VA by 3 list of mesh vertex positions of A
+      //   FA  #FA by 3 list of mesh triangle indices into VA
+      //   VB  #VB by 3 list of mesh vertex positions of B
+      //   FB  #FB by 3 list of mesh triangle indices into VB
+      // Outputs:
+      //   V  #V by 3 list of mesh vertex positions of output
+      //   F  #F by 3 list of mesh triangle indices into V
+      //   D  #F list of bools whether face is inside B
+      //   J  #F list of indices into FA revealing birth parent
+      //
+      template <
+        typename DerivedVA,
+        typename DerivedFA,
+        typename DerivedVB,
+        typename DerivedFB,
+        typename DerivedV,
+        typename DerivedF,
+        typename DerivedD,
+        typename DerivedJ>
+      IGL_INLINE void trim_with_solid(
+        const Eigen::PlainObjectBase<DerivedVA> & VA,
+        const Eigen::PlainObjectBase<DerivedFA> & FA,
+        const Eigen::PlainObjectBase<DerivedVB> & VB,
+        const Eigen::PlainObjectBase<DerivedFB> & FB,
+        Eigen::PlainObjectBase<DerivedV> & Vd,
+        Eigen::PlainObjectBase<DerivedF> & F,
+        Eigen::PlainObjectBase<DerivedD> & D,
+        Eigen::PlainObjectBase<DerivedJ> & J);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "trim_with_solid.cpp"
+#endif
+#endif

include/igl/extract_manifold_patches.cpp

@@ -1,18 +1,20 @@
 #include "extract_manifold_patches.h"
+#include "unique_edge_map.h"
 #include <cassert>
 #include <limits>
 #include <queue>
 
 template<
-typename DerivedF,
-typename DerivedEMAP,
-typename uE2EType,
-typename DerivedP>
+  typename DerivedF,
+  typename DerivedEMAP,
+  typename uE2EType,
+  typename DerivedP>
 IGL_INLINE size_t igl::extract_manifold_patches(
-        const Eigen::PlainObjectBase<DerivedF>& F,
-        const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
-        const std::vector<std::vector<uE2EType> >& uE2E,
-        Eigen::PlainObjectBase<DerivedP>& P) {
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
+  const std::vector<std::vector<uE2EType> >& uE2E,
+  Eigen::PlainObjectBase<DerivedP>& P) 
+{
     assert(F.cols() == 3);
     const size_t num_faces = F.rows();
 
@@ -67,6 +69,20 @@ IGL_INLINE size_t igl::extract_manifold_patches(
     return num_patches;
 }
 
+template<
+  typename DerivedF,
+  typename DerivedP>
+IGL_INLINE size_t igl::extract_manifold_patches(
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedP>& P) 
+{
+  Eigen::MatrixXi E, uE;
+  Eigen::VectorXi EMAP;
+  std::vector<std::vector<size_t> > uE2E;
+  igl::unique_edge_map(F, E, uE, EMAP, uE2E);
+  return igl::extract_manifold_patches(F, EMAP, uE2E, P);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 template unsigned long igl::extract_manifold_patches<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, unsigned long, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, std::vector<std::vector<unsigned long, std::allocator<unsigned long> >, std::allocator<std::vector<unsigned long, std::allocator<unsigned long> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/extract_manifold_patches.h

@@ -30,6 +30,12 @@ namespace igl {
       const Eigen::PlainObjectBase<DerivedEMAP>& EMAP,
       const std::vector<std::vector<uE2EType> >& uE2E,
       Eigen::PlainObjectBase<DerivedP>& P);
+    template <
+      typename DerivedF,
+      typename DerivedP>
+    IGL_INLINE size_t extract_manifold_patches(
+      const Eigen::PlainObjectBase<DerivedF>& F,
+      Eigen::PlainObjectBase<DerivedP>& P);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "extract_manifold_patches.cpp"