Merge pull request #1169 from hankstag/pr-predicates

Ear clipping triangulation and robust 2d segment intersection test

include/igl/predicates/ear_clipping.cpp

@@ -0,0 +1,124 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Hanxiao Shen <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <igl/slice.h>
+#include "ear_clipping.h"
+#include "point_inside_convex_polygon.h"
+#include "predicates.h"
+
+template <typename DerivedP, typename DerivedRT,
+          typename DerivedF, typename DerivedI>
+IGL_INLINE void igl::predicates::ear_clipping(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::MatrixBase<DerivedRT>& RT,
+  Eigen::PlainObjectBase<DerivedI>& I,
+  Eigen::PlainObjectBase<DerivedF>& eF,
+  Eigen::PlainObjectBase<DerivedP>& nP
+){
+  typedef typename DerivedF::Scalar Index;
+  typedef typename DerivedP::Scalar Scalar;
+  static_assert(std::is_same<typename DerivedI::Scalar,
+                             typename DerivedF::Scalar>::value,
+                "index type should be consistent");
+  
+  // check whether vertex i is an ear
+  auto is_ear = [](
+    const Eigen::MatrixBase<DerivedP>& P,
+    const Eigen::MatrixBase<DerivedRT>& RT,
+    const Eigen::Matrix<Index,Eigen::Dynamic,1>& L,
+    const Eigen::Matrix<Index,Eigen::Dynamic,1>& R,
+    const Index i
+  ){
+    
+    Index a = L(i), b = R(i);
+    if(RT(i) != 0 || RT(a) != 0 || RT(b) != 0) return false;
+    Eigen::Matrix<Scalar,1,2> pa = P.row(a);
+    Eigen::Matrix<Scalar,1,2> pb = P.row(b);
+    Eigen::Matrix<Scalar,1,2> pi = P.row(i);
+    auto r = igl::predicates::orient2d(pa,pi,pb);
+    if(r == igl::predicates::Orientation::NEGATIVE || 
+       r == igl::predicates::Orientation::COLLINEAR) return false;
+    
+    // check if any vertex is lying inside triangle (a,b,i);
+    Index k=R(b);
+    while(k!=a){
+      Eigen::Matrix<Scalar,-1,2> T(3,2);
+      T<<P.row(a),P.row(i),P.row(b);
+      Eigen::Matrix<Scalar,1,2> q=P.row(k);
+      if(igl::predicates::point_inside_convex_polygon(T,q))
+        return false;
+      k=R(k);
+    }
+    return true;
+  };
+
+  Eigen::Matrix<Index,Eigen::Dynamic,1> L(P.rows());
+  Eigen::Matrix<Index,Eigen::Dynamic,1> R(P.rows());
+  for(int i=0;i<P.rows();i++){
+    L(i) = Index((i-1+P.rows())%P.rows());
+    R(i) = Index((i+1)%P.rows());
+  }
+
+  Eigen::Matrix<Index,Eigen::Dynamic,1> ears; // mark ears
+  Eigen::Matrix<Index,Eigen::Dynamic,1> X; // clipped vertices
+  ears.setZero(P.rows());
+  X.setZero(P.rows());
+
+  // initialize ears
+  for(int i=0;i<P.rows();i++){
+    ears(i) = is_ear(P,RT,L,R,i);
+  }
+
+  // clip ears until none left
+  while(ears.maxCoeff()==1){
+    
+    // find the first ear
+    Index e = 0;
+    while(e<ears.rows()&&ears(e)!=1) e++;
+    
+    // find valid neighbors
+    Index a = L(e), b = R(e);
+    if(a == b) break;
+
+    // clip ear and store face
+    eF.conservativeResize(eF.rows()+1,3);
+    eF.bottomRows(1)<<a,e,b;
+    L(b) = a;
+    L(e) = -1;
+    R(a) = b;
+    R(e) = -1;
+    ears(e) = 0; // mark vertex e as non-ear
+
+    // update neighbor's ear status
+    ears(a) = is_ear(P,RT,L,R,a);
+    ears(b) = is_ear(P,RT,L,R,b);
+    X(e) = 1;
+
+    // when only one edge left
+    // mark the endpoints as clipped
+    if(L(a)==b && R(b)==a){
+      X(a) = 1;
+      X(b) = 1;
+    }
+  }
+  
+  // collect remaining vertices if theres any
+  for(int i=0;i<X.rows();i++)
+    X(i) = 1-X(i);
+  I.resize(X.sum());
+  Index j=0;
+  for(Index i=0;i<X.rows();i++)
+    if(X(i)==1){
+      I(j++) = i;
+    }
+  igl::slice(P,I,1,nP);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::predicates::ear_clipping<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/predicates/ear_clipping.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Hanxiao Shen <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_PREDICATES_EAR_CLIPPING_H
+#define IGL_PREDICATES_EAR_CLIPPING_H
+
+#include <Eigen/Core>
+#include "../igl_inline.h"
+
+namespace igl
+{
+  namespace predicates
+  {
+
+    // Implementation of ear clipping triangulation algorithm for a 2D polygon.
+    // https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
+    // If the polygon is simple, all vertices will be clipped and the result mesh is (P,eF)
+    // Otherwise, the function will try to clip as many ears as possible.
+    //
+    // Input:
+    // P : n*2, size n 2D polygon
+    // RT: n*1, preserved vertices (do not clip) marked as 1, otherwise 0
+    // Output:
+    // I : size #nP vector, maps index from nP to P, e.g. nP's ith vertex is origianlly I(i) in P
+    // eF: clipped ears, in original index of P
+    // nP: leftover vertices after clipping
+
+    template <typename DerivedP, typename DerivedRT,
+              typename DerivedF, typename DerivedI>
+    IGL_INLINE void ear_clipping(
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::MatrixBase<DerivedRT>& RT,
+      Eigen::PlainObjectBase<DerivedI>& I,
+      Eigen::PlainObjectBase<DerivedF>& eF, 
+      Eigen::PlainObjectBase<DerivedP>& nP
+    );
+
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "ear_clipping.cpp"
+#endif
+
+
+#endif

include/igl/predicates/point_inside_convex_polygon.cpp

@@ -0,0 +1,33 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Hanxiao Shen <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "point_inside_convex_polygon.h"
+
+template <typename DerivedP, typename DerivedQ>
+IGL_INLINE bool igl::predicates::point_inside_convex_polygon(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::MatrixBase<DerivedQ>& q
+){
+  EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(DerivedP, Eigen::Dynamic, 2);
+  EIGEN_STATIC_ASSERT_MATRIX_SPECIFIC_SIZE(DerivedQ, 1, 2);
+  typedef typename DerivedP::Scalar Scalar;
+  for(int i=0;i<P.rows();i++){
+    int i_1 = (i+1) % P.rows();
+    Eigen::Matrix<Scalar,1,2> a = P.row(i);
+    Eigen::Matrix<Scalar,1,2> b = P.row(i_1);
+    auto r = igl::predicates::orient2d(a,b,q);
+    if(r == igl::predicates::Orientation::COLLINEAR || 
+       r == igl::predicates::Orientation::NEGATIVE)
+      return false;
+  }
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::predicates::point_inside_convex_polygon<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&);
+#endif

include/igl/predicates/point_inside_convex_polygon.h

@@ -0,0 +1,39 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Hanxiao Shen <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_PREDICATES_POINT_INSIDE_CONVEX_POLYGON_H
+#define IGL_PREDICATES_POINT_INSIDE_CONVEX_POLYGON_H
+
+
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "predicates.h"
+
+namespace igl
+{
+  namespace predicates
+  {
+    // check whether 2d point lies inside 2d convex polygon
+    // Inputs:
+    //   P: n*2 polygon, n >= 3
+    //   q: 2d query point
+    // Returns true if point is inside polygon
+    template <typename DerivedP, typename DerivedQ>
+    IGL_INLINE bool point_inside_convex_polygon(
+        const Eigen::MatrixBase<DerivedP>& P,
+        const Eigen::MatrixBase<DerivedQ>& q
+    );
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_inside_convex_polygon.cpp"
+#endif
+
+#endif

include/igl/predicates/segment_segment_intersect.cpp

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Hanxiao Shen <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "segment_segment_intersect.h"
+
+// https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
+template <typename DerivedP>
+IGL_INLINE bool igl::predicates::segment_segment_intersect(
+  const Eigen::MatrixBase<DerivedP>& a,
+  const Eigen::MatrixBase<DerivedP>& b,
+  const Eigen::MatrixBase<DerivedP>& c,
+  const Eigen::MatrixBase<DerivedP>& d
+)
+{
+  typename DerivedP::Scalar Scalar;
+  
+  auto t1 = igl::predicates::orient2d(a,b,c);
+  auto t2 = igl::predicates::orient2d(b,c,d);
+  auto t3 = igl::predicates::orient2d(a,b,d);
+  auto t4 = igl::predicates::orient2d(a,c,d);
+
+  // assume m,n,p are colinear, check whether p is in range [m,n]
+  auto on_segment = [](
+    const Eigen::MatrixBase<DerivedP>& m,
+    const Eigen::MatrixBase<DerivedP>& n,
+    const Eigen::MatrixBase<DerivedP>& p
+  ){
+      return ((p(0) >= std::min(m(0),n(0))) &&
+              (p(0) <= std::max(m(0),n(0))) &&
+              (p(1) >= std::min(m(1),n(1))) &&
+              (p(1) <= std::max(m(1),n(1))));
+  };
+  
+  // colinear case
+  if((t1 == igl::predicates::Orientation::COLLINEAR && on_segment(a,b,c)) ||
+     (t2 == igl::predicates::Orientation::COLLINEAR && on_segment(c,d,b)) ||
+     (t3 == igl::predicates::Orientation::COLLINEAR && on_segment(a,b,d)) ||
+     (t4 == igl::predicates::Orientation::COLLINEAR && on_segment(c,d,a))) 
+     return true;
+  
+  // ordinary case
+  return (t1 != t3 && t2 != t4);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::predicates::segment_segment_intersect<Eigen::Matrix<double, 1, 2, 1, 1, 2> >(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&, Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2> > const&);
+#endif

include/igl/predicates/segment_segment_intersect.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2019 Hanxiao Shen <[email protected]>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef IGL_PREDICATES_SEGMENT_SEGMENT_INTERSECT_H
+#define IGL_PREDICATES_SEGMENT_SEGMENT_INTERSECT_H
+
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include "predicates.h"
+namespace igl
+{
+  namespace predicates
+  {
+  
+    // Given two segments in 2d test whether they intersect each other
+    // using predicates orient2d
+    // 
+    // Inputs:
+    //   A:   1st endpoint of segment 1
+    //   B:   2st endpoint of segment 1
+    //   C:   1st endpoint of segment 2
+    //   D:   2st endpoint of segment 2
+    // Returns true if they intersect
+    
+    template <typename DerivedP>
+    IGL_INLINE bool segment_segment_intersect(
+      // input:
+      const Eigen::MatrixBase<DerivedP>& A,
+      const Eigen::MatrixBase<DerivedP>& B,
+      const Eigen::MatrixBase<DerivedP>& C,
+      const Eigen::MatrixBase<DerivedP>& D
+    );
+
+  }
+}
+#ifndef IGL_STATIC_LIBRARY
+#   include "segment_segment_intersect.cpp"
+#endif
+#endif //IGL_PREDICATES_SEGMENT_SEGMENT_INTERSECT_H

tests/include/igl/predicates/ear_clipping.cpp

@@ -0,0 +1,16 @@
+#include <test_common.h>
+#include <igl/predicates/ear_clipping.h>
+
+TEST_CASE("ear_clipping: boolean", "[igl/predicates]")
+{
+  // Example1: simple polygon 
+  Eigen::MatrixXd polygon(10,2);
+  polygon<<2,-3,4,1,5.5,-2,6,2.5,5,1,4,5,3,0,1,1,1,5,0,0;
+  Eigen::VectorXi RT,nR,M;
+  Eigen::MatrixXi eF;
+  Eigen::MatrixXd nP;
+  RT.setZero(polygon.rows());
+  igl::predicates::ear_clipping(polygon,RT,M,eF,nP);
+  REQUIRE(nP.rows() == 0);
+
+}

tests/include/igl/predicates/segment_segment_intersect.cpp

@@ -0,0 +1,43 @@
+#include <test_common.h>
+#include <igl/predicates/segment_segment_intersect.h>
+#include <iomanip>
+
+TEST_CASE("segment_segment_intersect: robust", "[igl/predicates]")
+{
+  // example 1: vanila intersecting case
+  auto A1 = Eigen::RowVector2d(0, 128.5);
+  auto B1 = Eigen::RowVector2d(-77.44,1.2);
+  auto C1 = Eigen::RowVector2d(-83.2,2.8);
+  auto D1 = Eigen::RowVector2d(-1.0,-1.0);
+
+  bool check1 = igl::predicates::segment_segment_intersect(A1,B1,C1,D1);
+  REQUIRE(check1 == true);
+  
+  // example 2: colinear overlapping
+  auto A2 = Eigen::RowVector2d(1.0,5.0);
+  auto B2 = Eigen::RowVector2d(1.0,9.0);
+  auto C2 = Eigen::RowVector2d(1.0,8.0);
+  auto D2 = Eigen::RowVector2d(1.0,12.0);
+  
+  bool check2 = igl::predicates::segment_segment_intersect(A2,B2,C2,D2);
+  REQUIRE(check2 == true);
+
+  // example 3: colinear touching endpoint
+  auto A3 = Eigen::RowVector2d(0.0,0.0);
+  auto B3 = Eigen::RowVector2d(1.5,1.5);
+  auto C3 = Eigen::RowVector2d(1.5,1.5);
+  auto D3 = Eigen::RowVector2d(2.0,2.0);
+  
+  bool check3 = igl::predicates::segment_segment_intersect(A3,B3,C3,D3);
+  REQUIRE(check3 == true);
+
+  // example 6: colinear not touching endpoint
+  double eps = 1e-14;
+  auto A4 = Eigen::RowVector2d(0.0,0.0);
+  auto B4 = Eigen::RowVector2d(1.5,1.5);
+  auto C4 = Eigen::RowVector2d(1.5+eps,1.5+eps);
+  auto D4 = Eigen::RowVector2d(2.0,2.0);
+  bool check4 = igl::predicates::segment_segment_intersect(A4,B4,C4,D4);
+  REQUIRE(check4 == false);
+  
+}