Made BLI_delaunay_2d_cdt_calc better at tiny feature elimination.

The 'random' unit tests and some examples from the new boolean code
triggered asserts and crashes. This fixes those.
There is a new flag in the input that optionally disables a pass
over input to snap segment edges to other segments.

4 files changed, 1489 insertions, 468 deletions

diff --git a/source/blender/blenlib/BLI_delaunay_2d.h b/source/blender/blenlib/BLI_delaunay_2d.h
index fe81de5fc5e..9d853dd9509 100644
--- a/source/blender/blenlib/BLI_delaunay_2d.h
+++ b/source/blender/blenlib/BLI_delaunay_2d.h
@@ -104,6 +104,11 @@
  * If zero is supplied for epsilon, an internal value of 1e-8 used
  * instead, since this code will not work correctly if it is not allowed
  * to merge "too near" vertices.
+ *
+ * Normally, if epsilon is non-zero, there is an "input modify" pass which
+ * checks to see if some vertices are within epsilon of other edges, and
+ * snapping them to those edges if so. You can skip this pass by setting
+ * skip_input_modify to true. (This is also useful in some unit tests.)
  */
 typedef struct CDT_input {
   int verts_len;
@@ -115,6 +120,7 @@ typedef struct CDT_input {
   int *faces_start_table;
   int *faces_len_table;
   float epsilon;
+  bool skip_input_modify;
 } CDT_input;
 
 /**
diff --git a/source/blender/blenlib/intern/delaunay_2d.c b/source/blender/blenlib/intern/delaunay_2d.c
index c222ed06aaa..235ea6a4da4 100644
--- a/source/blender/blenlib/intern/delaunay_2d.c
+++ b/source/blender/blenlib/intern/delaunay_2d.c
@@ -32,7 +32,7 @@
 #include "BLI_delaunay_2d.h"
 
 /* Uncomment this define to get helpful debugging functions etc. defined. */
-#define DEBUG_CDT
+// #define DEBUG_CDT
 
 struct CDTEdge;
 struct CDTFace;
@@ -52,6 +52,7 @@ typedef struct CDTVert {
   LinkNode *input_ids; /* List of corresponding vertex input ids. */
   int index;           /* Index into array that cdt keeps. */
   int merge_to_index;  /* Index of a CDTVert that this has merged to. -1 if no merge. */
+  int visit_index;     /* Which visit epoch has this been seen. */
 } CDTVert;
 
 typedef struct CDTEdge {
@@ -61,7 +62,7 @@ typedef struct CDTEdge {
 } CDTEdge;
 
 typedef struct CDTFace {
-  SymEdge *symedge;    /* A symedge in face; only used during output. */
+  SymEdge *symedge;    /* A symedge in face; only used during output, so only valid then. */
   LinkNode *input_ids; /* List of input face ids that this is part of. */
   int visit_index;     /* Which visit epoch has this been seen. */
   bool deleted;        /* Marks this face no longer used. */
@@ -69,35 +70,30 @@ typedef struct CDTFace {
 } CDTFace;
 
 typedef struct CDT_state {
-  LinkNode *edges;
-  LinkNode *faces;
-  CDTFace *outer_face;
-  CDTVert **vert_array;
-  int vert_array_len;
-  int vert_array_len_alloc;
-  int input_vert_tot;
-  double minx;
-  double miny;
-  double maxx;
-  double maxy;
-  double margin;
-  int visit_count;
-  int face_edge_offset;
-  MemArena *arena;
-  BLI_mempool *listpool;
-  double epsilon;
-  double epsilon_squared;
-  bool output_prepared;
+  LinkNode *edges;          /* List of CDTEdge pointer. */
+  LinkNode *faces;          /* List of CDTFace pointer. */
+  CDTFace *outer_face;      /* Which CDTFace is the outer face. */
+  CDTVert **vert_array;     /* Array of CDTVert pointer, grows. */
+  int vert_array_len;       /* Current length of vert_array. */
+  int vert_array_len_alloc; /* Allocated length of vert_array. */
+  int input_vert_tot;       /* How many verts were in input (will be first in vert_array). */
+  double minx;              /* Used for debug drawing. */
+  double miny;              /* Used for debug drawing. */
+  double maxx;              /* Used for debug drawing. */
+  double maxy;              /* Used for debug drawing. */
+  double margin;            /* Used for debug drawing. */
+  int visit_count; /* Used for visiting things without having to initialized their visit fields. */
+  int face_edge_offset;  /* Input edge id where we start numbering the face edges. */
+  MemArena *arena;       /* Most allocations are done from here, so can free all at once at end. */
+  BLI_mempool *listpool; /* Allocations of ListNodes done from this pool. */
+  double epsilon;        /* The user-specified nearness limit. */
+  double epsilon_squared; /* Square of epsilon. */
+  bool output_prepared;   /* Set after the mesh has been modified for output (may not be all
+                             triangles now). */
 } CDT_state;
 
 #define DLNY_ARENASIZE 1 << 14
 
-/**
- * This margin means that will only be a 1 degree possible
- * concavity on outside if remove all border touching triangles.
- */
-#define DLNY_MARGIN_PCT 2000.0
-
 #ifdef DEBUG_CDT
 #  ifdef __GNUC__
 #    define ATTU __attribute__((unused))
@@ -106,14 +102,22 @@ typedef struct CDT_state {
 #  endif
 #  define F2(p) p[0], p[1]
 #  define F3(p) p[0], p[1], p[2]
+struct CrossData;
 ATTU static void dump_se(const SymEdge *se, const char *lab);
+ATTU static void dump_se_short(const SymEdge *se, const char *lab);
 ATTU static void dump_v(const CDTVert *v, const char *lab);
 ATTU static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit);
 ATTU static void dump_id_list(const LinkNode *id_list, const char *lab);
+ATTU static void dump_cross_data(struct CrossData *cd, const char *lab);
 ATTU static void dump_cdt(const CDT_state *cdt, const char *lab);
 ATTU static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab);
 ATTU static void cdt_draw(CDT_state *cdt, const char *lab);
+ATTU static void cdt_draw_region(
+    CDT_state *cdt, const char *lab, double minx, double miny, double maxx, double maxy);
+
 ATTU static void cdt_draw_vertex_region(CDT_state *cdt, int v, double dist, const char *lab);
+ATTU static void cdt_draw_edge_region(
+    CDT_state *cdt, int v1, int v2, double dist, const char *lab);
 ATTU static void write_cdt_input_to_file(const CDT_input *inp);
 ATTU static void validate_cdt(CDT_state *cdt,
                               bool check_all_tris,
@@ -140,8 +144,7 @@ BLI_INLINE SymEdge *prev(const SymEdge *se)
 /**
  * Return true if a -- b -- c are in that order, assuming they are on a straight line according to
  * orient2d and we know the order is either `abc` or `bac`.
- * This means `ab . ac` and `bc . ac` must both be non-negative.
- */
+ * This means `ab . ac` and `bc . ac` must both be non-negative.  */
 static bool in_line(const double a[2], const double b[2], const double c[2])
 {
   double ab[2], bc[2], ac[2];
@@ -187,6 +190,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->merge_to_index = -1;
+  v->visit_index = 0;
   cdt->vert_array[cdt->vert_array_len++] = v;
   return v;
 }
@@ -290,20 +294,38 @@ BLI_INLINE bool is_original_vert(const CDTVert *v, CDT_state *cdt)
   return (v->index < cdt->input_vert_tot);
 }
 
-/** Is there already an edge between a and b? */
-static bool exists_edge(const CDTVert *a, const CDTVert *b)
+/** Return the Symedge that goes from v1 to v2, if it exists, else return NULL. */
+static SymEdge *find_symedge_between_verts(const CDTVert *v1, const CDTVert *v2)
 {
-  SymEdge *se, *ss;
-  se = a->symedge;
-  if (se->next->vert == b) {
-    return true;
-  }
-  for (ss = se->rot; ss != se; ss = ss->rot) {
-    if (ss->next->vert == b) {
-      return true;
+  SymEdge *tstart, *t;
+
+  t = tstart = v1->symedge;
+  do {
+    if (t->next->vert == v2) {
+      return t;
     }
-  }
-  return false;
+  } while ((t = t->rot) != tstart);
+  return NULL;
+}
+
+/** Return the SymEdge attached to v that has face f, if it exists, else return NULL. */
+static SymEdge *find_symedge_with_face(const CDTVert *v, const CDTFace *f)
+{
+  SymEdge *tstart, *t;
+
+  t = tstart = v->symedge;
+  do {
+    if (t->face == f) {
+      return t;
+    }
+  } while ((t = t->rot) != tstart);
+  return NULL;
+}
+
+/** Is there already an edge between a and b? */
+static inline bool exists_edge(const CDTVert *v1, const CDTVert *v2)
+{
+  return find_symedge_between_verts(v1, v2) != NULL;
 }
 
 /** Is the vertex v incident on face f? */
@@ -533,7 +555,7 @@ static void delete_edge(CDT_state *cdt, SymEdge *e)
   }
 }
 
-static CDT_state *new_cdt_init(const CDT_input *in)
+static CDT_state *cdt_init(const CDT_input *in)
 {
   int i;
   MemArena *arena = BLI_memarena_new(DLNY_ARENASIZE, __func__);
@@ -1121,496 +1143,711 @@ static void re_delaunay_triangulate(CDT_state *cdt, SymEdge *se)
   }
 }
 
+static double tri_orient(const SymEdge *t)
+{
+  return orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co);
+}
+
 /**
- * Add a constrained edge between v1 and v2 to cdt structure.
- * This may result in a number of #CDTEdges created, due to intersections
- * and partial overlaps with existing cdt vertices and edges.
- * Each created #CDTEdge will have input_id added to its input_ids list.
+ * The CrossData struct gives defines either an endpoint or an intermediate point
+ * in the path we will take to insert an edge constraint.
+ * Each such point will either be
+ * (a) a vertex or
+ * (b) a fraction lambda (0 < lambda < 1) along some SymEdge.]
  *
- * If \a r_edges is not NULL, the #CDTEdges generated or found that go from
- * v1 to v2 are put into that linked list, in order.
+ * In general, lambda=0 indicates case a and lambda != 0 indicates case be.
+ * The 'in' edge gives the destination attachment point of a diagonal from the previous crossing,
+ * and the 'out' edge gives the origin attachment point of a diagonal to the next crossing.
+ * But in some cases, 'in' and 'out' are undefined or not needed, and will be NULL.
  *
- * Assumes that #BLI_constrained_delaunay_get_output has not been called yet.
+ * For case (a), 'vert' will be the vertex, and lambda will be 0, and 'in' will be the SymEdge from
+ * 'vert' that has as face the one that you go through to get to this vertex. If you go exactly
+ * along an edge then we set 'in' to NULL, since it won't be needed. The first crossing will have
+ * 'in' = NULL. We set 'out' to the SymEdge that has the face we go though to get to the next
+ * crossing, or, if the next crossing is a case (a), then it is the edge that goes to that next
+ * vertex. 'out' wlll be NULL for the last one.
+ *
+ * For case (b), vert will be NULL at first, and later filled in with the created split vertex,
+ * and 'in' will be the SymEdge that we go through, and lambda will be between 0 and 1,
+ * the fraction from in's vert to in->next's vert to put the split vertex.
+ * 'out' is not needed in this case, since the attachment point will be the sym of the first
+ * half of the split edge.
  */
-static void add_edge_constraint(
-    CDT_state *cdt, CDTVert *v1, CDTVert *v2, int input_id, LinkNode **r_edges)
+typedef struct CrossData {
+  double lambda;
+  CDTVert *vert;
+  SymEdge *in;
+  SymEdge *out;
+} CrossData;
+
+static bool get_next_crossing_from_vert(CDT_state *cdt,
+                                        CrossData *cd,
+                                        CrossData *cd_next,
+                                        const CDTVert *v2);
+
+/**
+ * As part of finding crossings, we found a case where the next crossing goes through vert v.
+ * If it came from a previous vert in cd, then cd_out is the edge that leads from that to v.
+ * Else cd_out can be NULL, because it won't be used.
+ * Set *cd_next to indicate this. We can set 'in' but not 'out'.  We can set the 'out' of the
+ * current cd.
+ */
+static void fill_crossdata_for_through_vert(CDTVert *v,
+                                            SymEdge *cd_out,
+                                            CrossData *cd,
+                                            CrossData *cd_next)
 {
-  CDTVert *va, *vb, *vc;
-  SymEdge *vse1;
-#ifdef DEBUG_CDT
-  SymEdge *vse2;
-#endif
-  SymEdge *t, *tstart, *tout, *tnext;
   SymEdge *se;
-  CDTEdge *edge;
-  int isect;
-  double orient1, orient2;
-  int i, search_count;
-  double curco[2];
-  double lambda;
-  bool done, state_through_vert;
-  LinkNodePair edge_list = {NULL, NULL};
-  typedef struct CrossData {
-    double lambda;
-    CDTVert *vert;
-    SymEdge *in;
-    SymEdge *out;
-  } CrossData;
-  CrossData cdata;
-  CrossData *crossings = NULL;
-  CrossData *cd;
-  BLI_array_staticdeclare(crossings, 128);
 #ifdef DEBUG_CDT
   int dbg_level = 0;
+
+  if (global_dbg) {
+    dbg_level = 1;
+  }
 #endif
 
-  /* Find path through structure from v1 to v2 and record how we got there in crossings.
-   * In crossings array, each CrossData is populated as follows:
-   *
-   * If ray from previous node goes through a face, not along an edge:
-   *
-   *               _  B
-   *            /    |\
-   *         - -     | \
-   *     prev........X  \
-   *       \ d       |   \C
-   *        --       |   /
-   *          \     a| b/
-   *           - -   | /
-   *               \  A
-   *
-   *   lambda = fraction of way along AB where X is.
-   *   vert = NULL initially, will later get new node that splits AB
-   *   in = a (SymEdge from A->B, whose face the ray goes through)
-   *   out = b (SymEdge from A->C, whose face the ray goes through next
-   *
-   * If the ray from the previous node goes directly to an existing vertex, say A
-   * in the previous diagram, maybe along an existing edge like d in that diagram
-   * but if prev had lambda !=0 then there may be no such edge d, then:
-   *
-   *    lambda = 0
-   *    vert = A
-   *    in = a
-   *    out = b
-   *
-   * crossings[0] will have in = NULL, and crossings[last] will have out = NULL.
-   */
-  if (r_edges) {
-    *r_edges = NULL;
+  cd_next->lambda = 0.0;
+  cd_next->vert = v;
+  cd_next->in = NULL;
+  cd_next->out = NULL;
+  if (cd->lambda == 0.0) {
+    cd->out = cd_out;
+  }
+  else {
+    /* One of the edges in the triangle with edge sym(cd->in) contains v. */
+    se = sym(cd->in);
+    if (se->vert != v) {
+      se = se->next;
+      if (se->vert != v) {
+        se = se->next;
+      }
+    }
+    BLI_assert(se->vert == v);
+    cd_next->in = se;
   }
-  vse1 = v1->symedge;
 #ifdef DEBUG_CDT
   if (dbg_level > 0) {
-    vse2 = v2->symedge;
-    if (dbg_level > 2) {
-      // dump_cdt(cdt, "before insert_segment");
-      cdt_draw(cdt, "before insert segment");
-    }
-    fprintf(stderr, "\ninsert_segment %d\n", input_id);
-    dump_v(v1, "  1");
-    dump_v(v2, "  2");
-    if (dbg_level > 1) {
-      dump_se(vse1, "  se1");
-      dump_se(vse2, "  se2");
-    }
+    dump_cross_data(cd, "cd through vert, cd");
+    dump_cross_data(cd_next, "cd_next through vert, cd");
   }
 #endif
-  if (v1 == v2) {
+}
+
+/**
+ * As part of finding crossings, we found a case where orient tests say that the next crossing
+ * is on the SymEdge t, while intersecting with the ray from curco to v2.
+ * Find the intersection point and fill in the CrossData for that point.
+ * It may turn out that when doing the intersection, we get an answer that says that
+ * this case is better handled as through-vertex case instead, so we may do that.
+ * In the latter case, we want to avoid a situation where the current crossing is on an edge
+ * and the next will be an endpoint of the same edge. When that happens, we "rewrite history"
+ * and turn the current crossing into a vert one, and then extend from there.
+ *
+ * We cannot fill cd_next's 'out' edge yet, in the case that the next one ends up being a vert
+ * case. We need to fill in cd's 'out' edge if it was a vert case.
+ */
+static void fill_crossdata_for_intersect(CDT_state *cdt,
+                                         const double *curco,
+                                         const CDTVert *v2,
+                                         SymEdge *t,
+                                         CrossData *cd,
+                                         CrossData *cd_next)
+{
+  CDTVert *va, *vb, *vc;
+  double lambda, mu, len_ab;
+  SymEdge *se_vcva, *se_vcvb;
+  int isect;
 #ifdef DEBUG_CDT
-    if (dbg_level > 0) {
-      fprintf(stderr, "segment between same vertices, ignored\n");
-    }
+  int dbg_level = 0;
 #endif
-    return;
+
+  va = t->vert;
+  vb = t->next->vert;
+  vc = t->next->next->vert;
+  se_vcvb = sym(t->next);
+  se_vcva = t->next->next;
+  BLI_assert(se_vcva->vert == vc && se_vcva->next->vert == va);
+  BLI_assert(se_vcvb->vert == vc && se_vcvb->next->vert == vb);
+  isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, curco, v2->co, &lambda, &mu);
+#ifdef DEBUG_CDT
+  if (dbg_level > 0) {
+    double co[2];
+    fprintf(stderr, "crossdata for intersect gets lambda=%.17g, mu=%.17g\n", lambda, mu);
+    fprintf(stderr,
+            "isect=%s\n",
+            isect == 2 ? "cross" : (isect == 1 ? "exact" : (isect == 0 ? "none" : "colinear")));
+    fprintf(stderr,
+            "va=v%d=(%g,%g), vb=v%d=(%g,%g), vc=v%d, curco=(%g,%g), v2=(%g,%g)\n",
+            va->index,
+            F2(va->co),
+            vb->index,
+            F2(vb->co),
+            vc->index,
+            F2(curco),
+            F2(v2->co));
+    dump_se_short(se_vcva, "vcva=");
+    dump_se_short(se_vcvb, " vcvb=");
+    interp_v2_v2v2_db(co, va->co, vb->co, lambda);
+    fprintf(stderr, "\nco=(%.17g,%.17g)\n", F2(co));
   }
-  /* Quick check for special case where segment is one of edges from v1. */
-  tstart = t = vse1;
-  done = false;
-  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);
-      cdata.in = t->next->rot;
-      cdata.out = NULL;
-      cdata.lambda = 0.0;
-      cdata.vert = v2;
-      BLI_array_append(crossings, cdata);
+#endif
+  switch (isect) {
+    case ISECT_LINE_LINE_CROSS:
+      len_ab = len_v2v2_db(va->co, vb->co);
 #ifdef DEBUG_CDT
       if (dbg_level > 0) {
-        fprintf(stderr, "special one segment case\n");
-        dump_se(t, "  ");
+        fprintf(stderr,
+                "len_ab=%g, near a=%g, near b=%g\n",
+                len_ab,
+                lambda * len_ab,
+                (1.0 - lambda) * len_ab);
       }
 #endif
-      done = true;
+      if (lambda * len_ab <= cdt->epsilon) {
+        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
+      }
+      else if ((1.0 - lambda) * len_ab <= cdt->epsilon) {
+        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
+      }
+      else {
+        *cd_next = (CrossData){lambda, NULL, t, NULL};
+        if (cd->lambda == 0.0) {
+          cd->out = se_vcva;
+        }
+      }
+      break;
+    case ISECT_LINE_LINE_EXACT:
+      if (lambda == 0.0) {
+        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
+      }
+      else if (lambda == 1.0) {
+        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
+      }
+      else {
+        *cd_next = (CrossData){lambda, NULL, t, NULL};
+        if (cd->lambda == 0.0) {
+          cd->out = se_vcva;
+        }
+      }
+      break;
+    case ISECT_LINE_LINE_NONE:
+      /* It should be very near one end or other of segment. */
+      if (lambda <= 0.5) {
+        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
+      }
+      else {
+        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
+      }
+      break;
+    case ISECT_LINE_LINE_COLINEAR:
+      if (len_squared_v2v2_db(va->co, v2->co) <= len_squared_v2v2_db(vb->co, v2->co)) {
+        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
+      }
+      else {
+        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
+      }
       break;
+  }
+}
+
+/**
+ * As part of finding the crossings of a ray to v2, find the next crossing after 'cd', assuming
+ * 'cd' represents a crossing that goes through a vertex.
+ *
+ * We do a rotational scan around cd's vertex, looking for the triangle where the ray from cd->vert
+ * to v2 goes between the two arms from cd->vert, or where it goes along one of the edges.
+ */
+static bool get_next_crossing_from_vert(CDT_state *cdt,
+                                        CrossData *cd,
+                                        CrossData *cd_next,
+                                        const CDTVert *v2)
+{
+  SymEdge *t, *tstart;
+  CDTVert *vcur, *va, *vb;
+  double orient1, orient2;
+  bool ok;
+#ifdef DEBUG_CDT
+  int dbg_level = 0;
+
+  if (dbg_level > 0) {
+    fprintf(stderr, "\nget_next_crossing_from_vert\n");
+    dump_v(cd->vert, "  cd->vert");
+  }
+#endif
+
+  t = tstart = cd->vert->symedge;
+  vcur = cd->vert;
+  ok = false;
+  do {
+    /*
+     * The ray from vcur to v2 has to go either between two successive
+     * edges around vcur or exactly along them. This time through the
+     * loop, check to see if the ray goes along vcur-va
+     * or between vcur-va and vcur-vb, where va is the end of t
+     * and vb is the next vertex (on the next rot edge around vcur, but
+     * should also be the next vert of triangle starting with vcur-va.
+     */
+    if (t->face != cdt->outer_face && tri_orient(t) < 0.0) {
+      fprintf(stderr, "BAD TRI orientation %g\n", tri_orient(t)); /* TODO: handle this */
+#ifdef DEBUG_CDT
+      dump_se_cycle(t, "top of ccw scan loop", 4);
+#endif
     }
-    t = t->rot;
-  } while (t != tstart);
-  if (!done) {
-    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.
-       */
+    va = t->next->vert;
+    vb = t->next->next->vert;
+    orient1 = orient2d(t->vert->co, va->co, v2->co);
+#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, "orient1=%g\n", orient1);
+    }
+#endif
+    if (orient1 == 0.0 && in_line(vcur->co, va->co, v2->co)) {
 #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);
+        fprintf(stderr, "ray goes through va\n");
       }
 #endif
-      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) {
+      fill_crossdata_for_through_vert(va, t, cd, cd_next);
+      ok = true;
+      break;
+    }
+    else if (t->face != cdt->outer_face) {
+      orient2 = orient2d(vcur->co, vb->co, v2->co);
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "found v2, so done\n");
-          }
+      if (dbg_level > 1) {
+        fprintf(stderr, "orient2=%g\n", orient2);
+      }
 #endif
-          done = true;
+      /* Don't handle orient2 == 0.0 case here: next rotation will get it. */
+      if (orient1 > 0.0 && orient2 < 0.0) {
+#ifdef DEBUG_CDT
+        if (dbg_level > 1) {
+          fprintf(stderr, "segment intersects\n");
         }
-        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;
-            orient1 = orient2d(t->vert->co, va->co, v2->co);
+#endif
+        t = t->next;
+        fill_crossdata_for_intersect(cdt, vcur->co, v2, t, cd, cd_next);
+        ok = true;
+        break;
+      }
+    }
+  } while ((t = t->rot) != tstart);
 #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, "orient1=%g\n", orient1);
-            }
+  if (!ok) {
+    /* Didn't find the exit! Shouldn't happen. */
+    fprintf(stderr, "shouldn't happen: can't find next crossing from vert\n");
+  }
 #endif
-            if (orient1 == 0.0 && in_line(t->vert->co, va->co, v2->co)) {
+  return ok;
+}
+
+/**
+ * As part of finding the crossings of a ray to v2, find the next crossing after 'cd', assuming
+ * 'cd' represents a crossing that goes through a an edge, not at either end of that edge.
+ *
+ * We have the triangle vb-va-vc, where va and vb are the split edge and vc is the third vertex on
+ * that new side of the edge (should be closer to v2). The next crossing should be through vc or
+ * intersecting vb-vc or va-vc.
+ */
+static void get_next_crossing_from_edge(CDT_state *cdt,
+                                        CrossData *cd,
+                                        CrossData *cd_next,
+                                        const CDTVert *v2)
+{
+  double curco[2];
+  double orient;
+  CDTVert *va, *vb, *vc;
+  SymEdge *se_ac;
 #ifdef DEBUG_CDT
-              if (dbg_level > 1) {
-                fprintf(stderr, "ray goes through va\n");
-              }
+  int dbg_level = 0;
+
+  if (dbg_level > 0) {
+    fprintf(stderr, "\nget_next_crossing_from_edge\n");
+    fprintf(stderr, "  lambda=%.17g\n", cd->lambda);
+    dump_se_short(cd->in, "  cd->in");
+  }
 #endif
-              state_through_vert = true;
-              tout = t;
-              t = t->next;
-              break;
-            }
-            else if (t->face != cdt->outer_face) {
-              orient2 = orient2d(t->vert->co, vb->co, v2->co);
+
+  va = cd->in->vert;
+  vb = cd->in->next->vert;
+  interp_v2_v2v2_db(curco, va->co, vb->co, cd->lambda);
 #ifdef DEBUG_CDT
-              if (dbg_level > 1) {
-                fprintf(stderr, "orient2=%g\n", orient2);
-              }
+  if (dbg_level > 0) {
+    fprintf(stderr, "  curco=(%.17g,%.17g)\n", F2(curco));
+  }
 #endif
-              if (orient2 == 0.0 && in_line(t->vert->co, vb->co, v2->co)) {
+  se_ac = sym(cd->in)->next;
+  vc = se_ac->next->vert;
+  orient = orient2d(curco, v2->co, vc->co);
 #ifdef DEBUG_CDT
-                if (dbg_level > 1) {
-                  fprintf(stderr, "ray goes through vb\n");
-                }
+  if (dbg_level > 1) {
+    fprintf(stderr, "now searching with third vertex ");
+    dump_v(vc, "vc");
+    fprintf(stderr, "orient2d(cur, v2, vc) = %g\n", orient);
+  }
 #endif
-                state_through_vert = true;
-                t = t->next->next;
-                tout = sym(t);
-                break;
-              }
-              else if (orient1 > 0.0 && orient2 < 0.0) {
+  if (orient < 0.0) {
 #ifdef DEBUG_CDT
-                if (dbg_level > 1) {
-                  fprintf(stderr, "segment intersects\n");
-                }
+    if (dbg_level > 1) {
+      fprintf(stderr, "curco--v2 intersects vb--vc\n");
+    }
 #endif
-                state_through_vert = false;
-                tout = t;
-                t = t->next;
-                break;
-              }
-            }
-            t = t->rot;
+    fill_crossdata_for_intersect(cdt, curco, v2, se_ac->next, cd, cd_next);
+  }
+  else if (orient > 0.0) {
 #ifdef DEBUG_CDT
-            if (dbg_level > 1) {
-              dump_se_cycle(t, "next rot tri", 4);
-            }
+    if (dbg_level > 1) {
+      fprintf(stderr, "curco--v2 intersects va--vc\n");
+    }
 #endif
-          } while (t != tstart);
-          BLI_assert(tout != NULL);
-          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;
-        /* Get curco; cdata should have data from last time through the loop still. */
-        if (cdata.lambda == 0.0) {
-          copy_v2_v2_db(curco, cdata.vert->co);
-        }
-        else {
-          interp_v2_v2v2_db(curco, cdata.in->vert->co, cdata.in->next->vert->co, lambda);
-        }
+    fill_crossdata_for_intersect(cdt, curco, v2, se_ac, cd, cd_next);
+  }
+  else {
 #ifdef DEBUG_CDT
-        if (dbg_level > 1) {
-          fprintf(stderr, "through edge case, curco=(%f,%f)\n", F2(curco));
-          dump_v(va, " va");
-          dump_v(vb, " vb");
-        }
+    if (dbg_level > 1) {
+      fprintf(stderr, "orient==0 case, so going through or to vc\n");
+    }
 #endif
-        isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, curco, v2->co, &lambda, NULL);
-        if (isect == ISECT_LINE_LINE_NONE || isect == ISECT_LINE_LINE_EXACT) {
-          /* The orient tests say that there is an intersection between
-           * va and vb, but the inexact intersection routine has either put the
-           * intersection exactly on one of the endpoints or just outside
-           * one of them.
-           * Or this is an exact intersection at one of the curco / v2 ends.
-           * If lambda is outside of range, move the intersection to somewhere
-           * just inside the segment.
-           * Could also snap to an endpoint and redo this as a "through vert"
-           * case, but the short edge will be cleaned up later and this seems
-           * less risky to get into "impossible" cases.
-           */
-          if (lambda <= 0.0) {
-            lambda = 4.0 * (double)FLT_EPSILON;
-          }
-          else if (lambda >= 1.0) {
-            lambda = 1.0 - 4.0 * (double)FLT_EPSILON;
-          }
-        }
-        if (isect == ISECT_LINE_LINE_COLINEAR) {
+    *cd_next = (CrossData){0.0, vc, se_ac->next, NULL};
+  }
+}
+
+/**
+ * Add a constrained edge between v1 and v2 to cdt structure.
+ * This may result in a number of #CDTEdges created, due to intersections
+ * and partial overlaps with existing cdt vertices and edges.
+ * Each created #CDTEdge will have input_id added to its input_ids list.
+ *
+ * If \a r_edges is not NULL, the #CDTEdges generated or found that go from
+ * v1 to v2 are put into that linked list, in order.
+ *
+ * Assumes that #BLI_constrained_delaunay_get_output has not been called yet.
+ */
+static void add_edge_constraint(
+    CDT_state *cdt, CDTVert *v1, CDTVert *v2, int input_id, LinkNode **r_edges)
+{
+  SymEdge *t, *se, *tstart, *tnext;
+  int i, j, n, visit;
+  bool ok;
+  CrossData *crossings = NULL;
+  CrossData *cd, *cd_prev, *cd_next;
+  CDTVert *v;
+  CDTEdge *edge;
+  LinkNodePair edge_list = {NULL, NULL};
+  BLI_array_staticdeclare(crossings, 128);
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "intersect is collinear, treating as through vert\n");
-          }
-          dump_v(va, "va");
+  int dbg_level = 0;
+
+  if (dbg_level > 0) {
+    fprintf(stderr, "\nADD_EDGE_CONSTRAINT\n");
+    dump_v(v1, "  1");
+    dump_v(v2, "  2");
+  }
 #endif
-          state_through_vert = true;
-          continue;
-        }
+
+  if (r_edges) {
+    *r_edges = NULL;
+  }
+
+  /*
+   * Handle two special cases first:
+   * 1) The two end vertices are the same (can happen because of merging).
+   * 2) There is already an edge between v1 and v2.
+   */
+  if (v1 == v2) {
+    return;
+  }
+  if ((t = find_symedge_between_verts(v1, v2)) != NULL) {
 #ifdef DEBUG_CDT
-        interp_v2_v2v2_db(curco, va->co, vb->co, lambda);
-        if (dbg_level > 0) {
-          fprintf(stderr, "intersect point at lambda=%.17g along va--vb\n", lambda);
-          fprintf(stderr, "which is (%g,%g)\n", F2(curco));
-          if (dbg_level > 1) {
-            dump_v(va, " va");
-            dump_v(vb, " vb");
-          }
-        }
+    if (dbg_level > 0) {
+      fprintf(stderr, "segment already there\n");
+    }
 #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);
+    add_to_input_ids(&t->edge->input_ids, input_id, cdt);
+    if (r_edges != NULL) {
+      BLI_linklist_append_pool(&edge_list, t->edge, cdt->listpool);
+      *r_edges = edge_list.list;
+    }
+    return;
+  }
+
+  /*
+   * Fill crossings array with CrossData points for intersection path from v1 to v2.
+   *
+   * At every point, the crossings array has the path so far, except that
+   * the .out field of the last element of it may not be known yet -- if that
+   * last element is a vertex, then we won't know the output edge until we
+   * find the next crossing.
+   *
+   * To protect against infinite loops, we keep track of which vertices
+   * we have visited by setting their visit_index to a new visit epoch.
+   *
+   * We check a special case first: where the segment is already there in
+   * one hop. Saves a bunch of orient2d tests in that common case.
+   */
+  visit = ++cdt->visit_count;
+  BLI_array_grow_one(crossings);
+  crossings[0] = (CrossData){0.0, v1, NULL, NULL};
+  while (!((n = BLI_array_len(crossings)) > 0 && crossings[n - 1].vert == v2)) {
+    BLI_array_grow_one(crossings);
+    cd = &crossings[n - 1];
+    cd_next = &crossings[n];
+    if (crossings[n - 1].lambda == 0.0) {
+      ok = get_next_crossing_from_vert(cdt, cd, cd_next, v2);
+    }
+    else {
+      get_next_crossing_from_edge(cdt, cd, cd_next, v2);
+    }
+    if (!ok || BLI_array_len(crossings) == 100000) {
+      /* Shouldn't happen but if does, just bail out. */
 #ifdef DEBUG_CDT
-        if (dbg_level > 0) {
-          dump_se_cycle(tout, "next search tri", 4);
-          dump_se(tout, "tout");
-        }
+      fprintf(stderr, "FAILURE adding segment, bailing out\n");
 #endif
-        /* 'tout' is 'symedge' from 'va' to third vertex, 'vc'. */
-        BLI_assert(tout->vert == va);
-        vc = tout->next->vert;
-        orient1 = orient2d(curco, v2->co, vc->co);
+      BLI_array_free(crossings);
+      return;
+    }
 #ifdef DEBUG_CDT
-        if (dbg_level > 1) {
-          fprintf(stderr, "now searching with third vertex ");
-          dump_v(vc, "vc");
-          fprintf(stderr, "orient2d(vcur, v2, vc) = %g\n", orient1);
-        }
+    if (dbg_level > 1) {
+      fprintf(stderr, "crossings[%d]: ", n);
+      dump_cross_data(&crossings[n], "");
+    }
 #endif
-        if (orient1 < 0.0) {
-          /* vcur--v2 should intersect vb--vc. */
+    if (crossings[n].lambda == 0.0) {
+      if (crossings[n].vert->visit_index == visit) {
+        fprintf(stderr, "WHOOPS, REVISIT. Bailing out.\n"); /*TODO: desperation search here. */
+        BLI_array_free(crossings);
+        return;
+      }
+      crossings[n].vert->visit_index = visit;
+    }
+  }
+
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "curco--v2 intersects vb--vc\n");
+  if (dbg_level > 0) {
+    fprintf(stderr, "\ncrossings found\n");
+    for (i = 0; i < BLI_array_len(crossings); i++) {
+      fprintf(stderr, "%d: ", i);
+      dump_cross_data(&crossings[i], "");
+      if (crossings[i].lambda == 0.0) {
+        if (i == 0 || crossings[i - 1].lambda == 0.0) {
+          BLI_assert(crossings[i].in == NULL);
+        }
+        else {
+          BLI_assert(crossings[i].in != NULL && crossings[i].in->vert == crossings[i].vert);
+          BLI_assert(crossings[i].in->face == sym(crossings[i - 1].in)->face);
+        }
+        if (i == BLI_array_len(crossings) - 1) {
+          BLI_assert(crossings[i].vert == v2);
+          BLI_assert(crossings[i].out == NULL);
+        }
+        else {
+          BLI_assert(crossings[i].out->vert == crossings[i].vert);
+          if (crossings[i + 1].lambda == 0.0) {
+            BLI_assert(crossings[i].out->next->vert == crossings[i + 1].vert);
           }
+          else {
+            BLI_assert(crossings[i].out->face == crossings[i + 1].in->face);
+          }
+        }
+      }
+      else {
+        if (i > 0 && crossings[i - 1].lambda == 0.0) {
+          BLI_assert(crossings[i].in->face == crossings[i - 1].out->face);
+        }
+        BLI_assert(crossings[i].out == NULL);
+      }
+    }
+  }
 #endif
-          t = tout->next;
-          state_through_vert = false;
+
+  /*
+   * Postprocess crossings.
+   * Some crossings may have an intersection crossing followed
+   * by a vertex crossing that is on the same edge that was just
+   * intersected. We prefer to go directly from the previous
+   * crossing directly to the vertex. This may chain backwards.
+   *
+   * This loop marks certain crossings as "deleted", by setting
+   * their lambdas to -1.0.
+   */
+  for (i = 2; i < BLI_array_len(crossings); i++) {
+    cd = &crossings[i];
+    if (cd->lambda == 0.0) {
+      v = cd->vert;
+      for (j = i - 1; j > 0; j--) {
+        cd_prev = &crossings[j];
+        if ((cd_prev->lambda == 0.0 && cd_prev->vert != v) ||
+            (cd_prev->lambda != 0.0 && cd_prev->in->vert != v && cd_prev->in->next->vert != v)) {
+          break;
         }
-        else if (orient1 > 0.0) {
-          /* vcur--v2 should intersect va--vc. */
+        else {
+          cd_prev->lambda = -1.0; /* Mark cd_prev as 'deleted'. */
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "curco--v2 intersects va--vc\n");
+          if (dbg_level > 0) {
+            fprintf(stderr, "deleted crossing %d\n", j);
           }
 #endif
-          t = tout;
-          state_through_vert = false;
         }
-        else if (orient1 == 0.0) {
+      }
+      if (j < i - 1) {
+        /* Some crossings were deleted. Fix the in and out edges across gap. */
+        cd_prev = &crossings[j];
+        if (cd_prev->lambda == 0.0) {
+          se = find_symedge_between_verts(cd_prev->vert, v);
+          if (se == NULL) {
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "ccw==0 case, so going through or to vc\n");
-          }
+            fprintf(stderr, "FAILURE(a) in delete crossings, bailing out.\n");
 #endif
-          t = tout->next;
-          state_through_vert = true;
+            BLI_array_free(crossings);
+            return;
+          }
+          cd_prev->out = se;
+          cd->in = NULL;
         }
         else {
+          se = find_symedge_with_face(v, sym(cd_prev->in)->face);
+          if (se == NULL) {
 #ifdef DEBUG_CDT
-          fprintf(stderr, "add_edge_constraint desperation search needed\n");
+            fprintf(stderr, "FAILURE(b) in delete crossings, bailing out.\n");
 #endif
+            BLI_array_free(crossings);
+            return;
+          }
+          cd->in = se;
         }
-      }
-      if (++search_count > 1000000) {
-        fprintf(stderr, "infinite loop? bailing out\n");
-        BLI_assert(0); /* Catch these while developing. */
-        break;
+#ifdef DEBUG_CDT
+        if (dbg_level > 0) {
+          fprintf(stderr, "after deleting crossings:\n");
+          fprintf(stderr, "cross[%d]: ", j);
+          dump_cross_data(cd_prev, "");
+          fprintf(stderr, "cross[%d]: ", i);
+          dump_cross_data(cd, "");
+        }
+#endif
       }
     }
   }
+
+  /*
+   * Insert all intersection points on constrained edges.
+   */
+  for (i = 0; i < BLI_array_len(crossings); i++) {
+    cd = &crossings[i];
+    if (cd->lambda != 0.0 && cd->lambda != -1.0 && is_constrained_edge(cd->in->edge)) {
+      edge = split_edge(cdt, cd->in, cd->lambda);
+      cd->vert = edge->symedges[0].vert;
 #ifdef DEBUG_CDT
-  if (dbg_level > 0) {
-    fprintf(stderr, "Crossing info gathered:\n");
-    for (i = 0; i < BLI_array_len(crossings); i++) {
-      cd = &crossings[i];
-      fprintf(stderr, "%d:\n", i);
-      if (cd->vert != NULL) {
-        dump_v(cd->vert, "  vert: ");
-      }
-      else {
-        fprintf(stderr, "  lambda=%f along in\n", cd->lambda);
-      }
-      if (cd->in) {
-        dump_se(cd->in, "  in: ");
-      }
-      if (cd->out) {
-        dump_se(cd->out, "  out: ");
+      if (dbg_level > 1) {
+        fprintf(stderr, "insert vert for crossing %d: ", i);
+        dump_v(cd->vert, "inserted");
       }
+#endif
     }
   }
-#endif
 
-  if (BLI_array_len(crossings) == 2) {
-    /* For speed, handle special case of segment must have already been there. */
-    se = crossings[1].in;
-    if (se->next->vert != v1) {
-      se = prev(se);
-    }
-    BLI_assert(se->vert == v1 || se->next->vert == v1);
+  /*
+   * Remove any crossed, non-intersected edges.
+   */
+  for (i = 0; i < BLI_array_len(crossings); i++) {
+    cd = &crossings[i];
+    if (cd->lambda != 0.0 && cd->lambda != -1.0 && !is_constrained_edge(cd->in->edge)) {
+      delete_edge(cdt, cd->in);
 #ifdef DEBUG_CDT
-    if (dbg_level > 0) {
-      fprintf(stderr, "segment already there: ");
-      dump_se(se, "");
-    }
+      if (dbg_level > 1) {
+        fprintf(stderr, "delete edge for crossing %d\n", i);
+      }
 #endif
-    add_to_input_ids(&se->edge->input_ids, input_id, cdt);
-    if (r_edges != NULL) {
-      BLI_linklist_append_pool(&edge_list, se->edge, cdt->listpool);
     }
   }
-  else {
-    /* Insert all intersection points. */
-    for (i = 0; i < BLI_array_len(crossings); i++) {
-      cd = &crossings[i];
-      if (cd->lambda != 0.0 && is_constrained_edge(cd->in->edge)) {
-        edge = split_edge(cdt, cd->in, cd->lambda);
-        cd->vert = edge->symedges[0].vert;
-#ifdef DEBUG_CDT
-        if (dbg_level > 1) {
-          fprintf(stderr, "insert vert for crossing %d: ", i);
-          dump_v(cd->vert, "inserted");
-        }
-#endif
+
+  /*
+   * Insert segments for v1->v2.
+   */
+  tstart = crossings[0].out;
+  for (i = 1; i < BLI_array_len(crossings); i++) {
+    cd = &crossings[i];
+    if (cd->lambda == -1.0) {
+      continue; /* This crossing was deleted. */
+    }
+    t = tnext = NULL;
+    if (cd->lambda != 0.0) {
+      if (is_constrained_edge(cd->in->edge)) {
+        t = cd->vert->symedge;
+        tnext = sym(t)->next;
       }
     }
-
-    /* Remove any crossed, non-intersected edges. */
-    for (i = 0; i < BLI_array_len(crossings); i++) {
-      cd = &crossings[i];
-      if (cd->lambda != 0.0 && !is_constrained_edge(cd->in->edge)) {
-        delete_edge(cdt, cd->in);
+    else if (cd->lambda == 0.0) {
+      t = cd->in;
+      tnext = cd->out;
+      if (t == NULL) {
+        /* Previous non-deleted crossing must also have been a vert, and segment should exist. */
+        for (j = i - 1; j >= 0; j--) {
+          cd_prev = &crossings[j];
+          if (cd_prev->lambda != -1.0) {
+            break;
+          }
+        }
+        BLI_assert(cd_prev->lambda == 0.0);
+        BLI_assert(cd_prev->out->next->vert == cd->vert);
 #ifdef DEBUG_CDT
         if (dbg_level > 1) {
-          fprintf(stderr, "delete edge for crossing %d\n", i);
+          fprintf(stderr, "edge to crossing %d already there\n", i);
         }
 #endif
+        edge = cd_prev->out->edge;
+        add_to_input_ids(&edge->input_ids, input_id, cdt);
       }
     }
-
-    /* Insert segments for v1->v2. */
-    tstart = crossings[0].out;
-    for (i = 1; i < BLI_array_len(crossings); i++) {
-      cd = &crossings[i];
-      t = tnext = NULL;
-      if (cd->lambda != 0.0) {
-        if (is_constrained_edge(cd->in->edge)) {
-          t = cd->vert->symedge;
-          tnext = sym(t)->next;
-        }
-      }
-      else if (cd->lambda == 0.0) {
-        t = cd->in;
-        tnext = cd->out;
+    if (t != NULL) {
+#ifdef DEBUG_CDT
+      if (dbg_level > 1) {
+        fprintf(stderr, "edge to crossing %d: insert diagonal between\n", i);
+        dump_se(tstart, "  ");
+        dump_se(t, " ");
+        dump_se_cycle(tstart, "tstart", 100);
+        dump_se_cycle(t, "t", 100);
       }
-      if (t) {
+#endif
+      if (tstart->next->vert == t->vert) {
+        edge = tstart->edge;
 #ifdef DEBUG_CDT
         if (dbg_level > 1) {
-          fprintf(stderr, "insert diagonal between\n");
-          dump_se(tstart, "  ");
-          dump_se(t, " ");
-          dump_se_cycle(tstart, "tstart", 100);
-          dump_se_cycle(t, "t", 100);
+          fprintf(stderr, "already there (b)\n");
         }
 #endif
-        if (tstart->next->vert == t->vert) {
-          edge = tstart->edge;
+      }
+      else {
+        edge = add_diagonal(cdt, tstart, t);
+      }
+      add_to_input_ids(&edge->input_ids, input_id, cdt);
 #ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "already there\n");
-          }
+      if (dbg_level > 1) {
+        fprintf(stderr, "added\n");
+      }
 #endif
-        }
-        else {
-          edge = add_diagonal(cdt, tstart, t);
-        }
-        add_to_input_ids(&edge->input_ids, input_id, cdt);
+      if (r_edges != NULL) {
+        BLI_linklist_append_pool(&edge_list, edge, cdt->listpool);
+      }
+      /* Now retriangulate upper and lower gaps. */
+      re_delaunay_triangulate(cdt, &edge->symedges[0]);
+      re_delaunay_triangulate(cdt, &edge->symedges[1]);
+    }
+    if (i < BLI_array_len(crossings) - 1) {
+      if (tnext != NULL) {
+        tstart = tnext;
 #ifdef DEBUG_CDT
         if (dbg_level > 1) {
-          fprintf(stderr, "added\n");
+          fprintf(stderr, "now tstart = ");
+          dump_se(tstart, "");
         }
 #endif
-        if (r_edges != NULL) {
-          BLI_linklist_append_pool(&edge_list, edge, cdt->listpool);
-        }
-        /* Now retriangulate upper and lower gaps. */
-        re_delaunay_triangulate(cdt, &edge->symedges[0]);
-        re_delaunay_triangulate(cdt, &edge->symedges[1]);
-      }
-      if (i < BLI_array_len(crossings) - 1) {
-        if (tnext != NULL) {
-          tstart = tnext;
-#ifdef DEBUG_CDT
-          if (dbg_level > 1) {
-            fprintf(stderr, "now tstart = ");
-            dump_se(tstart, "");
-          }
-#endif
-        }
       }
     }
   }
+
   if (r_edges) {
     *r_edges = edge_list.list;
   }
@@ -1699,6 +1936,475 @@ static void dissolve_symedge(CDT_state *cdt, SymEdge *se)
   delete_edge(cdt, se);
 }
 
+/* Slow way to get face and start vertex index within face for edge id e. */
+static bool get_face_edge_id_indices(const CDT_input *in, int e, int *r_f, int *r_fi)
+{
+  int f;
+  int id;
+
+  id = in->edges_len;
+  if (e < id) {
+    return false;
+  }
+  for (f = 0; f < in->faces_len; f++) {
+    if (e < id + in->faces_len_table[f]) {
+      *r_f = f;
+      *r_fi = e - id;
+      return true;
+    }
+    id += in->faces_len_table[f];
+  }
+  return false;
+}
+
+/* Is pt_co when snapped to segment seg1 seg2 all of:
+ * a) strictly within that segment
+ * b) within epsilon from original pt_co
+ * c) pt_co is not within epsilon of either seg1 or seg2.
+ * Return true if so, and return in *r_lambda the fraction of the way from seg1 to seg2 of the
+ * snapped point.
+ */
+static bool check_vert_near_segment(const double *pt_co,
+                                    const double *seg1,
+                                    const double *seg2,
+                                    double epsilon_squared,
+                                    double *r_lambda)
+{
+  double lambda, snap_co[2];
+
+  lambda = closest_to_line_v2_db(snap_co, pt_co, seg1, seg2);
+  *r_lambda = lambda;
+  if (lambda <= 0.0 || lambda >= 1.0) {
+    return false;
+  }
+  if (len_squared_v2v2_db(pt_co, snap_co) > epsilon_squared) {
+    return false;
+  }
+  if (len_squared_v2v2_db(pt_co, seg1) <= epsilon_squared ||
+      len_squared_v2v2_db(pt_co, seg2) <= epsilon_squared) {
+    return false;
+  }
+  return true;
+}
+
+typedef struct EdgeVertLambda {
+  int e_id;
+  int v_id;
+  double lambda;
+} EdgeVertLambda;
+
+/* For sorting first by edge id, then by lambda, then by vert id. */
+static int evl_cmp(const void *a, const void *b)
+{
+  const EdgeVertLambda *sa = a;
+  const EdgeVertLambda *sb = b;
+
+  if (sa->e_id < sb->e_id) {
+    return -1;
+  }
+  else if (sa->e_id > sb->e_id) {
+    return 1;
+  }
+  else if (sa->lambda < sb->lambda) {
+    return -1;
+  }
+  else if (sa->lambda > sb->lambda) {
+    return 1;
+  }
+  else if (sa->v_id < sb->v_id) {
+    return -1;
+  }
+  else if (sa->v_id > sb->v_id) {
+    return 1;
+  }
+  return 0;
+}
+
+/**
+ * If epsilon > 0, and input doesn't have skip_modify_input == true,
+ * check input to see if any constraint edge ends (including face edges) come
+ * within epsilon of another edge.
+ * For all such cases, we want to split the constraint edge at the point nearest to near vertex
+ * and move the vertex coordinates to be on that edge.
+ * But exclude cases where they come within epsilon of either end because those will be handled
+ * by vertex merging in the main triangulation algorithm.
+ *
+ * If any such splits are found, make a new CDT_input reflecting this change, and provide an
+ * edge map to map from edge ids in the new input space to edge ids in the old input space.
+ *
+ * TODO: replace naive O(n^2) algorithm with kdopbvh-based one.
+ */
+static const CDT_input *modify_input_for_near_edge_ends(const CDT_input *input, int **r_edge_map)
+{
+  CDT_input *new_input = NULL;
+  int e, eprev, e1, e2, f, fi, flen, start, i, j;
+  int i_new, i_old, i_evl;
+  int v11, v12, v21, v22;
+  double co11[2], co12[2], co21[2], co22[2];
+  double lambda;
+  double eps = (double)input->epsilon;
+  double eps_sq = eps * eps;
+  int tot_edge_constraints, edges_len, tot_face_edges;
+  int new_tot_face_edges, new_tot_con_edges;
+  int delta_con_edges, delta_face_edges, cur_e_cnt;
+  int *edge_map;
+  int evl_len;
+  EdgeVertLambda *edge_vert_lambda = NULL;
+  BLI_array_staticdeclare(edge_vert_lambda, 128);
+#ifdef DEBUG_CDT
+  EdgeVertLambda *evl;
+  int dbg_level = 0;
+
+  if (dbg_level > 0) {
+    fprintf(stderr, "\nMODIFY INPUT\n\n");
+  }
+#endif
+
+  *r_edge_map = NULL;
+  if (input->epsilon == 0.0 || input->skip_input_modify ||
+      (input->edges_len == 0 && input->faces_len == 0)) {
+    return input;
+  }
+
+  /* Edge constraints are union of the explicitly provided edges and the implicit edges
+   * that are part of the provided faces. We index constraints by have the first input->edges_len
+   * ints standing for the explicit edge with the same index, and the rest being face edges in
+   * the order that the faces appear and then edges within those faces, with indices offset by
+   * input->edges_len.
+   * Calculate tot_edge_constraints to be the sum of the two kinds of edges.
+   * We first have to count the number of face edges.
+   * That is the same as the number of vertices in the faces table, which
+   * we can find by adding the last length to the last start.
+   */
+  edges_len = input->edges_len;
+  tot_edge_constraints = edges_len;
+  if (input->faces_len > 0) {
+    tot_face_edges = input->faces_start_table[input->faces_len - 1] +
+                     input->faces_len_table[input->faces_len - 1];
+  }
+  else {
+    tot_face_edges = 0;
+  }
+  tot_edge_constraints = edges_len + tot_face_edges;
+
+  for (e1 = 0; e1 < tot_edge_constraints - 1; e1++) {
+    if (e1 < edges_len) {
+      v11 = input->edges[e1][0];
+      v12 = input->edges[e1][1];
+    }
+    else {
+      if (!get_face_edge_id_indices(input, e1, &f, &fi)) {
+        /* Must be bad input. Will be caught later so don't need to signal here. */
+        continue;
+      }
+      start = input->faces_start_table[f];
+      flen = input->faces_len_table[f];
+      v11 = input->faces[start + fi];
+      v12 = input->faces[(fi == flen - 1) ? start : start + fi + 1];
+    }
+    for (e2 = e1 + 1; e2 < tot_edge_constraints; e2++) {
+      if (e2 < edges_len) {
+        v21 = input->edges[e2][0];
+        v22 = input->edges[e2][1];
+      }
+      else {
+        if (!get_face_edge_id_indices(input, e2, &f, &fi)) {
+          continue;
+        }
+        start = input->faces_start_table[f];
+        flen = input->faces_len_table[f];
+        v21 = input->faces[start + fi];
+        v22 = input->faces[(fi == flen - 1) ? start : start + fi + 1];
+      }
+      copy_v2db_v2fl(co11, input->vert_coords[v11]);
+      copy_v2db_v2fl(co12, input->vert_coords[v12]);
+      copy_v2db_v2fl(co21, input->vert_coords[v21]);
+      copy_v2db_v2fl(co22, input->vert_coords[v22]);
+      if (check_vert_near_segment(co11, co21, co22, eps_sq, &lambda)) {
+
+        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e2, v11, lambda}));
+      }
+      if (check_vert_near_segment(co12, co21, co22, eps_sq, &lambda)) {
+        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e2, v12, lambda}));
+      }
+      if (check_vert_near_segment(co21, co11, co12, eps_sq, &lambda)) {
+        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e1, v21, lambda}));
+      }
+      if (check_vert_near_segment(co22, co11, co12, eps_sq, &lambda)) {
+        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e1, v22, lambda}));
+      }
+    }
+  }
+
+  evl_len = BLI_array_len(edge_vert_lambda);
+  if (evl_len > 0) {
+    /* Sort to bring splits for each edge together,
+     * and for each edge, to be in order of lambda. */
+    qsort(edge_vert_lambda, evl_len, sizeof(EdgeVertLambda), evl_cmp);
+#ifdef DEBUG_CDT
+    if (dbg_level > 0) {
+      fprintf(stderr, "\nafter sorting\n");
+      for (i = 0; i < evl_len; i++) {
+        evl = &edge_vert_lambda[i];
+        fprintf(stderr, "e%d, v%d, %g\n", evl->e_id, evl->v_id, evl->lambda);
+      }
+    }
+#endif
+
+    /* Remove dups in edge_vert_lambda, where dup means that the edge is the
+     * same, and the verts are either the same or will be merged by epsilon-nearness.
+     */
+    i = 0;
+    j = 0;
+    /* In loop, copy from position j to position i. */
+    for (j = 0; j < evl_len;) {
+      int k;
+      if (i != j) {
+        memmove(&edge_vert_lambda[i], &edge_vert_lambda[j], sizeof(EdgeVertLambda));
+      }
+      for (k = j + 1; k < evl_len; k++) {
+        int vj = edge_vert_lambda[j].v_id;
+        int vk = edge_vert_lambda[k].v_id;
+        if (vj != vk) {
+          if (len_squared_v2v2(input->vert_coords[vj], input->vert_coords[vk]) > (float)eps_sq) {
+            break;
+          }
+        }
+      }
+      j = k;
+      i++;
+    }
+
+    if (i != evl_len) {
+      evl_len = i;
+#ifdef DEBUG_CDT
+      if (dbg_level > 0) {
+        fprintf(stderr, "\nduplicates eliminated\n");
+        for (i = 0; i < evl_len; i++) {
+          evl = &edge_vert_lambda[i];
+          fprintf(stderr, "e%d, v%d, %g\n", evl->e_id, evl->v_id, evl->lambda);
+        }
+      }
+#endif
+    }
+    /* Find delta in number of constraint edges and face edges.
+     * This may be overestimates of true number, due to duplicates. */
+    delta_con_edges = 0;
+    delta_face_edges = 0;
+    cur_e_cnt = 0;
+    eprev = -1;
+    for (i = 0; i < evl_len; i++) {
+      e = edge_vert_lambda[i].e_id;
+      if (i > 0 && e > eprev) {
+        /* New edge group. Previous group had cur_e_cnt split vertices.
+         * That is the delta in the number of edges needed in input since
+         * there will be cur_e_cnt + 1 edges replacing one edge.
+         */
+        if (eprev < edges_len) {
+          delta_con_edges += cur_e_cnt;
+        }
+        else {
+          delta_face_edges += cur_e_cnt;
+        }
+        cur_e_cnt = 1;
+        ;
+      }
+      else {
+        cur_e_cnt++;
+      }
+      eprev = e;
+    }
+    if (eprev < edges_len) {
+      delta_con_edges += cur_e_cnt;
+    }
+    else {
+      delta_face_edges += cur_e_cnt;
+    }
+    new_tot_con_edges = input->edges_len + delta_con_edges;
+    if (input->faces_len > 0) {
+      new_tot_face_edges = input->faces_start_table[input->faces_len - 1] +
+                           input->faces_len_table[input->faces_len - 1] + delta_face_edges;
+    }
+    else {
+      new_tot_face_edges = 0;
+    }
+
+    /* Allocate new CDT_input, now we know sizes needed (perhaps overestimated a bit).
+     * Caller will be reponsible for freeing it and its arrays.
+     */
+    new_input = MEM_callocN(sizeof(CDT_input), __func__);
+    new_input->epsilon = input->epsilon;
+    new_input->verts_len = input->verts_len;
+    new_input->vert_coords = (float(*)[2])MEM_malloc_arrayN(
+        new_input->verts_len, 2 * sizeof(float), __func__);
+    /* We don't do it now, but may decide to change coords of snapped verts. */
+    memmove(new_input->vert_coords,
+            input->vert_coords,
+            (size_t)new_input->verts_len * sizeof(float) * 2);
+
+    if (edges_len > 0) {
+      new_input->edges_len = new_tot_con_edges;
+      new_input->edges = (int(*)[2])MEM_malloc_arrayN(
+          new_tot_con_edges, 2 * sizeof(int), __func__);
+    }
+
+    if (input->faces_len > 0) {
+      new_input->faces_len = input->faces_len;
+      new_input->faces_start_table = (int *)MEM_malloc_arrayN(
+          new_input->faces_len, sizeof(int), __func__);
+      new_input->faces_len_table = (int *)MEM_malloc_arrayN(
+          new_input->faces_len, sizeof(int), __func__);
+      new_input->faces = (int *)MEM_malloc_arrayN(new_tot_face_edges, sizeof(int), __func__);
+    }
+
+    edge_map = (int *)MEM_malloc_arrayN(
+        new_tot_con_edges + new_tot_face_edges, sizeof(int), __func__);
+    *r_edge_map = edge_map;
+
+    i_new = i_old = i_evl = 0;
+    e = edge_vert_lambda[0].e_id;
+    /* First do new constraint edges. */
+    for (i_old = 0; i_old < edges_len; i_old++) {
+      if (i_old < e) {
+        /* Edge for i_old not split; copy it into new_input. */
+        new_input->edges[i_new][0] = input->edges[i_old][0];
+        new_input->edges[i_new][1] = input->edges[i_old][1];
+        edge_map[i_new] = i_old;
+        i_new++;
+      }
+      else {
+        /* Edge for i_old is split. */
+        BLI_assert(i_old == e);
+        new_input->edges[i_new][0] = input->edges[i_old][0];
+        new_input->edges[i_new][1] = edge_vert_lambda[i_evl].v_id;
+        edge_map[i_new] = i_old;
+        i_new++;
+        i_evl++;
+        while (i_evl < evl_len && e == edge_vert_lambda[i_evl].e_id) {
+          new_input->edges[i_new][0] = new_input->edges[i_new - 1][1];
+          new_input->edges[i_new][1] = edge_vert_lambda[i_evl].v_id;
+          edge_map[i_new] = i_old;
+          i_new++;
+          i_evl++;
+        }
+        new_input->edges[i_new][0] = new_input->edges[i_new - 1][1];
+        new_input->edges[i_new][1] = input->edges[i_old][1];
+        edge_map[i_new] = i_old;
+        i_new++;
+        if (i_evl < evl_len) {
+          e = edge_vert_lambda[i_evl].e_id;
+        }
+        else {
+          e = INT_MAX;
+        }
+      }
+    }
+    BLI_assert(i_new <= new_tot_con_edges);
+    new_input->edges_len = i_new;
+
+    /* Now do face constraints. */
+    if (input->faces_len > 0) {
+      f = 0;
+      i_new = 0; /* Now will index cur place in new_input->faces. */
+      while (i_old < tot_edge_constraints) {
+        flen = input->faces_len_table[f];
+        BLI_assert(i_old - edges_len == input->faces_start_table[f]);
+        new_input->faces_start_table[f] = i_new;
+        if (i_old + flen - 1 < e) {
+          /* Face f is not split. */
+          for (j = 0; j < flen; j++) {
+            new_input->faces[i_new] = input->faces[i_old - edges_len + j];
+            edge_map[i_new + new_input->edges_len] = i_old + j;
+            i_new++;
+          }
+          i_old += flen;
+          new_input->faces_len_table[f] = flen;
+          f++;
+        }
+        else {
+          /* Face f has at least one split edge. */
+          int i_new_start = i_new;
+          for (j = 0; j < flen; j++) {
+            if (i_old + j < e) {
+              /* jth edge of f is not split. */
+              new_input->faces[i_new] = input->faces[i_old - edges_len + j];
+              edge_map[i_new + new_input->edges_len] = i_old + j;
+              i_new++;
+            }
+            else {
+              /* jth edge of f is split. */
+              BLI_assert(i_old + j == e);
+              new_input->faces[i_new] = input->faces[i_old - edges_len + j];
+              edge_map[i_new + new_input->edges_len] = i_old + j;
+              i_new++;
+              while (i_evl < evl_len && e == edge_vert_lambda[i_evl].e_id) {
+                new_input->faces[i_new] = edge_vert_lambda[i_evl].v_id;
+                edge_map[i_new + new_input->edges_len] = i_old + j;
+                i_new++;
+                i_evl++;
+              }
+              if (i_evl < evl_len) {
+                e = edge_vert_lambda[i_evl].e_id;
+              }
+              else {
+                e = INT_MAX;
+              }
+            }
+          }
+          new_input->faces_len_table[f] = i_new - i_new_start;
+          i_old += flen;
+          f++;
+        }
+      }
+    }
+
+#ifdef DEBUG_CDT
+    if (dbg_level > 0) {
+      fprintf(stderr, "\nnew constraint edges\n");
+      for (i = 0; i < new_input->edges_len; i++) {
+        fprintf(stderr, "  e%d: (%d,%d)\n", i, new_input->edges[i][0], new_input->edges[i][1]);
+      }
+      fprintf(stderr, "\nnew faces\n");
+      for (f = 0; f < new_input->faces_len; f++) {
+        flen = new_input->faces_len_table[f];
+        start = new_input->faces_start_table[f];
+        fprintf(stderr, "  f%d: start=%d, len=%d\n    ", f, start, flen);
+        for (i = start; i < start + flen; i++) {
+          fprintf(stderr, "%d ", new_input->faces[i]);
+        }
+        fprintf(stderr, "\n");
+      }
+      fprintf(stderr, "\nedge map (new->old)\n");
+      for (i = 0; i < new_tot_con_edges + new_tot_face_edges; i++) {
+        fprintf(stderr, "  %d->%d\n", i, edge_map[i]);
+      }
+    }
+#endif
+  }
+
+  BLI_array_free(edge_vert_lambda);
+  if (new_input != NULL) {
+    return (const CDT_input *)new_input;
+  }
+  else {
+    return input;
+  }
+}
+
+static void free_modified_input(CDT_input *input)
+{
+  MEM_freeN(input->vert_coords);
+  if (input->edges != NULL) {
+    MEM_freeN(input->edges);
+  }
+  if (input->faces != NULL) {
+    MEM_freeN(input->faces);
+    MEM_freeN(input->faces_len_table);
+    MEM_freeN(input->faces_start_table);
+  }
+  MEM_freeN(input);
+}
+
 /* Return true if we can merge se's vert into se->next's vert
  * without making the area of any new triangle formed by doing
  * that into a zero or negative area triangle.*/
@@ -2439,8 +3145,7 @@ static CDT_result *cdt_get_output(CDT_state *cdt,
  * face intersection.
  *
  * This routine also takes an epsilon parameter in the \a input. Input vertices closer than epsilon
- * will be merged, and we collapse tiny edges (less than epsilon length) and skinny triangles
- * (having an altitude of less than epsilon).
+ * will be merged, and we collapse tiny edges (less than epsilon length).
  *
  * The current initial Deluanay triangulation algorithm is the Guibas-Stolfi Divide and Conquer
  * algorithm (see "Primitives for the Manipulation of General Subdivisions and the Computation of
@@ -2463,13 +3168,15 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
   int nv = input->verts_len;
   int ne = input->edges_len;
   int nf = input->faces_len;
-  int i, iv1, iv2, f, fedge_start, fedge_end;
+  int i, iv1, iv2, f, fedge_start, fedge_end, ei;
   CDT_state *cdt;
   CDTVert *v1, *v2;
   CDTEdge *face_edge;
   SymEdge *face_symedge;
   LinkNode *edge_list;
   CDT_result *result;
+  const CDT_input *input_orig;
+  int *new_edge_map;
   static bool called_exactinit = false;
 #ifdef DEBUG_CDT
   int dbg_level = 0;
@@ -2504,11 +3211,26 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
     return NULL;
   }
 
-  cdt = new_cdt_init(input);
+  input_orig = input;
+  input = modify_input_for_near_edge_ends(input, &new_edge_map);
+  if (input != input_orig) {
+    nv = input->verts_len;
+    ne = input->edges_len;
+    nf = input->faces_len;
+#ifdef DEBUG_CDT
+    if (dbg_level > 4) {
+      fprintf(stderr, "input modified for near ends; now ne=%d\n", ne);
+    }
+#endif
+  }
+  cdt = cdt_init(input);
   initial_triangulation(cdt);
 #ifdef DEBUG_CDT
   if (dbg_level > 0) {
     validate_cdt(cdt, true, false, false);
+    if (dbg_level > 1) {
+      cdt_draw(cdt, "after initial triangulation");
+    }
   }
 #endif
 
@@ -2529,7 +3251,13 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
     if (v2->merge_to_index != -1) {
       v2 = cdt->vert_array[v2->merge_to_index];
     }
-    add_edge_constraint(cdt, v1, v2, i, NULL);
+    if (new_edge_map) {
+      ei = new_edge_map[i];
+    }
+    else {
+      ei = i;
+    }
+    add_edge_constraint(cdt, v1, v2, ei, NULL);
 #ifdef DEBUG_CDT
     if (dbg_level > 3) {
       char namebuf[60];
@@ -2553,6 +3281,9 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
     }
     for (i = 0; i < flen; i++) {
       int face_edge_id = cdt->face_edge_offset + fstart + i;
+      if (new_edge_map) {
+        face_edge_id = new_edge_map[face_edge_id];
+      }
       iv1 = input->faces[fstart + i];
       iv2 = input->faces[fstart + ((i + 1) % flen)];
       if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
@@ -2619,9 +3350,9 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
     remove_small_features(cdt);
 #ifdef DEBUG_CDT
     if (dbg_level > 2) {
-      cdt_draw(cdt, "after collapse skinny triangles\n");
+      cdt_draw(cdt, "after remove small features\n");
       if (dbg_level > 3) {
-        dump_cdt(cdt, "after collapse skinny triangles\n");
+        dump_cdt(cdt, "after remove small features\n");
       }
     }
 #endif
@@ -2634,6 +3365,9 @@ CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_ty
   }
 #endif
 
+  if (input != input_orig) {
+    free_modified_input((CDT_input *)input);
+  }
   new_cdt_free(cdt);
   return result;
 }
@@ -2708,14 +3442,14 @@ ATTU static const char *sename(const SymEdge *se)
 
 static void dump_v(const CDTVert *v, const char *lab)
 {
-  fprintf(stderr, "%s%s(%.3f,%.3f)\n", lab, vertname(v), F2(v->co));
+  fprintf(stderr, "%s%s(%.10f,%.10f)\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))",
+            "%s%s((%.10f,%.10f)->(%.10f,%.10f))",
             lab,
             vertname(se->vert),
             F2(se->vert->co),
@@ -2723,7 +3457,18 @@ static void dump_se(const SymEdge *se, const char *lab)
     fprintf(stderr, "%s\n", vertname(se->next->vert));
   }
   else {
-    fprintf(stderr, "%s%s((%.3f,%.3f)->NULL)\n", lab, vertname(se->vert), F2(se->vert->co));
+    fprintf(stderr, "%s%s((%.10f,%.10f)->NULL)\n", lab, vertname(se->vert), F2(se->vert->co));
+  }
+}
+
+static void dump_se_short(const SymEdge *se, const char *lab)
+{
+  if (se == NULL) {
+    fprintf(stderr, "%sNULL", lab);
+  }
+  else {
+    fprintf(stderr, "%s%s", lab, vertname(se->vert));
+    fprintf(stderr, "%s", se->next == NULL ? "[NULL]" : vertname(se->next->vert));
   }
 }
 
@@ -2754,6 +3499,20 @@ static void dump_id_list(const LinkNode *id_list, const char *lab)
   }
 }
 
+static void dump_cross_data(struct CrossData *cd, const char *lab)
+{
+  fprintf(stderr, "%s", lab);
+  if (cd->lambda == 0.0) {
+    fprintf(stderr, "v%d", cd->vert->index);
+  }
+  else {
+    fprintf(stderr, "lambda=%.17g", cd->lambda);
+  }
+  dump_se_short(cd->in, " in=");
+  dump_se_short(cd->out, " out=");
+  fprintf(stderr, "\n");
+}
+
 /* 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_filtered(const CDT_state *cdt,
@@ -2941,7 +3700,7 @@ static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const
 #  define MAX_DRAW_WIDTH 2000
 #  define MAX_DRAW_HEIGHT 1400
 #  define THIN_LINE 1
-#  define THICK_LINE 3
+#  define THICK_LINE 4
 #  define VERT_RADIUS 3
 #  define DRAW_VERT_LABELS 1
 #  define DRAW_EDGE_LABELS 0
@@ -3025,7 +3784,7 @@ static void cdt_draw_region(
             SX(v->co[0]),
             SY(v->co[1]),
             VERT_RADIUS);
-    fprintf(f, "  <title>%s(%.3f,%.3f)</title>\n", vertname(v), v->co[0], v->co[1]);
+    fprintf(f, "  <title>%s(%.10f,%.10f)</title>\n", vertname(v), v->co[0], v->co[1]);
     fprintf(f, "</circle>\n");
 #  if DRAW_VERT_LABELS
     fprintf(f,
@@ -3062,6 +3821,19 @@ static void cdt_draw_vertex_region(CDT_state *cdt, int v, double dist, const cha
   cdt_draw_region(cdt, lab, co[0] - dist, co[1] - dist, co[0] + dist, co[1] + dist);
 }
 
+static void cdt_draw_edge_region(CDT_state *cdt, int v1, int v2, double dist, const char *lab)
+{
+  const double *co1 = cdt->vert_array[v1]->co;
+  const double *co2 = cdt->vert_array[v2]->co;
+  double minx, miny, maxx, maxy;
+
+  minx = min_dd(co1[0], co2[0]);
+  miny = min_dd(co1[1], co2[1]);
+  maxx = max_dd(co1[0], co2[0]);
+  maxy = max_dd(co1[1], co2[1]);
+  cdt_draw_region(cdt, lab, minx - dist, miny - dist, maxx + dist, maxy + dist);
+}
+
 #  define CDTFILE "/tmp/cdtinput.txt"
 static void write_cdt_input_to_file(const CDT_input *inp)
 {
@@ -3087,7 +3859,7 @@ static void write_cdt_input_to_file(const CDT_input *inp)
 #  ifndef NDEBUG /* Only used in assert. */
 /*
  * Is a visible from b: i.e., ab crosses no edge of cdt?
- * If constrained is true, consider only constrained edges as possible crosser's.
+ * If constrained is true, consider only constrained edges as possible crossed edges.
  * In any case, don't count an edge ab itself.
  * Note: this is an expensive test if there are a lot of edges.
  */
@@ -3232,9 +4004,12 @@ static void validate_cdt(CDT_state *cdt,
         v1 = se->vert;
         v2 = se->next->vert;
         v3 = se->next->next->vert;
-        BLI_assert(orient2d(v1->co, v2->co, v3->co) >= 0.0);
-        BLI_assert(orient2d(v2->co, v3->co, v1->co) >= 0.0);
-        BLI_assert(orient2d(v3->co, v1->co, v2->co) >= 0.0);
+        /* The triangle should be positively oriented, but because
+         * the insertion of intersection vertices doesn't use exact
+         * arithmetic, this may not be true, so allow a little slop. */
+        BLI_assert(orient2d(v1->co, v2->co, v3->co) >= -FLT_EPSILON);
+        BLI_assert(orient2d(v2->co, v3->co, v1->co) >= -FLT_EPSILON);
+        BLI_assert(orient2d(v3->co, v1->co, v2->co) >= -FLT_EPSILON);
       }
       UNUSED_VARS_NDEBUG(limit);
       BLI_assert(se->next->next != se);
@@ -3466,7 +4241,7 @@ static double isperrboundA, isperrboundB, isperrboundC;
  *  error.  It is used for floating-point error analysis.
  *
  *  `splitter' is used to split floating-point numbers into two
- *  half-length significances for exact multiplication.
+ *  half-length significands for exact multiplication.
  *
  *  I imagine that a highly optimizing compiler might be too smart for its
  *  own good, and somehow cause this routine to fail, if it pretends that
diff --git a/source/blender/python/mathutils/mathutils_geometry.c b/source/blender/python/mathutils/mathutils_geometry.c
index 56bb60bf921..6474275ecaf 100644
--- a/source/blender/python/mathutils/mathutils_geometry.c
+++ b/source/blender/python/mathutils/mathutils_geometry.c
@@ -1649,6 +1649,7 @@ static PyObject *M_Geometry_delaunay_2d_cdt(PyObject *UNUSED(self), PyObject *ar
   in.faces_start_table = in_faces_start_table;
   in.faces_len_table = in_faces_len_table;
   in.epsilon = epsilon;
+  in.skip_input_modify = false;
 
   res = BLI_delaunay_2d_cdt_calc(&in, output_type);
 
diff --git a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc
index f377d5a8247..30be1a6df29 100644
--- a/tests/gtests/blenlib/BLI_delaunay_2d_test.cc
+++ b/tests/gtests/blenlib/BLI_delaunay_2d_test.cc
@@ -31,6 +31,7 @@ static void fill_input_verts(CDT_input *r_input, float (*vcos)[2], int nverts)
   r_input->faces_start_table = NULL;
   r_input->faces_len_table = NULL;
   r_input->epsilon = 1e-5f;
+  r_input->skip_input_modify = false;
 }
 
 static void add_input_edges(CDT_input *r_input, int (*edges)[2], int nedges)
@@ -1099,6 +1100,246 @@ TEST(delaunay, repeatededge)
   free_spec_arrays(&in);
   BLI_delaunay_2d_cdt_free(out);
 }
+
+TEST(delaunay, NearSeg)
+{
+  CDT_input in;
+  CDT_result *out;
+  int v[4], e0, e1, e2, i;
+  const char *spec = R"(4 2 0
+  0.0 0.0
+  1.0 0.0
+  0.25 0.09
+  0.25 1.0
+  0 1
+  2 3
+  )";
+
+  fill_input_from_string(&in, spec);
+  in.epsilon = 0.1;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 4);
+  EXPECT_EQ(out->edges_len, 3);
+  EXPECT_EQ(out->faces_len, 0);
+  if (out->edges_len == 3) {
+    for (i = 0; i < 4; i++) {
+      v[i] = get_output_vert_index(out, i);
+      EXPECT_NE(v[i], -1);
+    }
+    e0 = get_edge(out, v[0], v[2]);
+    e1 = get_edge(out, v[2], v[1]);
+    e2 = get_edge(out, v[2], v[3]);
+    EXPECT_TRUE(out_edge_has_input_id(out, e0, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e2, 1));
+  }
+  free_spec_arrays(&in);
+  BLI_delaunay_2d_cdt_free(out);
+}
+
+TEST(delaunay, OverlapSegs)
+{
+  CDT_input in;
+  CDT_result *out;
+  int v[4], e0, e1, e2, i;
+  const char *spec = R"(4 2 0
+  0.0 0.0
+  1.0 0.0
+  0.4 0.09
+  1.4 0.09
+  0 1
+  2 3
+  )";
+
+  fill_input_from_string(&in, spec);
+  in.epsilon = 0.1;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 4);
+  EXPECT_EQ(out->edges_len, 3);
+  EXPECT_EQ(out->faces_len, 0);
+  if (out->edges_len == 3) {
+    for (i = 0; i < 4; i++) {
+      v[i] = get_output_vert_index(out, i);
+      EXPECT_NE(v[i], -1);
+    }
+    e0 = get_edge(out, v[0], v[2]);
+    e1 = get_edge(out, v[2], v[1]);
+    e2 = get_edge(out, v[1], v[3]);
+    EXPECT_TRUE(out_edge_has_input_id(out, e0, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 1));
+    EXPECT_TRUE(out_edge_has_input_id(out, e2, 1));
+  }
+  free_spec_arrays(&in);
+  BLI_delaunay_2d_cdt_free(out);
+}
+
+TEST(delaunay, NearSegWithDup)
+{
+  CDT_input in;
+  CDT_result *out;
+  int v[5], e0, e1, e2, e3, i;
+  const char *spec = R"(5 3 0
+  0.0 0.0
+  1.0 0.0
+  0.25 0.09
+  0.25 1.0
+  0.75 0.09
+  0 1
+  2 3
+  2 4
+  )";
+
+  fill_input_from_string(&in, spec);
+  in.epsilon = 0.1;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 5);
+  EXPECT_EQ(out->edges_len, 4);
+  EXPECT_EQ(out->faces_len, 0);
+  if (out->edges_len == 5) {
+    for (i = 0; i < 5; i++) {
+      v[i] = get_output_vert_index(out, i);
+      EXPECT_NE(v[i], -1);
+    }
+    e0 = get_edge(out, v[0], v[2]);
+    e1 = get_edge(out, v[2], v[4]);
+    e2 = get_edge(out, v[4], v[1]);
+    e3 = get_edge(out, v[3], v[2]);
+    EXPECT_TRUE(out_edge_has_input_id(out, e0, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 2));
+    EXPECT_TRUE(out_edge_has_input_id(out, e2, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e3, 1));
+  }
+  free_spec_arrays(&in);
+  BLI_delaunay_2d_cdt_free(out);
+}
+
+TEST(delaunay, TwoNearSeg)
+{
+  CDT_input in;
+  CDT_result *out;
+  int v[5], e0, e1, e2, e3, e4, i;
+  const char *spec = R"(5 3 0
+  0.0 0.0
+  1.0 0.0
+  0.25 0.09
+  0.25 1.0
+  0.75 0.09
+  0 1
+  3 2
+  3 4
+  )";
+
+  fill_input_from_string(&in, spec);
+  in.epsilon = 0.1;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 5);
+  EXPECT_EQ(out->edges_len, 5);
+  EXPECT_EQ(out->faces_len, 1);
+  if (out->edges_len == 5) {
+    for (i = 0; i < 5; i++) {
+      v[i] = get_output_vert_index(out, i);
+      EXPECT_NE(v[i], -1);
+    }
+    e0 = get_edge(out, v[0], v[2]);
+    e1 = get_edge(out, v[2], v[4]);
+    e2 = get_edge(out, v[4], v[1]);
+    e3 = get_edge(out, v[3], v[2]);
+    e4 = get_edge(out, v[3], v[4]);
+    EXPECT_TRUE(out_edge_has_input_id(out, e0, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e2, 0));
+    EXPECT_TRUE(out_edge_has_input_id(out, e3, 1));
+    EXPECT_TRUE(out_edge_has_input_id(out, e4, 2));
+  }
+  free_spec_arrays(&in);
+  BLI_delaunay_2d_cdt_free(out);
+}
+
+TEST(delaunay, FaceNearSegs)
+{
+  CDT_input in;
+  CDT_result *out;
+  int v[9], e0, e1, e2, e3, i;
+  const char *spec = R"(8 1 2
+  0.0 0.0
+  2.0 0.0
+  1.0 1.0
+  0.21 0.2
+  1.79 0.2
+  0.51 0.5
+  1.49 0.5
+  1.0 0.19
+  2 7
+  0 1 2
+  3 4 6 5
+  )";
+
+  fill_input_from_string(&in, spec);
+  in.epsilon = 0.05;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 9);
+  EXPECT_EQ(out->edges_len, 13);
+  EXPECT_EQ(out->faces_len, 5);
+  if (out->verts_len == 9 && out->edges_len == 13) {
+    for (i = 0; i < 9; i++) {
+      v[i] = get_output_vert_index(out, i);
+      EXPECT_NE(v[i], -1);
+    }
+    e0 = get_edge(out, v[0], v[1]);
+    e1 = get_edge(out, v[4], v[6]);
+    e2 = get_edge(out, v[3], v[0]);
+    e3 = get_edge(out, v[2], v[8]);
+
+    EXPECT_TRUE(out_edge_has_input_id(out, e0, 1));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 2));
+    EXPECT_TRUE(out_edge_has_input_id(out, e1, 5));
+    EXPECT_TRUE(out_edge_has_input_id(out, e2, 3));
+    EXPECT_TRUE(out_edge_has_input_id(out, e3, 0));
+  }
+  free_spec_arrays(&in);
+  BLI_delaunay_2d_cdt_free(out);
+}
+
+TEST(delaunay, ChainNearIntersects)
+{
+  CDT_input in;
+  CDT_result *out;
+  const char *spec = R"(6 10 0
+  0.8 1.25
+  1.25 0.75
+  3.25 1.25
+  5.0 1.9
+  2.5 4.0
+  1.0 2.25
+  0 1
+  1 2
+  2 3
+  3 4
+  4 5
+  5 0
+  0 2
+  5 2
+  4 2
+  1 3
+  )";
+
+  fill_input_from_string(&in, spec);
+  in.epsilon = 0.05;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 9);
+  EXPECT_EQ(out->edges_len, 16);
+  BLI_delaunay_2d_cdt_free(out);
+  in.epsilon = 0.11;
+  /* The chaining we want to test happens prematurely if modify input. */
+  in.skip_input_modify = true;
+  out = BLI_delaunay_2d_cdt_calc(&in, CDT_CONSTRAINTS);
+  EXPECT_EQ(out->verts_len, 6);
+  EXPECT_EQ(out->edges_len, 9);
+  free_spec_arrays(&in);
+  BLI_delaunay_2d_cdt_free(out);
+}
 #endif
 
 #if DO_RANDOM_TESTS
@@ -1112,8 +1353,12 @@ enum {
 };
 
 #  define DO_TIMING
-static void rand_delaunay_test(
-    int test_kind, int max_lg_size, int reps_per_size, double param, CDT_output_type otype)
+static void rand_delaunay_test(int test_kind,
+                               int start_lg_size,
+                               int max_lg_size,
+                               int reps_per_size,
+                               double param,
+                               CDT_output_type otype)
 {
   CDT_input in;
   CDT_result *out;
@@ -1182,12 +1427,6 @@ static void rand_delaunay_test(
       fprintf(stderr, "unknown random delaunay test kind\n");
       return;
   }
-  fprintf(stderr,
-          "size_max=%d, npts_max=%d, nedges_max=%d, nfaces_max=%d\n",
-          size_max,
-          npts_max,
-          nedges_max,
-          nfaces_max); /*DEBUG!!*/
   p = (float(*)[2])MEM_malloc_arrayN(npts_max, 2 * sizeof(float), __func__);
   if (nedges_max > 0) {
     e = (int(*)[2])MEM_malloc_arrayN(nedges_max, 2 * sizeof(int), __func__);
@@ -1201,7 +1440,7 @@ static void rand_delaunay_test(
   times = (double *)MEM_malloc_arrayN(max_lg_size + 1, sizeof(double), __func__);
 
   /* For powers of 2 sizes up to max_lg_size power of 2. */
-  for (lg_size = 0; lg_size <= max_lg_size; lg_size++) {
+  for (lg_size = start_lg_size; lg_size <= max_lg_size; lg_size++) {
     size = 1 << lg_size;
     nedges = 0;
     nfaces = 0;
@@ -1344,62 +1583,62 @@ static void rand_delaunay_test(
 
 TEST(delaunay, randompts)
 {
-  rand_delaunay_test(RANDOM_PTS, 7, 1, 0.0, CDT_FULL);
+  rand_delaunay_test(RANDOM_PTS, 0, 7, 1, 0.0, CDT_FULL);
 }
 
 TEST(delaunay, randomsegs)
 {
-  rand_delaunay_test(RANDOM_SEGS, 7, 1, 0.0, CDT_FULL);
+  rand_delaunay_test(RANDOM_SEGS, 1, 7, 1, 0.0, CDT_FULL);
 }
 
 TEST(delaunay, randompoly)
 {
-  rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_FULL);
+  rand_delaunay_test(RANDOM_POLY, 1, 7, 1, 0.0, CDT_FULL);
 }
 
 TEST(delaunay, randompoly_inside)
 {
-  rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_INSIDE);
+  rand_delaunay_test(RANDOM_POLY, 1, 7, 1, 0.0, CDT_INSIDE);
 }
 
 TEST(delaunay, randompoly_constraints)
 {
-  rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_CONSTRAINTS);
+  rand_delaunay_test(RANDOM_POLY, 1, 7, 1, 0.0, CDT_CONSTRAINTS);
 }
 
 TEST(delaunay, randompoly_validbmesh)
 {
-  rand_delaunay_test(RANDOM_POLY, 7, 1, 0.0, CDT_CONSTRAINTS_VALID_BMESH);
+  rand_delaunay_test(RANDOM_POLY, 1, 7, 1, 0.0, CDT_CONSTRAINTS_VALID_BMESH);
 }
 
 TEST(delaunay, grid)
 {
-  rand_delaunay_test(RANDOM_TILTED_GRID, 6, 1, 0.0, CDT_FULL);
+  rand_delaunay_test(RANDOM_TILTED_GRID, 1, 6, 1, 0.0, CDT_FULL);
 }
 
 TEST(delaunay, tilted_grid_a)
 {
-  rand_delaunay_test(RANDOM_TILTED_GRID, 6, 1, 1.0, CDT_FULL);
+  rand_delaunay_test(RANDOM_TILTED_GRID, 1, 6, 1, 1.0, CDT_FULL);
 }
 
 TEST(delaunay, tilted_grid_b)
 {
-  rand_delaunay_test(RANDOM_TILTED_GRID, 6, 1, 0.01, CDT_FULL);
+  rand_delaunay_test(RANDOM_TILTED_GRID, 1, 6, 1, 0.01, CDT_FULL);
 }
 
 TEST(delaunay, randomcircle)
 {
-  rand_delaunay_test(RANDOM_CIRCLE, 7, 1, 0.0, CDT_FULL);
+  rand_delaunay_test(RANDOM_CIRCLE, 1, 7, 1, 0.0, CDT_FULL);
 }
 
 TEST(delaunay, random_tris_circle)
 {
-  rand_delaunay_test(RANDOM_TRI_BETWEEN_CIRCLES, 6, 1, 0.25, CDT_FULL);
+  rand_delaunay_test(RANDOM_TRI_BETWEEN_CIRCLES, 1, 6, 1, 0.25, CDT_FULL);
 }
 
 TEST(delaunay, random_tris_circle_b)
 {
-  rand_delaunay_test(RANDOM_TRI_BETWEEN_CIRCLES, 6, 1, 1e-4, CDT_FULL);
+  rand_delaunay_test(RANDOM_TRI_BETWEEN_CIRCLES, 1, 6, 1, 1e-4, CDT_FULL);
 }
 #endif