Ear clipping function fix (#1565) [ci skip]

* fix the is_ear check

* fix order of outputs; overload that does any orientation

---------

Co-authored-by: Alec Jacobson <[email protected]>

include/igl/predicates/ear_clipping.cpp

@@ -6,20 +6,23 @@
 // 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 "../slice.h"
 #include "ear_clipping.h"
-#include "point_inside_convex_polygon.h"
 #include "predicates.h"
+#include "segment_segment_intersect.h"
+#include "../slice.h"
+#include "../turning_number.h"
 
-template <typename DerivedP, typename DerivedRT,
-          typename DerivedF, typename DerivedI>
+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
-){
+  Eigen::PlainObjectBase<DerivedF>& eF, 
+  Eigen::PlainObjectBase<DerivedI>& I)
+{
   typedef typename DerivedF::Scalar Index;
   typedef typename DerivedP::Scalar Scalar;
   static_assert(std::is_same<typename DerivedI::Scalar,
@@ -35,43 +38,75 @@ IGL_INLINE void igl::predicates::ear_clipping(
     const Index i
   ){
     
+    // check if any edge in the leftover polygon intersects triangle (a,b,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 || 
+    auto r = igl::predicates::orient2d(pa, pi, pb);
+    if(r == igl::predicates::Orientation::NEGATIVE ||
        r == igl::predicates::Orientation::COLLINEAR) return false;
+
+    // check edge (b, R(b))
+    Index k = b;
+    Index l = R(b);
+    Eigen::Matrix<Scalar, 1, 2> pk = P.row(k);
+    Eigen::Matrix<Scalar, 1, 2> pl = P.row(l);
+    auto r1 = igl::predicates::orient2d(pb, pl, pa);
+    auto r2 = igl::predicates::orient2d(pi, pl, pb);
+    if (l == a) return true;  // single triangle
+    if (r1 != igl::predicates::Orientation::POSITIVE &&
+        r2 != igl::predicates::Orientation::POSITIVE) {
+      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))
+    // iterate through all edges do not have common endpoints with a, b
+    k = l;
+    l = R(k);
+    while (l != a) {
+      pk = P.row(k);
+      pl = P.row(l);
+      if (igl::predicates::segment_segment_intersect(pa, pb, pk, pl) ||
+          igl::predicates::segment_segment_intersect(pa, pi, pk, pl) ||
+          igl::predicates::segment_segment_intersect(pi, pb, pk, pl)
+      ){
         return false;
-      k=R(k);
+      }
+      k = l;
+      l = R(k);
+    }
+
+    // check edge (L(a), a)
+    pk = P.row(k);
+    pl = P.row(l);
+    r1 = igl::predicates::orient2d(pb, pk, pl);
+    r2 = igl::predicates::orient2d(pi, pa, pk);
+    if (r1 != igl::predicates::Orientation::POSITIVE &&
+        r2 != igl::predicates::Orientation::POSITIVE
+    ){
+      return false;
     }
+
     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());
+  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
+  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);
+  for(int i = 0; i < P.rows(); i++){
+    ears(i) = is_ear(P, RT, L, R, i);
   }
 
   // clip ears until none left
@@ -79,7 +114,7 @@ IGL_INLINE void igl::predicates::ear_clipping(
     
     // find the first ear
     Index e = 0;
-    while(e<ears.rows()&&ears(e)!=1) e++;
+    while(e < ears.rows() && ears(e) != 1) e++;
     
     // find valid neighbors
     Index a = L(e), b = R(e);
@@ -87,7 +122,7 @@ IGL_INLINE void igl::predicates::ear_clipping(
 
     // clip ear and store face
     eF.conservativeResize(eF.rows()+1,3);
-    eF.bottomRows(1)<<a,e,b;
+    eF.bottomRows(1)<<a, e, b;
     L(b) = a;
     L(e) = -1;
     R(a) = b;
@@ -95,30 +130,57 @@ IGL_INLINE void igl::predicates::ear_clipping(
     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);
+    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){
+    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);
+  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){
+  for(Index i = 0;i < X.rows(); i++)
+  {
+    if(X(i) == 1)
+    {
       I(j++) = i;
     }
-  igl::slice(P,I,1,nP);
+  }
+}
+
+template <typename DerivedP,
+          typename DerivedF>
+IGL_INLINE bool igl::predicates::ear_clipping(
+  const Eigen::MatrixBase<DerivedP>& P,
+  Eigen::PlainObjectBase<DerivedF>& eF)
+{
+  if(turning_number(P) < 0)
+  {
+    // reverse input
+    const auto rP = P.colwise().reverse().eval();
+    const bool ret = igl::predicates::ear_clipping(rP, eF);
+    // flip output
+    eF = eF.rowwise().reverse();
+    return ret;
+  }
+  const Eigen::Matrix<typename DerivedP::Scalar, Eigen::Dynamic, 1> RT = 
+    Eigen::Matrix<typename DerivedP::Scalar, Eigen::Dynamic, 1>::Zero(P.rows());
+  Eigen::Matrix<typename DerivedF::Scalar, Eigen::Dynamic, 1> I;
+  igl::predicates::ear_clipping(P, RT, eF, I);
+  return I.size() == 0;
 }
 
 #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> >&);
+template bool igl::predicates::ear_clipping<Eigen::Matrix<double, -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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>>&);
+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>>&);
 #endif

include/igl/predicates/ear_clipping.h

@@ -18,25 +18,39 @@ namespace igl
   {
     /// 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.
+    /// If the polygon is simple and oriented counter-clockwise, all vertices
+    /// will be clipped and the result mesh is (P,eF) Otherwise, the function
+    /// will try to clip as many ears as possible.
     ///
-    /// @param[in] P : n*2, size n 2D polygon
+    /// @param[in] P : n*2, size n 2D polygon 
     /// @param[in] RT: n*1, preserved vertices (do not clip) marked as 1, otherwise 0
-    /// @param[out] I : size #nP vector, maps index from nP to P, e.g. nP's ith vertex is origianlly I(i) in P
     /// @param[out] eF: clipped ears, in original index of P
-    /// @param[out] nP: leftover vertices after clipping
+    /// @param[out] I : size #nP vector, maps index from nP to P, e.g. nP's ith vertex is origianlly I(i) in P
+    ///
+    /// \pre To result in a proper mesh, P should be oriented counter-clockwise
+    /// with no self-intersections.
+    ///
+    /// \note This implementation does not handle polygons with holes.
     ///
     /// \bug https://github.com/libigl/libigl/issues/1563
-    template <typename DerivedP, typename DerivedRT,
-              typename DerivedF, typename DerivedI>
+    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
-    );
+      Eigen::PlainObjectBase<DerivedI>& I);
+    /// \overload
+    /// \brief Reverses P if necessary. Orientation of output will match input
+    /// \return true if mesh is proper (should correspond to input being a
+    /// simple polygon in either orientation), false otherwise.
+    template <typename DerivedP, typename DerivedF>
+    IGL_INLINE bool ear_clipping(
+      const Eigen::MatrixBase<DerivedP>& P,
+      Eigen::PlainObjectBase<DerivedF>& eF);
 
   }
 }

include/igl/predicates/polygons_to_triangles.cpp

@@ -84,10 +84,9 @@ IGL_INLINE void igl::predicates::polygons_to_triangles(
         // This is a really low quality triangulator and will contain nearly
         // degenerate elements which become degenerate or worse when unprojected
         // back to 3D.
-        igl::predicates::ear_clipping(S,RT,_I,pF,_nS);
         // igl::predicates::ear_clipping does not gracefully fail when the input
         // is not simple. Instead it (tends?) to output too few triangles.
-        if(pF.rows() < np-2)
+        if(! igl::predicates::ear_clipping(S,pF) )
         {
           // Fallback, use a fan
           //std::cout<<igl::matlab_format(S,"S")<<std::endl;

tests/include/igl/predicates/ear_clipping.cpp

@@ -1,16 +1,64 @@
 #include <test_common.h>
 #include <igl/predicates/ear_clipping.h>
 
-TEST_CASE("ear_clipping: boolean", "[igl/predicates]")
+TEST_CASE("ear_clipping: simple poly 1", "[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::VectorXi RT,nR,I;
   Eigen::MatrixXi eF;
-  Eigen::MatrixXd nP;
   RT.setZero(polygon.rows());
-  igl::predicates::ear_clipping(polygon,RT,M,eF,nP);
-  REQUIRE(nP.rows() == 0);
 
+  igl::predicates::ear_clipping(polygon,RT,eF,I);
+  REQUIRE(I.size() == 0);
+  Eigen::MatrixXi eF1;
+  REQUIRE(igl::predicates::ear_clipping(polygon,eF1));
+  REQUIRE(eF1.rows() == polygon.rows()-2);
+
+  // Check that reverse also works
+  Eigen::MatrixXd polygon2 = polygon.colwise().reverse();
+  Eigen::MatrixXi eF2;
+  REQUIRE(igl::predicates::ear_clipping(polygon2,eF2));
+  REQUIRE(eF2.rows() == polygon2.rows()-2);
+}
+
+TEST_CASE("ear_clipping: simple poly 2", "[igl/predicates]")
+{
+  using namespace std;
+  const Eigen::MatrixXd P = (Eigen::MatrixXd(6,2)<<
+                        0,-0.212132034355964,
+        0.212132034355964,                 0,
+        0.106066017177982, 0.106066017177982,
+                        0,                 0,
+       -0.106066017177982, 0.106066017177982,
+       -0.212132034355964,                 0
+      ).finished();
+  const Eigen::VectorXi RT = Eigen::VectorXi::Zero(P.rows(),1);
+  Eigen::VectorXi I;
+  Eigen::MatrixXi eF;
+  igl::predicates::ear_clipping(P,RT,eF,I);
+  Eigen::MatrixXi ans(4, 3);
+  ans << 0,1,2,0,2,3,5,0,3,5,3,4;
+  REQUIRE(ans == eF);
+  REQUIRE(I.size() ==0);
+
+}
+
+TEST_CASE("ear_clipping: simple poly 3", "[igl/predicates]")
+{
+  using namespace std;
+  const Eigen::MatrixXd P = (Eigen::MatrixXd(3,2)<<
+      0, 0, 
+      1, 0,
+      0, 1
+      ).finished();
+  const Eigen::VectorXi RT = Eigen::VectorXi::Zero(P.rows(),1);
+  Eigen::VectorXi I;
+  Eigen::MatrixXi eF;
+  igl::predicates::ear_clipping(P,RT,eF,I);
+  Eigen::MatrixXi ans(1, 3);
+  ans << 2, 0, 1;
+  REQUIRE(ans == eF);
+  REQUIRE(I.size() == 0);
 }