Fix crash in delaunay triangulation due to epsilon issues.

2 files changed, 647 insertions, 211 deletions

diff --git a/source/blender/blenlib/intern/delaunay_2d.c b/source/blender/blenlib/intern/delaunay_2d.c
index 27f7a2a4430..350f99ce434 100644
--- a/source/blender/blenlib/intern/delaunay_2d.c
+++ b/source/blender/blenlib/intern/delaunay_2d.c
@@ -53,6 +53,7 @@ typedef struct CDTVert {
   SymEdge *symedge;    /* Some edge attached to it. */
   LinkNode *input_ids; /* List of corresponding vertex input ids. */
   int index;           /* Index into array that cdt keeps. */
+  int visit_index;     /* Which visit epoch has this been seen. */
 } CDTVert;
 
 typedef struct CDTEdge {
@@ -110,7 +111,9 @@ static void dump_v(const CDTVert *v, const char *lab);
 static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit);
 static void dump_id_list(const LinkNode *id_list, const char *lab);
 static void dump_cdt(const CDT_state *cdt, const char *lab);
+static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab);
 static void cdt_draw(CDT_state *cdt, const char *lab);
+static void write_cdt_input_to_file(const CDT_input *inp);
 static void validate_face_centroid(SymEdge *se);
 static void validate_cdt(CDT_state *cdt, bool check_all_tris);
 #endif
@@ -151,14 +154,11 @@ static int CCW_test(const double a[2], const double b[2], const double c[2], con
   return 0;
 }
 
-/** return true if a -- b -- c are in that order, assuming they are on a straight line. */
-static bool in_line(const double a[2], const double b[2], const double c[2])
+/** return true if a -- b -- c are in that order, assuming they are on a straight line according to
+ * CCW_test. */
+static bool in_line(const double a[2], const double b[2], const double c[2], double eps)
 {
-  double dir_ab[2], dir_ac[2];
-
-  sub_v2_v2v2_db(dir_ab, a, b);
-  sub_v2_v2v2_db(dir_ac, a, c);
-  return dot_v2v2_db(dir_ab, dir_ac) >= 0.0;
+  return fabs(len_v2v2_db(a, c) - (len_v2v2_db(a, b) + len_v2v2_db(b, c))) <= eps;
 }
 
 #ifndef NDEBUG
@@ -208,6 +208,7 @@ static CDTVert *add_cdtvert(CDT_state *cdt, double x, double y)
   }
   BLI_assert(cdt->vert_array_len < cdt->vert_array_len_alloc);
   v->index = cdt->vert_array_len;
+  v->visit_index = 0;
   cdt->vert_array[cdt->vert_array_len++] = v;
   return v;
 }
@@ -365,8 +366,8 @@ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2)
   CDTFace *fold, *fnew;
   SymEdge *sdiag, *sdiagsym;
   SymEdge *s1prev, *s1prevsym, *s2prev, *s2prevsym, *se;
-  BLI_assert(reachable(s1, s2, 20));
-  BLI_assert(reachable(s2, s1, 20));
+  BLI_assert(reachable(s1, s2, 2000));
+  BLI_assert(reachable(s2, s1, 2000));
   fold = s1->face;
   fnew = add_cdtface(cdt);
   s1prev = prev(s1);
@@ -385,7 +386,7 @@ static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2)
   s2->rot = sdiagsym;
   sdiagsym->rot = s2prevsym;
 #ifdef DEBUG_CDT
-  BLI_assert(reachable(s2, sdiag, 20));
+  BLI_assert(reachable(s2, sdiag, 2000));
 #endif
   for (se = s2; se != sdiag; se = se->next) {
     se->face = fnew;
@@ -1505,7 +1506,8 @@ static void add_edge_constraint(
   SymEdge *se;
   CDTEdge *edge;
   int ccw1, ccw2, isect;
-  int i, search_count;
+  int i, search_count, visit;
+  double curco[2];
   double lambda;
   const double epsilon = cdt->epsilon;
   bool done, state_through_vert;
@@ -1553,7 +1555,7 @@ static void add_edge_constraint(
    *    in = a
    *    out = b
    *
-   * crossings[0] will have in = NULL, and crossings[last] will have out = NULL
+   * crossings[0] will have in = NULL, and crossings[last] will have out = NULL.
    */
   if (r_edges) {
     *r_edges = NULL;
@@ -1569,6 +1571,9 @@ static void add_edge_constraint(
       dump_se(vse1, "  se1");
       dump_se(vse2, "  se2");
     }
+    if (dbg_level > 2) {
+      dump_cdt(cdt, "before insert_segment");
+    }
   }
 #endif
   if (v1 == v2) {
@@ -1579,181 +1584,292 @@ static void add_edge_constraint(
 #endif
     return;
   }
-  state_through_vert = true;
+  /* Quick check for special case where segment is one of edges from v1. */
+  tstart = t = vse1;
   done = false;
-  t = vse1;
-  search_count = 0;
-  while (!done) {
-    /* Invariant: crossings[0 .. BLI_array_len(crossings)] has crossing info for path up to
-     * but not including the crossing of edge t, which will either be through a vert
-     * (if state_through_vert is true) or through edge t not at either end.
-     * In the latter case, t->face is the face that ray v1--v2 goes through after path-so-far.
-     */
-#ifdef DEBUG_CDT
-    if (dbg_level > 1) {
-      fprintf(
-          stderr, "top of insert_segment main loop, state_through_vert=%d\n", state_through_vert);
-      dump_se_cycle(t, "current t ", 4);
-    }
-#endif
-    if (state_through_vert) {
-      /* Invariant: ray v1--v2 contains t->vert. */
-      cdata.in = (BLI_array_len(crossings) == 0) ? NULL : t;
-      cdata.out = NULL; /* To be filled in if this isn't final. */
+  do {
+    if (t->next->vert == v2) {
+      cdata.in = NULL;
+      cdata.out = t->next;
       cdata.lambda = 0.0;
       cdata.vert = t->vert;
       BLI_array_append(crossings, cdata);
-      if (t->vert == v2) {
+      cdata.in = t->next->rot;
+      cdata.out = NULL;
+      cdata.lambda = 0.0;
+      cdata.vert = v2;
+      BLI_array_append(crossings, cdata);
 #ifdef DEBUG_CDT
-        if (dbg_level > 0) {
-          fprintf(stderr, "found v2, so done\n");
-        }
+      if (dbg_level > 1) {
+        fprintf(stderr, "special one segment case\n");
+        dump_se(t, "  ");
+      }
 #endif
-        done = true;
+      done = true;
+      break;
+    }
+    t = t->rot;
+  } while (t != tstart);
+  if (!done) {
+    /* To prevent infinite loop in the face of epsilon tests that might lead us back to
+     * an already-visited (vertex, face) pair, use visit indices.
+     */
+    visit = ++cdt->visit_count;
+    state_through_vert = true;
+    done = false;
+    t = vse1;
+    search_count = 0;
+    while (!done) {
+      /* Invariant: crossings[0 .. BLI_array_len(crossings)] has crossing info for path up to
+       * but not including the crossing of edge t, which will either be through a vert
+       * (if state_through_vert is true) or through edge t not at either end.
+       * In the latter case, t->face is the face that ray v1--v2 goes through after path-so-far.
+       * Rather than continually looking for intersection of v1--v2, we keep track of
+       * last vertex or intersection point in curco,because that may be slightly off the ray
+       * v1--v2.
+       */
+#ifdef DEBUG_CDT
+      if (dbg_level > 1) {
+        fprintf(stderr,
+                "top of insert_segment main loop, state_through_vert=%d\n",
+                state_through_vert);
+        dump_se_cycle(t, "current t ", 4);
       }
-      else {
-        /* Do ccw scan of triangles around t->vert to find exit triangle for ray v1--v2. */
-        tstart = t;
-        tout = NULL;
-        do {
-          va = t->next->vert;
-          vb = t->next->next->vert;
-          ccw1 = CCW_test(t->vert->co, va->co, v2->co, epsilon);
-          ccw2 = CCW_test(t->vert->co, vb->co, v2->co, epsilon);
+#endif
+      BLI_assert(t->vert->visit_index != visit || t->face->visit_index != visit);
+      t->vert->visit_index = visit;
+      t->face->visit_index = visit;
+      if (state_through_vert) {
+        /* Invariant: ray vcur--v2 contains t->vert. */
+        cdata.in = (BLI_array_len(crossings) == 0) ? NULL : t;
+        cdata.out = NULL; /* To be filled in if this isn't final. */
+        cdata.lambda = 0.0;
+        cdata.vert = t->vert;
+        BLI_array_append(crossings, cdata);
+        if (t->vert == v2) {
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "non-final through vert case\n");
-            dump_v(va, " va");
-            dump_v(vb, " vb");
-            fprintf(stderr, "ccw1=%d, ccw2=%d\n", ccw1, ccw2);
+          if (dbg_level > 0) {
+            fprintf(stderr, "found v2, so done\n");
           }
 #endif
-          if (ccw1 == 0 && in_line(t->vert->co, va->co, v2->co)) {
+          done = true;
+        }
+        else {
+          /* Do ccw scan of triangles around t->vert to find exit triangle for ray vcur--v2. */
+          tstart = t;
+          tout = NULL;
+          do {
+            va = t->next->vert;
+            vb = t->next->next->vert;
+            ccw1 = CCW_test(t->vert->co, va->co, v2->co, epsilon);
+            ccw2 = CCW_test(t->vert->co, vb->co, v2->co, epsilon);
+#ifdef DEBUG_CDT
+            if (dbg_level > 1) {
+              fprintf(stderr, "non-final through vert case\n");
+              dump_v(va, " va");
+              dump_v(vb, " vb");
+              fprintf(stderr, "ccw1=%d, ccw2=%d\n", ccw1, ccw2);
+            }
+#endif
+            if (ccw1 == 0 && in_line(t->vert->co, va->co, v2->co, epsilon) &&
+                va->visit_index != visit) {
 #ifdef DEBUG_CDT
-            if (dbg_level > 0) {
-              fprintf(stderr, "ray goes through va\n");
+              if (dbg_level > 0) {
+                fprintf(stderr, "ray goes through va\n");
+              }
+#endif
+              state_through_vert = true;
+              tout = t;
+              t = t->next;
+              break;
             }
+            else if (ccw2 == 0 && in_line(t->vert->co, vb->co, v2->co, epsilon) &&
+                     vb->visit_index != visit) {
+#ifdef DEBUG_CDT
+              if (dbg_level > 0) {
+                fprintf(stderr, "ray goes through vb\n");
+              }
 #endif
-            state_through_vert = true;
-            tout = t;
-            t = t->next;
-            break;
-          }
-          else if (ccw2 == 0 && in_line(t->vert->co, vb->co, v2->co)) {
+              state_through_vert = true;
+              t = t->next->next;
+              tout = sym(t);
+              break;
+            }
+            else if (ccw1 > 0 && ccw2 < 0 &&
+                     (t->next->vert->visit_index != visit ||
+                      t->next->face->visit_index != visit)) {
+#ifdef DEBUG_CDT
+              if (dbg_level > 0) {
+                fprintf(stderr, "segment intersects\n");
+              }
+#endif
+              state_through_vert = false;
+              tout = t;
+              t = t->next;
+              break;
+            }
+            t = t->rot;
 #ifdef DEBUG_CDT
-            if (dbg_level > 0) {
-              fprintf(stderr, "ray goes through vb\n");
+            if (dbg_level > 1) {
+              dump_se_cycle(t, "next rot tri", 4);
             }
 #endif
-            state_through_vert = true;
-            t = t->next->next;
-            tout = sym(t);
-            break;
-          }
-          else if (ccw1 > 0 && ccw2 < 0) {
+          } while (t != tstart);
+          if (tout == NULL) {
+            /* With exact arithmetic this shouldn't happen, but maybe the epsilon tests made it so
+             * that we want to go back to a previous vertex.
+             * As desperation measure, pick unvisited vertex that is closest in line with
+             * destination.
+             */
 #ifdef DEBUG_CDT
             if (dbg_level > 0) {
-              fprintf(stderr, "segment intersects\n");
+              fprintf(stderr, "add_edge_constraint desperation search\n");
             }
 #endif
-            state_through_vert = false;
-            tout = t;
-            t = t->next;
-            break;
-          }
-          t = t->rot;
+            SymEdge *bestt = NULL;
+            double dot, bestdot = -2.0;
+            double dir_tv_v2[2], dir_tvnext_v2[2];
+            sub_v2_v2v2_db(dir_tv_v2, v2->co, t->vert->co);
+            do {
+              if (t->next->vert->visit_index != visit) {
+                sub_v2_v2v2_db(dir_tvnext_v2, v2->co, t->next->vert->co);
+                dot = dot_v2v2_db(dir_tv_v2, dir_tvnext_v2);
+                if (dot > bestdot) {
+                  bestdot = dot;
+                  bestt = t->next;
+                }
+              }
+              t = t->rot;
+            } while (t != tstart);
+            if (bestt == NULL) {
+              /* No unvisited place to go! Give up on adding this edge constraint. */
+#ifdef DEBUG_CDT
+              fprintf(stderr, "could not add edge constraint\n");
+#endif
+              return;
+            }
+            else {
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            dump_se_cycle(t, "next rot tri", 4);
-          }
+              if (dbg_level > 0) {
+                fprintf(stderr, "add_edge_constraint desperation search chose to go through\n");
+                dump_v(bestt->vert, "desperation vert");
+              }
 #endif
-        } while (t != tstart);
-        BLI_assert(tout != NULL); /* TODO: something sensible for "this can't happen" */
-        crossings[BLI_array_len(crossings) - 1].out = tout;
-      }
-    }
-    else { /* State is "through edge", not "through vert" */
-      /* Invariant: ray v1--v2 intersects segment t->edge, not at either end.
-       * and t->face is the face we have just passed through. */
-      va = t->vert;
-      vb = t->next->vert;
-#ifdef DEBUG_CDT
-      if (dbg_level > 1) {
-        fprintf(stderr, "through edge case\n");
-        dump_v(va, " va");
-        dump_v(vb, " vb");
+              tout = bestt;
+              t = t->next;
+            }
+          }
+          crossings[BLI_array_len(crossings) - 1].out = tout;
+        }
       }
-#endif
-      isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, v1->co, v2->co, &lambda, NULL);
-      /* TODO: something sensible for "this can't happen" */
-      BLI_assert(isect == ISECT_LINE_LINE_CROSS);
-      UNUSED_VARS_NDEBUG(isect);
+      else { /* State is "through edge", not "through vert" */
+        /* Invariant: ray v1--v2 intersects segment t->edge, not at either end.
+         * and t->face is the face we have just passed through.
+         * Whatever we make t next should not have both vert and face visited.
+         */
+        va = t->vert;
+        vb = t->next->vert;
+        /* Get curco; cdata should have data from last time through the loop still. */
+        if (cdata.lambda == 0.0) {
+          copy_v3_v3_db(curco, cdata.vert->co);
+        }
+        else {
+          interp_v2_v2v2_db(curco, cdata.in->vert->co, cdata.in->next->vert->co, lambda);
+        }
 #ifdef DEBUG_CDT
-      if (dbg_level > 0) {
-        fprintf(stderr, "intersect point at %f along va--vb\n", lambda);
-        if (dbg_level == 1) {
+        if (dbg_level > 1) {
+          fprintf(stderr, "through edge case, curco=(%f,%f)\n", F2(curco));
           dump_v(va, " va");
           dump_v(vb, " vb");
         }
-      }
 #endif
-      tout = sym(t)->next;
-      cdata.in = t;
-      cdata.out = tout;
-      cdata.lambda = lambda;
-      cdata.vert = NULL; /* To be filled in with edge split vertex later. */
-      BLI_array_append(crossings, cdata);
+        isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, curco, v2->co, &lambda, NULL);
+        if (isect != ISECT_LINE_LINE_CROSS) {
+          /* Shouldn't happen. Just pick something. */
 #ifdef DEBUG_CDT
-      if (dbg_level > 0) {
-        dump_se_cycle(tout, "next search tri", 4);
-      }
+          if (dbg_level > 1) {
+            fprintf(stderr, "add_edge_constraint no intersect found, using lambda = 0.5\n");
+          }
 #endif
-      /* 'tout' is 'symedge' from 'vb' to third vertex, 'vc'. */
-      BLI_assert(tout->vert == va);
-      vc = tout->next->vert;
-      ccw1 = CCW_test(v1->co, v2->co, vc->co, epsilon);
+          lambda = 0.5;
+        }
 #ifdef DEBUG_CDT
-      if (dbg_level > 1) {
-        fprintf(stderr, "now searching with third vertex ");
-        dump_v(vc, "vc");
-        fprintf(stderr, "ccw(v1, v2, vc) = %d\n", ccw1);
-      }
+        if (dbg_level > 0) {
+          fprintf(stderr, "intersect point at %f along va--vb\n", lambda);
+          if (dbg_level == 1) {
+            dump_v(va, " va");
+            dump_v(vb, " vb");
+          }
+        }
 #endif
-      if (ccw1 == -1) {
-        /* v1--v2 should intersect vb--vc. */
+        tout = sym(t)->next;
+        cdata.in = t;
+        cdata.out = tout;
+        cdata.lambda = lambda;
+        cdata.vert = NULL; /* To be filled in with edge split vertex later. */
+        BLI_array_append(crossings, cdata);
 #ifdef DEBUG_CDT
-        if (dbg_level > 1) {
-          fprintf(stderr, "v1--v2 intersects vb--vc\n");
+        if (dbg_level > 0) {
+          dump_se_cycle(tout, "next search tri", 4);
+          dump_se(tout, "tout");
         }
 #endif
-        t = tout->next;
-        state_through_vert = false;
-      }
-      else if (ccw1 == 1) {
-        /* v1--v2 should intersect va--vc. */
+        /* 'tout' is 'symedge' from 'va' to third vertex, 'vc'. */
+        BLI_assert(tout->vert == va);
+        vc = tout->next->vert;
+        ccw1 = CCW_test(curco, v2->co, vc->co, epsilon);
 #ifdef DEBUG_CDT
         if (dbg_level > 1) {
-          fprintf(stderr, "v1--v2 intersects va--vc\n");
+          fprintf(stderr, "now searching with third vertex ");
+          dump_v(vc, "vc");
+          fprintf(stderr, "ccw(vcur, v2, vc) = %d\n", ccw1);
         }
 #endif
-        t = tout;
-        state_through_vert = false;
-      }
-      else {
-        /* ccw1 == 0. */
+        if (ccw1 == -1 && (vc->visit_index != visit || tout->next->face->visit_index != visit)) {
+          /* vcur--v2 should intersect vb--vc. */
 #ifdef DEBUG_CDT
-        if (dbg_level > 1) {
-          fprintf(stderr, "ccw==0 case, so going through or to vc\n");
+          if (dbg_level > 1) {
+            fprintf(stderr, "curco--v2 intersects vb--vc\n");
+          }
+#endif
+          t = tout->next;
+          state_through_vert = false;
         }
+        else if (ccw1 == 1 && tout->face->visit_index != visit) {
+          /* vcur--v2 should intersect va--vc. */
+#ifdef DEBUG_CDT
+          if (dbg_level > 1) {
+            fprintf(stderr, "curco--v2 intersects va--vc\n");
+          }
+#endif
+          t = tout;
+          state_through_vert = false;
+        }
+        else if (ccw1 == 0 && vc->visit_index != visit) {
+#ifdef DEBUG_CDT
+          if (dbg_level > 1) {
+            fprintf(stderr, "ccw==0 case, so going through or to vc\n");
+          }
 #endif
-        t = tout->next;
-        state_through_vert = true;
+          t = tout->next;
+          state_through_vert = true;
+        }
+        else {
+#ifdef DEBUG_CDT
+          fprintf(stderr, "add_edge_constraint desperation search 2\n");
+#endif
+          if (tout->face->visit_index != visit) {
+            /* Treat as if an intersection of va-vc. */
+            t = tout;
+            state_through_vert = false;
+          }
+        }
+      }
+      if (++search_count > 10000) {
+        fprintf(stderr, "infinite loop? bailing out\n");
+        BLI_assert(0); /* Catch these while developing. */
+        break;
       }
-    }
-    if (++search_count > 10000) {
-      fprintf(stderr, "infinite loop? bailing out\n");
-      BLI_assert(0); /* Catch these while developing. */
-      break;
     }
   }
 #ifdef DEBUG_CDT
@@ -1968,7 +2084,7 @@ static void dissolve_symedge(CDT_state *cdt, SymEdge *se)
     if (se->face->symedge == se) {
       se->face->symedge = se->next;
     }
-    if (symse->face->symedge == se) {
+    if (symse->face->symedge == symse) {
       symse->face->symedge = symse->next;
     }
   }
@@ -2380,7 +2496,13 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
   CDTEdge *face_edge;
   SymEdge *face_symedge;
 #ifdef DEBUG_CDT
-  int dbg_level = 1;
+  int dbg_level = 0;
+#endif
+
+#ifdef DEBUG_CDT
+  if (dbg_level == -1) {
+    write_cdt_input_to_file(input);
+  }
 #endif
 
   if ((nv > 0 && input->vert_coords == NULL) || (ne > 0 && input->edges == NULL) ||
@@ -2446,6 +2568,15 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
       continue;
     }
     add_edge_constraint(cdt, verts[v1], verts[v2], i, NULL);
+#ifdef DEBUG_CDT
+    if (dbg_level > 3) {
+      char namebuf[60];
+      sprintf(namebuf, "after edge constraint %d = (%d,%d)\n", i, v1, v2);
+      cdt_draw(cdt, namebuf);
+      dump_cdt(cdt, namebuf);
+      validate_cdt(cdt, true);
+    }
+#endif
   }
 #ifdef DEBUG_CDT
   if (dbg_level > 2) {
@@ -2579,20 +2710,37 @@ void BLI_delaunay_2d_cdt_free(CDT_result *result)
 
 #ifdef DEBUG_CDT
 
-static void dump_se(const SymEdge *se, const char *lab)
+static const char *vertname(const CDTVert *v)
 {
-  if (se->next) {
-    fprintf(
-        stderr, "%s((%.3f,%.3f)->(%.3f,%.3f))\n", lab, F2(se->vert->co), F2(se->next->vert->co));
+  static char vertnamebuf[20];
+
+  if (v->index < 4) {
+    sprintf(vertnamebuf, "[%c]", "ABCD"[v->index]);
   }
   else {
-    fprintf(stderr, "%s((%.3f,%.3f)->NULL)\n", lab, F2(se->vert->co));
+    sprintf(vertnamebuf, "[%d]", v->index - 4);
   }
+  return vertnamebuf;
 }
 
 static void dump_v(const CDTVert *v, const char *lab)
 {
-  fprintf(stderr, "%s(%.3f,%.3f)\n", lab, F2(v->co));
+  fprintf(stderr, "%s%s(%.3f,%.3f)\n", lab, vertname(v), F2(v->co));
+}
+
+static void dump_se(const SymEdge *se, const char *lab)
+{
+  if (se->next) {
+    fprintf(stderr,
+            "%s%s((%.3f,%.3f)->(%.3f,%.3f))\n",
+            lab,
+            vertname(se->vert),
+            F2(se->vert->co),
+            F2(se->next->vert->co));
+  }
+  else {
+    fprintf(stderr, "%s%s((%.3f,%.3f)->NULL)\n", lab, vertname(se->vert), F2(se->vert->co));
+  }
 }
 
 static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit)
@@ -2622,21 +2770,12 @@ static void dump_id_list(const LinkNode *id_list, const char *lab)
   }
 }
 
-static const char *vertname(CDTVert *v)
-{
-  static char vertnamebuf[20];
-
-  if (v->index < 4) {
-    sprintf(vertnamebuf, "[%c]", "ABCD"[v->index]);
-  }
-  else {
-    sprintf(vertnamebuf, "[%d]", v->index - 4);
-  }
-  return vertnamebuf;
-}
-
+/* If filter_fn != NULL, only dump vert v its edges when filter_fn(cdt, v, filter_data) is true. */
 #  define PL(p) (POINTER_AS_UINT(p) & 0xFFFF)
-static void dump_cdt(const CDT_state *cdt, const char *lab)
+static void dump_cdt_filtered(const CDT_state *cdt,
+                              bool (*filter_fn)(const CDT_state *, int, void *),
+                              void *filter_data,
+                              const char *lab)
 {
   LinkNode *ln;
   CDTVert *v, *vother;
@@ -2648,6 +2787,9 @@ static void dump_cdt(const CDT_state *cdt, const char *lab)
   fprintf(stderr, "\nCDT %s\n", lab);
   fprintf(stderr, "\nVERTS\n");
   for (i = 0; i < cdt->vert_array_len; i++) {
+    if (filter_fn && !filter_fn(cdt, i, filter_data)) {
+      continue;
+    }
     v = cdt->vert_array[i];
     fprintf(stderr, "%s %x: (%f,%f) symedge=%x\n", vertname(v), PL(v), F2(v->co), PL(v->symedge));
     dump_id_list(v->input_ids, "  ");
@@ -2664,6 +2806,9 @@ static void dump_cdt(const CDT_state *cdt, const char *lab)
       fprintf(stderr, "\n");
     }
   }
+  if (filter_fn) {
+    return;
+  }
   fprintf(stderr, "\nEDGES\n");
   for (ln = cdt->edges; ln; ln = ln->next) {
     e = (CDTEdge *)ln->link;
@@ -2709,6 +2854,55 @@ static void dump_cdt(const CDT_state *cdt, const char *lab)
 }
 #  undef PL
 
+static void dump_cdt(const CDT_state *cdt, const char *lab)
+{
+  dump_cdt_filtered(cdt, NULL, NULL, lab);
+}
+
+typedef struct ReachableFilterData {
+  int vstart_index;
+  int maxdist;
+} ReachableFilterData;
+
+/* Stupid algorithm will repeat itself. Don't use for large n. */
+static bool reachable_filter(const CDT_state *cdt, int v_index, void *filter_data)
+{
+  CDTVert *v;
+  SymEdge *se;
+  ReachableFilterData *rfd_in = (ReachableFilterData *)filter_data;
+  ReachableFilterData rfd_next;
+
+  if (v_index == rfd_in->vstart_index) {
+    return true;
+  }
+  if (rfd_in->maxdist <= 0 || v_index < 0 || v_index >= cdt->vert_array_len) {
+    return false;
+  }
+  else {
+    v = cdt->vert_array[v_index];
+    se = v->symedge;
+    if (se != NULL) {
+      rfd_next.vstart_index = rfd_in->vstart_index;
+      rfd_next.maxdist = rfd_in->maxdist - 1;
+      do {
+        if (reachable_filter(cdt, se->next->vert->index, &rfd_next)) {
+          return true;
+        }
+        se = se->rot;
+      } while (se != v->symedge);
+    }
+  }
+  return false;
+}
+
+static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab)
+{
+  ReachableFilterData rfd;
+  rfd.vstart_index = v;
+  rfd.maxdist = maxdist;
+  dump_cdt_filtered(cdt, reachable_filter, &rfd, lab);
+}
+
 /**
  * Make an html file with svg in it to display the argument cdt.
  * Mouse-overs will reveal the coordinates of vertices and edges.
@@ -2716,8 +2910,12 @@ static void dump_cdt(const CDT_state *cdt, const char *lab)
  * The first call creates DRAWFILE; subsequent calls append to it.
  */
 #  define DRAWFILE "/tmp/debug_draw.html"
-#  define MAX_DRAW_WIDTH 1000
-#  define MAX_DRAW_HEIGHT 700
+#  define MAX_DRAW_WIDTH 2000
+#  define MAX_DRAW_HEIGHT 1400
+#  define THIN_LINE 1
+#  define THICK_LINE 3
+#  define VERT_RADIUS 3
+#  define DRAW_VERT_LABELS 1
 static void cdt_draw(CDT_state *cdt, const char *lab)
 {
   static bool append = false;
@@ -2727,9 +2925,7 @@ static void cdt_draw(CDT_state *cdt, const char *lab)
   double maxx = cdt->maxx + draw_margin;
   double miny = cdt->miny - draw_margin;
   double maxy = cdt->maxy + draw_margin;
-  double width = maxx - minx;
-  double height = maxy - miny;
-  double aspect = height / width;
+  double width, height, aspect;
   int view_width, view_height;
   double scale;
   LinkNode *ln;
@@ -2737,6 +2933,10 @@ static void cdt_draw(CDT_state *cdt, const char *lab)
   CDTEdge *e;
   int i, strokew;
 
+  /* Note: to debug a small area: assign custom min's/max's here. */
+  width = maxx - minx;
+  height = maxy - miny;
+  aspect = height / width;
   view_width = MAX_DRAW_WIDTH;
   view_height = (int)(view_width * aspect);
   if (view_height > MAX_DRAW_HEIGHT) {
@@ -2767,28 +2967,36 @@ static void cdt_draw(CDT_state *cdt, const char *lab)
     }
     u = e->symedges[0].vert;
     v = e->symedges[1].vert;
-    strokew = is_constrained_edge(e) ? 5 : 2;
+    strokew = is_constrained_edge(e) ? THICK_LINE : THIN_LINE;
     fprintf(f,
             "<line fill=\"none\" stroke=\"black\" stroke-width=\"%d\" "
-            "x1=\"%.3f\" y1=\"%.3f\" x2=\"%.3f\" y2=\"%.3f\">\n",
+            "x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\">\n",
             strokew,
             SX(u->co[0]),
             SY(u->co[1]),
             SX(v->co[0]),
             SY(v->co[1]));
-    fprintf(
-        f, "  <title>(%.3f,%.3f)(%.3f,%.3f)</title>\n", u->co[0], u->co[1], v->co[0], v->co[1]);
+    fprintf(f, "  <title>%s", vertname(u));
+    fprintf(f, "%s</title>\n", vertname(v));
     fprintf(f, "</line>\n");
   }
   i = cdt->output_prepared ? NUM_BOUND_VERTS : 0;
   for (; i < cdt->vert_array_len; i++) {
     v = cdt->vert_array[i];
     fprintf(f,
-            "<circle fill=\"black\" cx=\"%.3f\" cy=\"%.3f\" r=\"5\">\n",
+            "<circle fill=\"black\" cx=\"%f\" cy=\"%f\" r=\"%d\">\n",
             SX(v->co[0]),
-            SY(v->co[1]));
-    fprintf(f, "  <title>(%.3f,%.3f)</title>\n", v->co[0], v->co[1]);
+            SY(v->co[1]),
+            VERT_RADIUS);
+    fprintf(f, "  <title>%s(%.3f,%.3f)</title>\n", vertname(v), v->co[0], v->co[1]);
     fprintf(f, "</circle>\n");
+#  if DRAW_VERT_LABELS
+    fprintf(f,
+            "<text x=\"%f\" y=\"%f\" font-size=\"small\">%s</text>\n",
+            SX(v->co[0]) + VERT_RADIUS,
+            SY(v->co[1]) - VERT_RADIUS,
+            vertname(v));
+#  endif
   }
 
   fprintf(f, "</svg>\n</div>\n");
@@ -2798,6 +3006,28 @@ static void cdt_draw(CDT_state *cdt, const char *lab)
 #  undef SY
 }
 
+#  define CDTFILE "/tmp/cdtinput.txt"
+static void write_cdt_input_to_file(const CDT_input *inp)
+{
+  int i, j;
+  FILE *f = fopen(CDTFILE, "w");
+
+  fprintf(f, "%d %d %d\n", inp->verts_len, inp->edges_len, inp->faces_len);
+  for (i = 0; i < inp->verts_len; i++) {
+    fprintf(f, "%.17f %.17f\n", inp->vert_coords[i][0], inp->vert_coords[i][1]);
+  }
+  for (i = 0; i < inp->edges_len; i++) {
+    fprintf(f, "%d %d\n", inp->edges[i][0], inp->edges[i][1]);
+  }
+  for (i = 0; i < inp->faces_len; i++) {
+    for (j = 0; j < inp->faces_len_table[i]; j++) {
+      fprintf(f, "%d ", inp->faces[j + inp->faces_start_table[i]]);
+    }
+    fprintf(f, "\n");
+  }
+  fclose(f);
+}
+
 #  ifndef NDEBUG /* Only used in assert. */
 /**
  * Is a visible from b: i.e., ab crosses no edge of cdt?
diff --git a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc
index 8b29128eb0d..85e8ae3f141 100644
--- a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc
+++ b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc
@@ -15,6 +15,10 @@ extern "C" {
 #include <fstream>
 #include <sstream>
 
+#define DO_REGULAR_TESTS 1
+#define DO_RANDOM_TESTS 0
+#define DO_FILE_TESTS 0
+
 static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts)
 {
   r_input->verts_len = nverts;
@@ -43,6 +47,117 @@ static void add_input_faces(
   r_input->faces_len_table = faces_len_table;
 }
 
+/* The spec should have the form:
+ * #verts #edges #faces
+ * <float> <float>   [#verts lines)
+ * <int> <int>   [#edges lines]
+ * <int> <int> ... <int>   [#faces lines]
+ */
+static void fill_input_from_string(CDT_input *r_input, const char *spec)
+{
+  std::string line;
+  std::vector<std::vector<int>> faces;
+  int i, j;
+  int nverts, nedges, nfaces;
+  float(*p)[2];
+  int(*e)[2];
+  int *farr;
+  int *flen;
+  int *fstart;
+
+  std::istringstream ss(spec);
+  getline(ss, line);
+  std::istringstream hdrss(line);
+  hdrss >> nverts >> nedges >> nfaces;
+  if (nverts == 0) {
+    return;
+  }
+  p = (float(*)[2])MEM_malloc_arrayN(nverts, 2 * sizeof(float), __func__);
+  if (nedges > 0) {
+    e = (int(*)[2])MEM_malloc_arrayN(nedges, 2 * sizeof(int), __func__);
+  }
+  if (nfaces > 0) {
+    flen = (int *)MEM_malloc_arrayN(nfaces, sizeof(int), __func__);
+    fstart = (int *)MEM_malloc_arrayN(nfaces, sizeof(int), __func__);
+  }
+  i = 0;
+  while (i < nverts && getline(ss, line)) {
+    std::istringstream iss(line);
+    iss >> p[i][0] >> p[i][1];
+    i++;
+  }
+  i = 0;
+  while (i < nedges && getline(ss, line)) {
+    std::istringstream ess(line);
+    ess >> e[i][0] >> e[i][1];
+    i++;
+  }
+  i = 0;
+  while (i < nfaces && getline(ss, line)) {
+    std::istringstream fss(line);
+    int v;
+    faces.push_back(std::vector<int>());
+    while (fss >> v) {
+      faces[i].push_back(v);
+    }
+    i++;
+  }
+  fill_input_verts(r_input, p, nverts);
+  if (nedges > 0) {
+    add_input_edges(r_input, e, nedges);
+  }
+  if (nfaces > 0) {
+    for (i = 0; i < nfaces; i++) {
+      flen[i] = (int)faces[i].size();
+      if (i == 0) {
+        fstart[i] = 0;
+      }
+      else {
+        fstart[i] = fstart[i - 1] + flen[i - 1];
+      }
+    }
+    farr = (int *)MEM_malloc_arrayN(fstart[nfaces - 1] + flen[nfaces - 1], sizeof(int), __func__);
+    for (i = 0; i < nfaces; i++) {
+      for (j = 0; j < (int)faces[i].size(); j++) {
+        farr[fstart[i] + j] = faces[i][j];
+      }
+    }
+    add_input_faces(r_input, farr, fstart, flen, nfaces);
+  }
+}
+
+static void fill_input_from_file(CDT_input *in, const char *filename)
+{
+  std::FILE *fp = std::fopen(filename, "rb");
+  if (fp) {
+    std::string contents;
+    std::fseek(fp, 0, SEEK_END);
+    contents.resize(std::ftell(fp));
+    std::rewind(fp);
+    std::fread(&contents[0], 1, contents.size(), fp);
+    std::fclose(fp);
+    fill_input_from_string(in, contents.c_str());
+  }
+  else {
+    printf("couldn't open file %s\n", filename);
+  }
+}
+
+static void free_spec_arrays(CDT_input *in)
+{
+  if (in->vert_coords) {
+    MEM_freeN(in->vert_coords);
+  }
+  if (in->edges) {
+    MEM_freeN(in->edges);
+  }
+  if (in->faces_len_table) {
+    MEM_freeN(in->faces_len_table);
+    MEM_freeN(in->faces_start_table);
+    MEM_freeN(in->faces);
+  }
+}
+
 /* which output vert index goes with given input vertex? -1 if not found */
 static int get_output_vert_index(const CDT_result *r, int in_index)
 {
@@ -187,6 +302,7 @@ static void dump_result(CDT_result *r)
   }
 }
 
+#if DO_REGULAR_TESTS
 TEST(delaunay, Empty)
 {
   CDT_input in;
@@ -692,7 +808,108 @@ TEST(delaunay, TriInTri)
   BLI_delaunay_2d_cdt_free(out);
 }
 
-#if 0
+TEST(delaunay, DiamondInSquare)
+{
+  CDT_input in;
+  CDT_result *out;
+  const char *spec = R"(8 0 2
+  0.0 0.0
+  1.0 0.0
+  1.0 1.0
+  0.0 1.0
+  0.14644660940672627 0.5
+  0.5 0.14644660940672627
+  0.8535533905932737 0.5
+  0.5 0.8535533905932737
+  0 1 2 3
+  4 5 6 7
+  )";
+  fill_input_from_string(&in, spec);
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS_VALID_BMESH);
+  EXPECT_EQ(out->verts_len, 8);
+  EXPECT_EQ(out->edges_len, 10);
+  EXPECT_EQ(out->faces_len, 3);
+  free_spec_arrays(&in);
+}
+
+TEST(delaunay, DiamondInSquareWire)
+{
+  CDT_input in;
+  CDT_result *out;
+  const char *spec = R"(8 8 0
+  0.0 0.0
+  1.0 0.0
+  1.0 1.0
+  0.0 1.0
+  0.14644660940672627 0.5
+  0.5 0.14644660940672627
+  0.8535533905932737 0.5
+  0.5 0.8535533905932737
+  0 1
+  1 2
+  2 3
+  3 0
+  4 5
+  5 6
+  6 7
+  7 4
+  )";
+  fill_input_from_string(&in, spec);
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 8);
+  EXPECT_EQ(out->edges_len, 8);
+  EXPECT_EQ(out->faces_len, 2);
+  free_spec_arrays(&in);
+}
+
+TEST(delaunay, ClosePts)
+{
+  CDT_input in;
+  CDT_result *out;
+  const char *spec = R"(7 2 1
+  0.46876350045204163 0.06087132915854454
+  0.46865847706794739 0.03632887825369835
+  0.49176687002182007 0.03632888197898865
+  0.49166208505630493 0.06087132543325424
+  0.49171400070190430 0.04841339960694313
+  0.49171534180641174 0.04839951172471046
+  0.49045535922050476 0.06087132915854454
+  4 5
+  6 4
+  0 1 2 3
+  )";
+  fill_input_from_string(&in, spec);
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL);
+  EXPECT_EQ(out->verts_len, 7);
+  EXPECT_EQ(out->edges_len, 12);
+  EXPECT_EQ(out->faces_len, 6);
+  free_spec_arrays(&in);
+}
+
+TEST(delaunay, ClosePts2)
+{
+  CDT_input in;
+  CDT_result *out;
+  const char *spec = R"(6 1 1
+  -0.17878936231136322 -0.44374340772628784
+  -0.17871695756912231 -0.45601493120193481
+  -0.17544384300708771 -0.45601493120193481
+  -0.17537136375904083 -0.44374340772628784
+  -0.17544738948345184 -0.45602506399154663
+  -0.17872454226016998 -0.45472940802574158
+  4 5
+  0 1 2 3
+  )";
+  fill_input_from_string(&in, spec);
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL);
+  EXPECT_EQ(out->verts_len, 6);
+  EXPECT_EQ(out->edges_len, 10);
+  EXPECT_EQ(out->faces_len, 5);
+  free_spec_arrays(&in);
+}
+#endif
+
+#if DO_RANDOM_TESTS
 enum {
   RANDOM_PTS,
   RANDOM_SEGS,
@@ -780,17 +997,17 @@ static void rand_delaunay_test(int test_kind,
 
 TEST(delaunay, randompts)
 {
-  rand_delaunay_test(RANDOM_PTS, 7, 1, CDT_FULL);
+  rand_delaunay_test(RANDOM_PTS, 7, 100, CDT_FULL);
 }
 
 TEST(delaunay, randomsegs)
 {
-  rand_delaunay_test(RANDOM_SEGS, 7, 1, CDT_FULL);
+  rand_delaunay_test(RANDOM_SEGS, 7, 100, CDT_FULL);
 }
 
 TEST(delaunay, randompoly)
 {
-  rand_delaunay_test(RANDOM_POLY, 7, 1, CDT_FULL);
+  rand_delaunay_test(RANDOM_POLY, 7, 100, CDT_FULL);
 }
 
 TEST(delaunay, randompoly_inside)
@@ -809,48 +1026,36 @@ TEST(delaunay, randompoly_validbmesh)
 }
 #endif
 
-#if 0
-/* For debugging or timing large examples.
- * The given file should have one point per line, as a space-separated pair of floats
+#if DO_FILE_TESTS
+/* For timing large examples of points only.
+ * The given file should have one point per line, as a space-separated pair of floats.
  */
 static void points_from_file_test(const char *filename)
 {
-  std::ifstream f;
-  std::string line;
-  struct XY {
-    float x;
-    float y;
-  } xy;
-  std::vector<XY> pts;
-  int i, n;
   CDT_input in;
   CDT_result *out;
   double tstart;
-  float (*p)[2];
 
-  f.open(filename);
-  while (getline(f, line)) {
-    std::istringstream iss(line);
-    iss >> xy.x >> xy.y;
-    pts.push_back(xy);
-  }
-  n = (int)pts.size();
-  fprintf(stderr, "read %d points\n", (int)pts.size());
-  p = (float (*)[2])MEM_malloc_arrayN(n, 2 * sizeof(float), "delaunay");
-  for (i = 0; i < n; i++) {
-    xy = pts[i];
-    p[i][0] = xy.x;
-    p[i][1] = xy.y;
-  }
+  fill_input_from_file(&in, filename);
   tstart = PIL_check_seconds_timer();
-  fill_input_verts(&in, p, n);
   out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL);
   fprintf(stderr, "time to triangulate=%f seconds\n", PIL_check_seconds_timer() - tstart);
   BLI_delaunay_2d_cdt_free(out);
-  MEM_freeN(p);
+  free_spec_arrays(&in);
+}
+
+TEST(delaunay, debug)
+{
+  CDT_input in;
+  CDT_result *out;
+  fill_input_from_file(&in, "/tmp/cdtinput.txt");
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_FULL);
+  BLI_delaunay_2d_cdt_free(out);
+  free_spec_arrays(&in);
 }
 
-#  define POINTFILEROOT "/tmp/"
+#  if 0
+#    define POINTFILEROOT "/tmp/"
 
 TEST(delaunay, terrain1)
 {
@@ -866,4 +1071,5 @@ TEST(delaunay, terrain3)
 {
   points_from_file_test(POINTFILEROOT "points3.txt");
 }
+#  endif
 #endif