snap rounding, segment-segment intersection resolution

Former-commit-id: 1a094f6f7b2396d5a0c13e59055a48a2df6e2fb2

include/igl/copyleft/cgal/point_segment_squared_distance.cpp

@@ -0,0 +1,34 @@
+#include "point_segment_squared_distance.h"
+
+template < typename Kernel>
+IGL_INLINE void igl::copyleft::cgal::point_segment_squared_distance(
+  const CGAL::Point_3<Kernel> & P1,
+  const CGAL::Segment_3<Kernel> & S2,
+  CGAL::Point_3<Kernel> & P2,
+  typename Kernel::FT & d)
+{
+  if(S2.is_degenerate())
+  {
+    P2 = S2.source();
+    d = (P1-P2).squared_length();
+    return;
+  }
+  // http://stackoverflow.com/a/1501725/148668
+  const auto sqr_len = S2.squared_length();
+  assert(sqr_len != 0);
+  const auto & V = S1.source();
+  const auto & W = S2.target();
+  const auto t = (P1-V).dot(W-V)/sqr_len;
+  if(t<0)
+  {
+    P2 = V;
+  }else if(t>1)
+  {
+    P2 = W;
+  }else
+  {
+    P2 = V  + t*(W-V);
+  }
+  d = (P1-P2).squared_length();
+}
+

include/igl/copyleft/cgal/point_segment_squared_distance.h

@@ -0,0 +1,37 @@
+#ifndef IGL_COPYLEFT_CGAL_POINT_SEGMENT_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_POINT_SEGMENT_SQUARED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <CGAL/Segment_3.h>
+#include <CGAL/Point_3.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given a point P1 and segment S2 find the points on each of closest
+      // approach and the squared distance thereof.
+      // 
+      // Inputs:
+      //   P1  point
+      //   S2  segment
+      // Outputs:
+      //   P2  point on S2 closest to P1
+      //   d  distance betwee P1 and S2
+      template < typename Kernel>
+      IGL_INLINE void point_segment_squared_distance(
+          const CGAL::Point_3<Kernel> & P1,
+          const CGAL::Segment_3<Kernel> & S2,
+          CGAL::Point_3<Kernel> & P2,
+          typename Kernel::FT & d
+          );
+
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_segment_squared_distance.cpp"
+#endif
+
+#endif
+

include/igl/copyleft/cgal/point_triangle_squared_distance.cpp

@@ -0,0 +1,39 @@
+#include "point_triangle_squared_distance.h"
+template < typename Kernel>
+IGL_INLINE void point_triangle_squared_distance(
+  const CGAL::Point_3<Kernel> & P1,
+  const CGAL::Triangle_3<Kernel> & T2,
+  CGAL::Point_3<Kernel> & P2,
+  typename Kernel::FT & d)
+{
+  assert(!T2.is_degenerate());
+  if(T2.has_on(P1))
+  {
+    P2 = P1;
+    d = 0;
+    return;
+  }
+  const auto proj_1 = T2.supporting_plane().projection(P2);
+  if(T2.has_on(proj_1))
+  {
+    P2 = proj_1;
+    d = (proj_1-P1).squared_length();
+    return;
+  }
+  // closest point must be on the boundary
+  bool first = true;
+  // loop over edges
+  for(int i=0;i<3;i++)
+  {
+    Point_3<Kernel> P2i;
+    typename Kernel::FT di;
+    const Segment_3<Kernel> si( T2.vertex(i+1), T2.vertex(i+2));
+    point_segment_squared_distance(P1,si,P2i,di);
+    if(first || di < d)
+    {
+      first = false;
+      d = di;
+      P2 = P2i;
+    }
+  }
+}

include/igl/copyleft/cgal/point_triangle_squared_distance.h

@@ -0,0 +1,38 @@
+#ifndef IGL_COPYLEFT_CGAL_POINT_TRIANGLE_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_POINT_TRIANGLE_SQUARED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <CGAL/Triangle_3.h>
+#include <CGAL/Point_3.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given a point P1 and triangle T2 find the points on each of closest
+      // approach and the squared distance thereof.
+      // 
+      // Inputs:
+      //   P1  point
+      //   T2  triangle
+      // Outputs:
+      //   P2  point on T2 closest to P1
+      //   d  distance betwee P1 and T2
+      template < typename Kernel>
+      IGL_INLINE void point_triangle_squared_distance(
+        const CGAL::Point_3<Kernel> & P1,
+        const CGAL::Triangle_3<Kernel> & T2,
+        CGAL::Point_3<Kernel> & P2,
+        typename Kernel::FT & d
+        );
+
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_triangle_squared_distance.cpp"
+#endif
+
+#endif
+
+

include/igl/copyleft/cgal/resolve_intersections.cpp

@@ -0,0 +1,84 @@
+#include "resolve_intersections.h"
+#include "subdivide_segments.h"
+#include "row_to_point.h"
+#include <igl/unique.h>
+#include <igl/list_to_matrix.h>
+#include <igl/copyleft/cgal/assign_scalar.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Segment_2.h>
+#include <CGAL/Point_2.h>
+#include <algorithm>
+#include <vector>
+
+template <
+  typename DerivedV, 
+  typename DerivedE, 
+  typename DerivedVI, 
+  typename DerivedEI,
+  typename DerivedJ,
+  typename DerivedIM>
+IGL_INLINE void igl::copyleft::cgal::resolve_intersections(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedVI> & VI,
+  Eigen::PlainObjectBase<DerivedEI> & EI,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<DerivedIM> & IM)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  // Exact scalar type
+  typedef CGAL::Epeck K;
+  typedef K::FT EScalar;
+  typedef CGAL::Segment_2<K> Segment_2;
+  typedef CGAL::Point_2<K> Point_2;
+  typedef Matrix<EScalar,Dynamic,Dynamic>  MatrixXE;
+
+  // Convert vertex positions to exact kernel
+  MatrixXE VE(V.rows(),V.cols());
+  for(int i = 0;i<V.rows();i++)
+  {
+    for(int j = 0;j<V.cols();j++)
+    {
+      VE(i,j) = V(i,j);
+    }
+  }
+
+  const int m = E.rows();
+  // resolve all intersections: silly O(n²) implementation
+  std::vector<std::vector<Point_2> > steiner(m);
+  for(int i = 0;i<m;i++)
+  {
+    Segment_2 si(row_to_point<K>(VE,E(i,0)),row_to_point<K>(VE,E(i,1)));
+    steiner[i].push_back(si.vertex(0));
+    steiner[i].push_back(si.vertex(1));
+    for(int j = i+1;j<m;j++)
+    {
+      Segment_2 sj(row_to_point<K>(VE,E(j,0)),row_to_point<K>(VE,E(j,1)));
+      // do they intersect?
+      if(CGAL::do_intersect(si,sj))
+      {
+        CGAL::Object result = CGAL::intersection(si,sj);
+        if(const Point_2 * p = CGAL::object_cast<Point_2 >(&result))
+        {
+          steiner[i].push_back(*p);
+          steiner[j].push_back(*p);
+          // add intersection point
+        }else if(const Segment_2 * s = CGAL::object_cast<Segment_2 >(&result))
+        {
+          // add both endpoints
+          steiner[i].push_back(s->vertex(0));
+          steiner[j].push_back(s->vertex(0));
+          steiner[i].push_back(s->vertex(1));
+          steiner[j].push_back(s->vertex(1));
+        }else
+        {
+          assert(false && "Unknown intersection type");
+        }
+      }
+    }
+  }
+
+  subdivide_segments(V,E,steiner,VI,EI,J,IM);
+}

include/igl/copyleft/cgal/resolve_intersections.h

@@ -0,0 +1,44 @@
+#ifndef IGL_COPYLEFT_CGAL_RESOLVE_INTERSECTIONS_H
+#define IGL_COPYLEFT_CGAL_RESOLVE_INTERSECTIONS_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    // RESOLVE_INTERSECTIONS Given a list of possible intersecting segments with
+    // endpoints, split segments to overlap only at endpoints
+    //
+    // Inputs:
+    //   V  #V by 2 list of vertex positions
+    //   E  #E by 2 list of segment indices into V
+    // Outputs:
+    //   VI  #VI by 2 list of output vertex positions, copies of V are always
+    //     the first #V vertices
+    //   EI  #EI by 2 list of segment indices into V, #EI ≥ #E
+    //   J  #EI list of indices into E revealing "parent segments"
+    //   IM  #VI list of indices into VV of unique vertices.
+    namespace cgal
+    {
+      template <
+        typename DerivedV, 
+        typename DerivedE, 
+        typename DerivedVI, 
+        typename DerivedEI,
+        typename DerivedJ,
+        typename DerivedIM>
+      IGL_INLINE void resolve_intersections(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedE> & E,
+        Eigen::PlainObjectBase<DerivedVI> & VI,
+        Eigen::PlainObjectBase<DerivedEI> & EI,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<DerivedIM> & IM);
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "resolve_intersections.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/row_to_point.cpp

@@ -0,0 +1,11 @@
+#include "row_to_point.h"
+
+template <
+  typename Kernel,
+  typename DerivedV>
+IGL_INLINE CGAL::Point_2<Kernel> igl::copyleft::cgal::row_to_point(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const typename DerivedV::Index & i)
+{
+  return CGAL::Point_2<Kernel>(V(i,0),V(i,1));
+}

include/igl/copyleft/cgal/row_to_point.h

@@ -0,0 +1,31 @@
+#ifndef IGL_COPYLEFT_CGAL_ROW_TO_POINT_H
+#define IGL_COPYLEFT_CGAL_ROW_TO_POINT_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <CGAL/Point_2.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Extract a row from V and treat as a 2D cgal point (only first two
+      // columns of V are used).
+      // 
+      // Inputs:
+      //   V  #V by 2 list of vertex positions
+      //   i  row index
+      // Returns 2D cgal point
+      template <
+        typename Kernel,
+        typename DerivedV>
+      IGL_INLINE CGAL::Point_2<Kernel> row_to_point(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const typename DerivedV::Index & i);
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "row_to_point.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/segment_segment_squared_distance.cpp

@@ -0,0 +1,121 @@
+#include "segment_segment_squared_distance.h"
+
+// http://geomalgorithms.com/a07-_distance.html
+template < typename Kernel>
+IGL_INLINE bool igl::copyleft::cgal::segment_segment_squared_distance(
+    const CGAL::Segment_3<Kernel> & S1,
+    const CGAL::Segment_3<Kernel> & S2,
+    CGAL::Point_3<Kernel> & P1,
+    CGAL::Point_3<Kernel> & P2,
+    typename Kernel::FT & d)
+{
+  typedef CGAL::Point_3<Kernel> Point_3;
+  typedef CGAL::Vector_3<Kernel> Vector_3;
+  typedef Kernel::FT EScalar;
+  if(S1.is_degenerate())
+  {
+    // All points on S1 are the same
+    P1 = S1.source();
+    point_segment_squared_distance(P1,S2,P2,d);
+    return true;
+  }else if(S2.is_degenerate())
+  {
+    assert(!S1.is_degenerate());
+    // All points on S2 are the same
+    P2 = S2.source();
+    point_segment_squared_distance(P2,S1,P1,d);
+    return true;
+  }
+
+  assert(!S1.is_degenerate());
+  assert(!S2.is_degenerate());
+
+  Vector_3 u = S1.target() - S1.source();
+  Vector_3 v = S2.target() - S2.source();
+  Vector_3 w = S1.source() - S2.source();
+
+  const auto a = u.dot(u);         // always >= 0
+  const auto b = u.dot(v);
+  const auto c = v.dot(v);         // always >= 0
+  const auto d = u.dot(w);
+  const auto e = v.dot(w);
+  const auto D = a*c - b*b;        // always >= 0
+  assert(D>=0);
+  const auto sc, sN, sD = D;       // sc = sN / sD, default sD = D >= 0
+  const auto tc, tN, tD = D;       // tc = tN / tD, default tD = D >= 0
+
+  bool parallel = false;
+  // compute the line parameters of the two closest points
+  if (D==0) 
+  { 
+    // the lines are almost parallel
+    parallel = true;
+    sN = 0.0;         // force using source point on segment S1
+    sD = 1.0;         // to prevent possible division by 0.0 later
+    tN = e;
+    tD = c;
+  } else
+  {
+    // get the closest points on the infinite lines
+    sN = (b*e - c*d);
+    tN = (a*e - b*d);
+    if (sN < 0.0) 
+    { 
+      // sc < 0 => the s=0 edge is visible
+      sN = 0.0;
+      tN = e;
+      tD = c;
+    } else if (sN > sD) 
+    {  // sc > 1  => the s=1 edge is visible
+      sN = sD;
+      tN = e + b;
+      tD = c;
+    }
+  }
+
+  if (tN < 0.0) 
+  {
+    // tc < 0 => the t=0 edge is visible
+    tN = 0.0;
+    // recompute sc for this edge
+    if (-d < 0.0)
+    {
+      sN = 0.0;
+    }else if (-d > a)
+    {
+      sN = sD;
+    }else 
+    {
+      sN = -d;
+      sD = a;
+    }
+  }else if (tN > tD) 
+  {
+    // tc > 1  => the t=1 edge is visible
+    tN = tD;
+    // recompute sc for this edge
+    if ((-d + b) < 0.0)
+    {
+      sN = 0;
+    }else if ((-d + b) > a)
+    {
+      sN = sD;
+    }else
+    {
+      sN = (-d +  b);
+      sD = a;
+    }
+  }
+  // finally do the division to get sc and tc
+  sc = (abs(sN) == 0 ? 0.0 : sN / sD);
+  tc = (abs(tN) == 0 ? 0.0 : tN / tD);
+
+  // get the difference of the two closest points
+  P1 = S1.source() + sc*(S1.target()-S1.source());
+  P2 = S2.source() + sc*(S2.target()-S2.source());
+  Vector_3   dP = w + (sc * u) - (tc * v);  // =  S1(sc) - S2(tc)
+  assert(dP == (P1-P2));
+
+  d = dP.squared_length();   // return the closest distance
+  return parallel;
+}

include/igl/copyleft/cgal/segment_segment_squared_distance.h

@@ -0,0 +1,39 @@
+#ifndef IGL_COPYLEFT_CGAL_SEGMENT_SEGMENT_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_SEGMENT_SEGMENT_SQUARED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <CGAL/Segment_3.h>
+#include <CGAL/Point_3.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given two segments S1 and S2 find the points on each of closest
+      // approach and the squared distance thereof.
+      // 
+      // Inputs:
+      //   S1  first segment
+      //   S2  second segment
+      // Outputs:
+      //   P1  point on S1 closest to S2
+      //   P2  point on S2 closest to S1
+      //   d  distance betwee P1 and S2
+      // Returns true if the closest approach is unique.
+      template < typename Kernel>
+      IGL_INLINE bool segment_segment_squared_distance(
+          const CGAL::Segment_3<Kernel> & S1,
+          const CGAL::Segment_3<Kernel> & S2,
+          CGAL::Point_3<Kernel> & P1,
+          CGAL::Point_3<Kernel> & P2,
+          typename Kernel::FT & d
+          );
+
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "segment_segment_squared_distance.cpp"
+#endif
+
+#endif

include/igl/copyleft/cgal/snap_rounding.cpp

@@ -0,0 +1,199 @@
+#include "snap_rounding.h"
+#include "resolve_intersections.h"
+#include "subdivide_segments.h"
+#include <igl/remove_unreferenced.h>
+#include <igl/unique.h>
+#include <CGAL/Segment_2.h>
+#include <CGAL/Point_2.h>
+#include <CGAL/Vector_2.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <algorithm>
+
+template <
+  typename DerivedV, 
+  typename DerivedE, 
+  typename DerivedVI, 
+  typename DerivedEI,
+  typename DerivedJ>
+IGL_INLINE void igl::copyleft::cgal::snap_rounding(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedVI> & VI,
+  Eigen::PlainObjectBase<DerivedEI> & EI,
+  Eigen::PlainObjectBase<DerivedJ> & J)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace igl::copyleft::cgal;
+  using namespace std;
+  // Exact scalar type
+  typedef CGAL::Epeck Kernel;
+  typedef Kernel::FT EScalar;
+  typedef CGAL::Segment_2<Kernel> Segment_2;
+  typedef CGAL::Point_2<Kernel> Point_2;
+  typedef CGAL::Vector_2<Kernel> Vector_2;
+  typedef Matrix<EScalar,Dynamic,Dynamic>  MatrixXE;
+  // Convert vertex positions to exact kernel
+
+  MatrixXE VE;
+  {
+    VectorXi IM;
+    resolve_intersections(V,E,VE,EI,J,IM);
+    for_each(
+      EI.data(),
+      EI.data()+EI.size(),
+      [&IM](typename DerivedEI::Scalar& i){i=IM(i);});
+    VectorXi _;
+    remove_unreferenced( MatrixXE(VE), DerivedEI(EI), VE,EI,_);
+  }
+
+  // find all hot pixels
+  //// southwest and north east corners
+  //const RowVector2i SW(
+  //  round(VE.col(0).minCoeff()),
+  //  round(VE.col(1).minCoeff()));
+  //const RowVector2i NE(
+  //  round(VE.col(0).maxCoeff()),
+  //  round(VE.col(1).maxCoeff()));
+
+  // https://github.com/CGAL/cgal/issues/548
+  // Round an exact scalar to the nearest integer. A priori can't just round
+  // double. Suppose e=0.5+ε but double(e)<0.5
+  //
+  // Inputs:
+  //   e  exact number
+  // Outputs:
+  //   i  closest integer
+  const auto & round = [](const EScalar & e)->int
+  {
+    const double d = CGAL::to_double(e);
+    // get an integer that's near the closest int
+    int i = std::round(d);
+    EScalar di_sqr = CGAL::square((e-EScalar(i)));
+    const auto & search = [&i,&di_sqr,&e](const int dir)
+    {
+      while(true)
+      {
+        const int j = i+dir;
+        const EScalar dj_sqr = CGAL::square((e-EScalar(j)));
+        if(dj_sqr < di_sqr)
+        {
+          i = j;
+          di_sqr = dj_sqr;
+        }else
+        {
+          break;
+        }
+      }
+    };
+    // Try to increase/decrease int
+    search(1);
+    search(-1);
+    return i;
+  };
+  vector<Point_2> hot;
+  for(int i = 0;i<VE.rows();i++)
+  {
+    hot.emplace_back(round(VE(i,0)),round(VE(i,1)));
+  }
+  {
+    std::vector<size_t> _1,_2;
+    igl::unique(vector<Point_2>(hot),hot,_1,_2);
+  }
+
+  // find all segments intersecting hot pixels
+  //   split edge at closest point to hot pixel center
+  vector<vector<Point_2>>  steiner(EI.rows());
+  // initialize each segment with endpoints
+  for(int i = 0;i<EI.rows();i++)
+  {
+    steiner[i].emplace_back(VE(EI(i,0),0),VE(EI(i,0),1));
+    steiner[i].emplace_back(VE(EI(i,1),0),VE(EI(i,1),1));
+  }
+  // silly O(n²) implementation
+  for(const Point_2 & h : hot)
+  {
+    // North, East, South, West
+    Segment_2 wall[4] = 
+    {
+      {h+Vector_2(-0.5, 0.5),h+Vector_2( 0.5, 0.5)},
+      {h+Vector_2( 0.5, 0.5),h+Vector_2( 0.5,-0.5)},
+      {h+Vector_2( 0.5,-0.5),h+Vector_2(-0.5,-0.5)},
+      {h+Vector_2(-0.5,-0.5),h+Vector_2(-0.5, 0.5)}
+    };
+    // consider all segments
+    for(int i = 0;i<EI.rows();i++)
+    {
+      // endpoints
+      const Point_2 s(VE(EI(i,0),0),VE(EI(i,0),1));
+      const Point_2 d(VE(EI(i,1),0),VE(EI(i,1),1));
+      // if either end-point is in h's pixel then ignore
+      const Point_2 rs(round(s.x()),round(s.y()));
+      const Point_2 rd(round(d.x()),round(d.y()));
+      if(h == rs || h == rd)
+      {
+        continue;
+      }
+      // otherwise check for intersections with walls consider all walls
+      const Segment_2 si(s,d);
+      vector<Point_2> hits;
+      for(int j = 0;j<4;j++)
+      {
+        const Segment_2 & sj = wall[j];
+        if(CGAL::do_intersect(si,sj))
+        {
+          CGAL::Object result = CGAL::intersection(si,sj);
+          if(const Point_2 * p = CGAL::object_cast<Point_2 >(&result))
+          {
+            hits.push_back(*p);
+          }else if(const Segment_2 * s = CGAL::object_cast<Segment_2 >(&result))
+          {
+            // add both endpoints
+            hits.push_back(s->vertex(0));
+            hits.push_back(s->vertex(1));
+          }
+        }
+      }
+      if(hits.size() == 0)
+      {
+        continue;
+      }
+      // centroid of hits
+      Vector_2 cen(0,0);
+      for(const Point_2 & hit : hits)
+      {
+        cen = Vector_2(cen.x()+hit.x(), cen.y()+hit.y());
+      }
+      cen = Vector_2(cen.x()/EScalar(hits.size()),cen.y()/EScalar(hits.size()));
+      const Point_2 rcen(round(cen.x()),round(cen.y()));
+      // after all of that, don't add as a steiner unless it's going to round
+      // to h
+      if(rcen == h)
+      {
+        steiner[i].emplace_back(cen.x(),cen.y());
+      }
+    }
+  }
+  {
+    DerivedJ prevJ = J;
+    VectorXi IM;
+    subdivide_segments(MatrixXE(VE),MatrixXi(EI),steiner,VE,EI,J,IM);
+    for_each(J.data(),J.data()+J.size(),[&prevJ](typename DerivedJ::Scalar & j){j=prevJ(j);});
+    for_each(
+      EI.data(),
+      EI.data()+EI.size(),
+      [&IM](typename DerivedEI::Scalar& i){i=IM(i);});
+    VectorXi _;
+    remove_unreferenced( MatrixXE(VE), DerivedEI(EI), VE,EI,_);
+  }
+
+
+  VI.resize(VE.rows(),VE.cols());
+  for(int i = 0;i<VE.rows();i++)
+  {
+    for(int j = 0;j<VE.cols();j++)
+    {
+      VI(i,j) = round(CGAL::to_double(VE(i,j)));
+    }
+  }
+}

include/igl/copyleft/cgal/snap_rounding.h

@@ -0,0 +1,44 @@
+#ifndef IGL_COPYLEFT_CGAL_SNAP_ROUNDING_H
+#define IGL_COPYLEFT_CGAL_SNAP_ROUNDING_H
+
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // SNAP_ROUNDING Snap a list of possible intersecting segments with
+      // endpoints in any precision to _the_ integer grid.
+      //
+      // Inputs:
+      //   V  #V by 2 list of vertex positions
+      //   E  #E by 2 list of segment indices into V
+      // Outputs:
+      //   VI  #VI by 2 list of output integer vertex positions, rounded copies
+      //     of V are always the first #V vertices
+      //   EI  #EI by 2 list of segment indices into V, #EI ≥ #E
+      //   J  #EI list of indices into E revealing "parent segments"
+      template <
+        typename DerivedV, 
+        typename DerivedE, 
+        typename DerivedVI, 
+        typename DerivedEI,
+        typename DerivedJ>
+      IGL_INLINE void snap_rounding(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedE> & E,
+        Eigen::PlainObjectBase<DerivedVI> & VI,
+        Eigen::PlainObjectBase<DerivedEI> & EI,
+        Eigen::PlainObjectBase<DerivedJ> & J);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "snap_rounding.cpp"
+#endif
+
+#endif

include/igl/copyleft/cgal/subdivide_segments.cpp

@@ -0,0 +1,127 @@
+#include "subdivide_segments.h"
+#include "row_to_point.h"
+#include <igl/unique.h>
+#include <igl/list_to_matrix.h>
+#include <igl/copyleft/cgal/assign_scalar.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/Segment_2.h>
+#include <CGAL/Point_2.h>
+#include <algorithm>
+#include <vector>
+
+template <
+  typename DerivedV, 
+  typename DerivedE,
+  typename Kernel, 
+  typename DerivedVI, 
+  typename DerivedEI,
+  typename DerivedJ,
+  typename DerivedIM>
+IGL_INLINE void igl::copyleft::cgal::subdivide_segments(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedE> & E,
+  const std::vector<std::vector<CGAL::Point_2<Kernel> > > & _steiner,
+  Eigen::PlainObjectBase<DerivedVI> & VI,
+  Eigen::PlainObjectBase<DerivedEI> & EI,
+  Eigen::PlainObjectBase<DerivedJ> & J,
+  Eigen::PlainObjectBase<DerivedIM> & IM)
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+
+  // Exact scalar type
+  typedef Kernel K;
+  typedef typename Kernel::FT EScalar;
+  typedef CGAL::Segment_2<Kernel> Segment_2;
+  typedef CGAL::Point_2<Kernel> Point_2;
+  typedef Matrix<EScalar,Dynamic,Dynamic>  MatrixXE;
+
+  // non-const copy
+  std::vector<std::vector<CGAL::Point_2<Kernel> > > steiner = _steiner;
+
+  // Convert vertex positions to exact kernel
+  MatrixXE VE(V.rows(),V.cols());
+  for(int i = 0;i<V.rows();i++)
+  {
+    for(int j = 0;j<V.cols();j++)
+    {
+      VE(i,j) = V(i,j);
+    }
+  }
+
+  // number of original vertices
+  const int n = V.rows();
+  // number of original segments
+  const int m = E.rows();
+  // now steiner contains lists of points (unsorted) for each edge. Sort them
+  // and count total number of vertices and edges
+  int ni = 0;
+  int mi = 0;
+  // new steiner points
+  std::vector<Point_2> S;
+  std::vector<std::vector<typename DerivedE::Scalar> > vEI;
+  std::vector<typename DerivedJ::Scalar> vJ;
+  for(int i = 0;i<m;i++)
+  {
+    {
+      const Point_2 s = row_to_point<K>(VE,E(i,0));
+      std::sort(
+        steiner[i].begin(),
+        steiner[i].end(),
+        [&s](const Point_2 & A, const Point_2 & B)->bool
+        {
+          return (A-s).squared_length() < (B-s).squared_length();
+        });
+    }
+    // remove duplicates
+    steiner[i].erase(
+      std::unique(steiner[i].begin(), steiner[i].end()), 
+      steiner[i].end());
+    {
+      int s = E(i,0);
+      // legs to each steiner in order
+      for(int j = 1;j<steiner[i].size()-1;j++)
+      {
+        int d = n+S.size();
+        S.push_back(steiner[i][j]);
+        vEI.push_back({s,d});
+        vJ.push_back(i);
+        s = d;
+      }
+      // don't forget last (which might only) leg
+      vEI.push_back({s,E(i,1)});
+      vJ.push_back(i);
+    }
+  }
+  // potentially unnecessary copying ...
+  VI.resize(n+S.size(),2);
+  for(int i = 0;i<V.rows();i++)
+  {
+    for(int j = 0;j<V.cols();j++)
+    {
+      assign_scalar(V(i,j),VI(i,j));
+    }
+  }
+  for(int i = 0;i<S.size();i++)
+  {
+    assign_scalar(S[i].x(),VI(n+i,0));
+    assign_scalar(S[i].y(),VI(n+i,1));
+  }
+  list_to_matrix(vEI,EI);
+  list_to_matrix(vJ,J);
+  {
+    // Find unique mapping
+    std::vector<Point_2> vVES,_1;
+    for(int i = 0;i<n;i++)
+    {
+      vVES.push_back(row_to_point<K>(VE,i));
+    }
+    vVES.insert(vVES.end(),S.begin(),S.end());
+    std::vector<size_t> vA,vIM;
+    igl::unique(vVES,_1,vA,vIM);
+    // Push indices back into vVES
+    for_each(vIM.data(),vIM.data()+vIM.size(),[&vA](size_t & i){i=vA[i];});
+    list_to_matrix(vIM,IM);
+  }
+}

include/igl/copyleft/cgal/subdivide_segments.h

@@ -0,0 +1,49 @@
+#ifndef IGL_COPYLEFT_CGAL_SUBDIVIDE_SEGMENTS_H
+#define IGL_COPYLEFT_CGAL_SUBDIVIDE_SEGMENTS_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include <CGAL/Segment_2.h>
+#include <CGAL/Point_2.h>
+#include <vector>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Insert steiner points to subdivide a given set of line segments
+      // 
+      // Inputs:
+      //   V  #V by 2 list of vertex positions
+      //   E  #E by 2 list of segment indices into V
+      //   steiner  #E list of lists of unsorted steiner points (including
+      //     endpoints) along the #E original segments
+      // Outputs:
+      //   VI  #VI by 2 list of output vertex positions, copies of V are always
+      //     the first #V vertices
+      //   EI  #EI by 2 list of segment indices into V, #EI ≥ #E
+      //   J  #EI list of indices into E revealing "parent segments"
+      //   IM  #VI list of indices into VV of unique vertices.
+      template <
+        typename DerivedV, 
+        typename DerivedE,
+        typename Kernel, 
+        typename DerivedVI, 
+        typename DerivedEI,
+        typename DerivedJ,
+        typename DerivedIM>
+      IGL_INLINE void subdivide_segments(
+        const Eigen::PlainObjectBase<DerivedV> & V,
+        const Eigen::PlainObjectBase<DerivedE> & E,
+        const std::vector<std::vector<CGAL::Point_2<Kernel> > > & steiner,
+        Eigen::PlainObjectBase<DerivedVI> & VI,
+        Eigen::PlainObjectBase<DerivedEI> & EI,
+        Eigen::PlainObjectBase<DerivedJ> & J,
+        Eigen::PlainObjectBase<DerivedIM> & IM);
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "subdivide_segments.cpp"
+#endif
+#endif

include/igl/copyleft/cgal/triangle_triangle_squared_distance.cpp

@@ -0,0 +1,104 @@
+#include "triangle_triangle_squared_distance.h"
+#include "point_triangle_squared_distance.h"
+
+template < typename Kernel>
+IGL_INLINE bool igl::copyleft:cgal::triangle_triangle_squared_distance(
+  const CGAL::Triangle_3<Kernel> & T1,
+  const CGAL::Triangle_3<Kernel> & T2,
+  CGAL::Point_3<Kernel> & P1,
+  CGAL::Point_3<Kernel> & P2,
+  typename Kernel::FT & d)
+{
+  typedef CGAL::Point_3<Kernel> Point_3;
+  typedef CGAL::Vector_3<Kernel> Vector_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef CGAL::Segment_3<Kernel> Segment_3;
+  typedef Kernel::FT EScalar;
+  assert(!T1.is_degenerate());
+  assert(!T2.is_degenerate());
+
+  bool unique = true;
+  if(CGAL::do_intersect(T1,T2))
+  {
+    // intersecting triangles have zero (squared) distance
+    CGAL::Object result = CGAL::intersection(T1,T2);
+    // Some point on the intersection result
+    CGAL::Point_3<Kernel> Q;
+    if(const Point_3 * p = CGAL::object_cast<Point_3 >(&result))
+    {
+      Q = *p;
+    }else if(const Segment_3 * s = CGAL::object_cast<Segment_3 >(&result))
+    {
+      unique = false;
+      Q = s->source();
+    }else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj))
+    {
+      Q = s->vertex(0);
+      unique = false;
+    }else if(const std::vector<Point_3 > *polyp = 
+      CGAL::object_cast< std::vector<Point_3 > >(&obj))
+    {
+      assert(polyp->size() > 0 && "intersection poly should not be empty");
+      Q = polyp[0];
+      unique = false;
+    }else
+    {
+      assert(false && "Unknown intersection result");
+    }
+    P1 = Q;
+    P2 = Q;
+    d = 0;
+    return unique;
+  }
+  // triangles do not intersect: the points of closest approach must be on the
+  // boundary of one of the triangles
+  d = std::numeric_limits<double>::infinity();
+  const auto & vertices_face = [](
+    const Triangle_3 & T1,
+    const Triangle_3 & T2,
+    Point_3 & P1,
+    Point_3 & P2,
+    EScalar & d)
+  {
+    for(int i = 0;i<3;i++)
+    {
+      const Point_3 vi = T1.vertex(i);
+      Point_3 P2i;
+      EScalar di;
+      point_triangle_squared_distance(vi,T2,P2i,di);
+      if(di<d)
+      {
+        d = di;
+        P1 = vi;
+        P2 = P2i;
+        unique = true;
+      }else if(d == di)
+      {
+        // edge of T1 floating parallel above T2
+        unique = false;
+      }
+    }
+  };
+  vertices_face(T1,T2,P1,P2,d);
+  vertices_face(T2,T1,P1,P2,d);
+  for(int i = 0;i<3;i++)
+  {
+    const Segment_3<Kernel> s1( T1.vertex(i+1), T1.vertex(i+2));
+    for(int j = 0;j<3;j++)
+    {
+      const Segment_3<Kernel> s2( T2.vertex(i+1), T2.vertex(i+2));
+      Point_3 P1ij;
+      Point_3 P2ij;
+      EScalar dij;
+      bool uniqueij = segment_segment_squared_distance(s1,s2,P1ij,P2ij,dij);
+      if(dij < d)
+      {
+        P1 = P1ij;
+        P2 = P2ij;
+        d = dij;
+        unique = uniqueij;
+      }
+    }
+  }
+  return unique;
+}

include/igl/copyleft/cgal/triangle_triangle_squared_distance.h

@@ -0,0 +1,38 @@
+#ifndef IGL_COPYLEFT_CGAL_TRIANGLE_TRIANGLE_SQUARED_DISTANCE_H
+#define IGL_COPYLEFT_CGAL_TRIANGLE_TRIANGLE_SQUARED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <CGAL/Triangle_3.h>
+#include <CGAL/Point_3.h>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Given two triangles T1 and T2 find the points on each of closest
+      // approach and the squared distance thereof.
+      // 
+      // Inputs:
+      //   T1  first triangle
+      //   T2  second triangle
+      // Outputs:
+      //   P1  point on T1 closest to T2
+      //   P2  point on T2 closest to T1
+      //   d  distance betwee P1 and T2
+      // Returns true if the closest approach is unique.
+      template < typename Kernel>
+      IGL_INLINE bool triangle_triangle_squared_distance(
+        const CGAL::Triangle_3<Kernel> & T1,
+        const CGAL::Triangle_3<Kernel> & T2,
+        CGAL::Point_3<Kernel> & P1,
+        CGAL::Point_3<Kernel> & P2,
+        typename Kernel::FT & d);
+    }
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "triangle_triangle_squared_distance.cpp"
+#endif
+
+#endif
+