Issue 521 (#1327)

* Added unit test for issue #521

Unit test adapted from https://github.com/qnzhou/libigl-unit-tests/blob/master/include/igl/delaunay_triangulation.cpp

* Fix iterator out of range error in lexicographic_triangulation (closes #521)

The iterator is now always decreased, so it points to the current left
side of the triangulation again after erasing the current item.
When the iterator is at the beginning of the boundary, it needs to point
to the last element (boundary.end()-1) after decreasing, so we have to
decrease the iterator in all cases.

In addition, the line
    if (itr == left_itr || itr == right_itr) continue;
is simplified to
    if (itr == left_itr) continue;

because the condition itr == right_itr will never be hit, because the
line before breaks the loop, when itr == right_itr:
    if (itr == right_itr) break;

include/igl/lexicographic_triangulation.cpp

@@ -95,13 +95,12 @@ IGL_INLINE void igl::lexicographic_triangulation(
       for (auto itr=left_itr; itr!=right_itr; itr++) {
         if (itr == boundary.end()) itr = boundary.begin();
         if (itr == right_itr) break;
-        if (itr == left_itr || itr == right_itr) continue;
+        if (itr == left_itr) continue;
         itr = boundary.erase(itr);
         if (itr == boundary.begin()) {
             itr = boundary.end();
-        } else {
-            itr--;
         }
+        itr--;
       }
 
       if (right_itr == boundary.begin()) {

tests/include/igl/delaunay_triangulation.cpp

@@ -0,0 +1,103 @@
+#ifdef IGL_STATIC_LIBRARY
+#undef IGL_STATIC_LIBRARY
+#endif
+
+#include <test_common.h>
+#include <igl/delaunay_triangulation.h>
+
+namespace git_issue {
+
+constexpr static double EPSILON_LENGTH = 0.0005; // Sketchup support 1/1000 precision
+constexpr static double EPSILON_ANGLE = 0.0000000000005;
+constexpr static double INV_EPSILON_LENGTH = 2000.0;
+
+template<typename T>
+inline int orient2d(const T &pa, const T &pb, const T &pc) {
+    double acx, bcx, acy, bcy;
+    acx = pa[0] - pc[0];
+    bcx = pb[0] - pc[0];
+    acy = pa[1] - pc[1];
+    bcy = pb[1] - pc[1];
+
+    double val = acx * bcy - acy * bcx;
+    if (val < -EPSILON_LENGTH) {
+        return -1;
+    } else if (val > EPSILON_LENGTH) {
+        return 1;
+    } else {
+        return 0;
+    }
+    return 0;
+}
+
+template <typename T>
+inline int incircle(const T &pa, const T &pb, const T &pc, const T &pd) {
+    double adx, ady, bdx, bdy, cdx, cdy;
+    double abdet, bcdet, cadet;
+    double alift, blift, clift;
+
+    adx = pa[0] - pd[0];
+    ady = pa[1] - pd[1];
+    bdx = pb[0] - pd[0];
+    bdy = pb[1] - pd[1];
+    cdx = pc[0] - pd[0];
+    cdy = pc[1] - pd[1];
+
+    abdet = adx * bdy - bdx * ady;
+    bcdet = bdx * cdy - cdx * bdy;
+    cadet = cdx * ady - adx * cdy;
+    alift = adx * adx + ady * ady;
+    blift = bdx * bdx + bdy * bdy;
+    clift = cdx * cdx + cdy * cdy;
+
+    double val = alift * bcdet + blift * cadet + clift * abdet;
+    if (val < -EPSILON_LENGTH) {
+        return -1;
+    } else if (val > EPSILON_LENGTH) {
+        return 1;
+    } else {
+        return 0;
+    }
+    return 0;
+}
+}
+
+
+TEST_CASE("delaunay_triangulation_issue_521", "[igl]") {
+    using namespace Eigen;
+    using namespace git_issue;
+    MatrixXd V(16, 2);
+    MatrixXi F;
+
+    V << 4.55E-13, 2.33E-12,
+        248.718, 249.939,
+        463.602, 36.1059,
+        764.953, 338.937,
+        -1002.42, 2097.68,
+        -1303.78, 1794.85,
+        -1120.79, 1612.75,
+        -1369.5, 1362.81,
+        -1552.49, 1544.91,
+        -1843.04, 1252.94,
+        -75.6625, -505.806,
+        214.883, -213.834,
+        -191.721, 190.784,
+        -1166.14, 1160.44,
+        -917.424, 1410.38,
+        56.9975, 440.723;
+
+    const auto &orient2d_predicates = [](const double * pa, const double * pb, const double * pc) {
+        return orient2d(pa, pb, pc);
+    };
+    const auto &incircle_predicates = [](const double * pa, const double * pb, const double * pc, const double * pd) {
+        return incircle(pa, pb, pc, pd);
+    };
+
+    igl::delaunay_triangulation(V, orient2d_predicates, incircle_predicates, F);
+
+    REQUIRE(F.rows() > 0);
+    REQUIRE(F.cols() == 3);
+    REQUIRE(F.maxCoeff() < 16);
+    REQUIRE(F.minCoeff() == 0);
+}
+