ear clipping triangulation and robust 2d segment intersection detect

include/igl/predicates/ear_clipping.cpp

@@ -0,0 +1,113 @@
+
+#include <igl/slice.h>
+#include "ear_clipping.h"
+#include "point_inside_polygon.h"
+#include "predicates.h"
+
+template <typename DerivedP>
+IGL_INLINE bool igl::predicates::is_ear(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::VectorXi& RT,
+  const Eigen::VectorXi& L,
+  const Eigen::VectorXi& R,
+  const int i
+){
+  typedef typename DerivedP::Scalar Scalar;
+
+  int a = L(i), b = R(i);
+  if(RT(i) != 0 || RT(a) != 0 || RT(b) != 0) return false;
+  Eigen::RowVector2d pa = P.row(a);
+  Eigen::RowVector2d pb = P.row(b);
+  Eigen::RowVector2d 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);
+  int k=R(b);
+  while(k!=a){
+    Eigen::MatrixX2d T(3,2);
+    T<<P.row(a),P.row(i),P.row(b);
+    Eigen::RowVector2d q=P.row(k);
+    if(point_inside_polygon(T,q)) 
+      return false;
+    k=R(k);
+  }
+  return true;
+}
+
+// https://www.geometrictools.com/Documentation/TriangulationByEarClipping.pdf
+template <typename DerivedP>
+IGL_INLINE void igl::predicates::ear_clipping(
+  const Eigen::MatrixBase<DerivedP>& P,
+  const Eigen::VectorXi& RT,
+  Eigen::VectorXi& I,
+  Eigen::MatrixXi& eF, 
+  Eigen::PlainObjectBase<DerivedP>& nP
+){
+
+  Eigen::VectorXi L(P.rows());
+  Eigen::VectorXi R(P.rows());
+  for(int i=0;i<P.rows();i++){
+    L(i) = (i-1+P.rows())%P.rows();
+    R(i) = (i+1)%P.rows();
+  }
+
+  Eigen::VectorXi ears; // mark ears
+  Eigen::VectorXi 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
+    int e = 0;
+    while(e<ears.rows()&&ears(e)!=1) e++;
+    
+    // find valid neighbors
+    int 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());
+  int j=0;
+  for(int 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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&, 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,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 EAR_CLIPPING
+#define EAR_CLIPPING
+
+#include <Eigen/Core>
+#include "../igl_inline.h"
+
+namespace igl
+{
+  namespace predicates
+  {
+    template <typename DerivedP>
+    IGL_INLINE void ear_clipping(
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::VectorXi& RT,
+      Eigen::VectorXi& I,
+      Eigen::MatrixXi& eF, 
+      Eigen::PlainObjectBase<DerivedP>& nP
+    );
+
+    template <typename DerivedP>
+    IGL_INLINE bool is_ear(
+      const Eigen::MatrixBase<DerivedP>& P,
+      const Eigen::VectorXi& RT,
+      const Eigen::VectorXi& L,
+      const Eigen::VectorXi& R, 
+      const int i
+    );
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "ear_clipping.cpp"
+#endif
+
+
+#endif

include/igl/predicates/point_inside_polygon.cpp

@@ -0,0 +1,30 @@
+// 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_polygon.h"
+
+template <typename Scalar>
+bool igl::predicates::point_inside_polygon(
+    const Eigen::Matrix<Scalar,-1,2>& P,
+    const Eigen::Matrix<Scalar,1,2>& q
+){
+  for(int i=0;i<P.rows();i++){
+    int i_1 = (i+1) % P.rows();
+    Eigen::RowVector2d a = P.row(i);
+    Eigen::RowVector2d 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_polygon<double>(Eigen::Matrix<double, -1, 2, 0, -1, 2> const&, Eigen::Matrix<double, 1, 2, 1, 1, 2> const&);
+#endif

include/igl/predicates/point_inside_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 PREDICATE_POINT_INSIDE_POLYGON_H
+#define PREDICATE_POINT_INSIDE_POLYGON_H
+
+
+
+#include "../igl_inline.h"
+#include <Eigen/Core>
+#include "predicates.h"
+
+namespace igl
+{
+  namespace predicates
+  {
+      // check whether 2d point lies inside 2d polygon
+      // Inputs:
+      //   P: n*2 polygon, n >= 3
+      //   q: 2d query point
+      // Returns true if point is inside polygon
+      template <typename Scalar>
+      bool point_inside_polygon(
+          const Eigen::Matrix<Scalar,-1,2>& P,
+          const Eigen::Matrix<Scalar,1,2>& q
+      );
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_inside_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,46 @@
+// 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 <igl/predicates/predicates.h>
+namespace igl
+{
+namespace predicates
+{
+  
+    // Given two segments in 2d test whether they intersect each other
+    // using 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
+    //   eps: precision used for colinear case
+    // 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);
+  
+}