Merge newboolean branch into master.

This is for design task T67744, Boolean Redesign.
It adds a choice of solver to the Boolean modifier and the
Intersect (Boolean) and Intersect (Knife) tools.
The 'Fast' choice is the current Bmesh boolean.
The new 'Exact' choice is a more advanced algorithm that supports
overlapping geometry and uses more robust calculations, but is
slower than the Fast choice.
The default with this commit is set to 'Exact'. We can decide before
the 2.91 release whether or not this is the right choice, but this
choice now will get us more testing and feedback on the new code.

1 files changed, 2500 insertions, 0 deletions

diff --git a/source/blender/blenlib/intern/delaunay_2d.cc b/source/blender/blenlib/intern/delaunay_2d.cc
new file mode 100644
index 00000000000..7b0f6a658ce
--- /dev/null
+++ b/source/blender/blenlib/intern/delaunay_2d.cc
@@ -0,0 +1,2500 @@
+/*
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+/** \file
+ * \ingroup bli
+ */
+
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+
+#include "BLI_array.hh"
+#include "BLI_double2.hh"
+#include "BLI_linklist.h"
+#include "BLI_math_boolean.hh"
+#include "BLI_math_mpq.hh"
+#include "BLI_mpq2.hh"
+#include "BLI_vector.hh"
+
+#include "BLI_delaunay_2d.h"
+
+namespace blender::meshintersect {
+
+/* Throughout this file, template argument T will be an
+ * arithmetic-like type, like float, double, or mpq_class. */
+
+template<typename T> T math_abs(const T v)
+{
+  return (v < 0) ? -v : v;
+}
+
+#ifdef WITH_GMP
+template<> mpq_class math_abs<mpq_class>(const mpq_class v)
+{
+  return abs(v);
+}
+#endif
+
+template<> double math_abs<double>(const double v)
+{
+  return fabs(v);
+}
+
+template<typename T> double math_to_double(const T UNUSED(v))
+{
+  BLI_assert(false); /* Need implementation for other type. */
+  return 0.0;
+}
+
+#ifdef WITH_GMP
+template<> double math_to_double<mpq_class>(const mpq_class v)
+{
+  return v.get_d();
+}
+#endif
+
+template<> double math_to_double<double>(const double v)
+{
+  return v;
+}
+
+/**
+ * Define a templated 2D arrangement of vertices, edges, and faces.
+ * The #SymEdge data structure is the basis for a structure that allows
+ * easy traversal to neighboring (by topology) geometric elements.
+ * Each of #CDTVert, #CDTEdge, and #CDTFace have an input_id linked list,
+ * whose nodes contain integers that keep track of which input verts, edges,
+ * and faces, respectively, that the element was derived from.
+ *
+ * While this could be cleaned up some, it is usable by other routines in Blender
+ * that need to keep track of a 2D arrangement, with topology.
+ */
+template<typename Arith_t> struct CDTVert;
+template<typename Arith_t> struct CDTEdge;
+template<typename Arith_t> struct CDTFace;
+
+template<typename Arith_t> struct SymEdge {
+  /** Next #SymEdge in face, doing CCW traversal of face. */
+  SymEdge<Arith_t> *next{nullptr};
+  /** Next #SymEdge CCW around vert. */
+  SymEdge<Arith_t> *rot{nullptr};
+  /** Vert at origin. */
+  CDTVert<Arith_t> *vert{nullptr};
+  /** Un-directed edge this is for. */
+  CDTEdge<Arith_t> *edge{nullptr};
+  /** Face on left side. */
+  CDTFace<Arith_t> *face{nullptr};
+
+  SymEdge() = default;
+};
+
+/**
+ * Return other #SymEdge for same #CDTEdge as \a se.
+ */
+template<typename T> inline SymEdge<T> *sym(const SymEdge<T> *se)
+{
+  return se->next->rot;
+}
+
+/** Return #SymEdge whose next is \a se. */
+template<typename T> inline SymEdge<T> *prev(const SymEdge<T> *se)
+{
+  return se->rot->next->rot;
+}
+
+template<typename Arith_t> struct CDTVert {
+  /** Coordinate. */
+  vec2<Arith_t> co;
+  /** Some edge attached to it. */
+  SymEdge<Arith_t> *symedge{nullptr};
+  /** List of corresponding vertex input ids. */
+  LinkNode *input_ids{nullptr};
+  /** Index into array that #CDTArrangement keeps. */
+  int index{-1};
+  /** Index of a CDTVert that this has merged to. -1 if no merge. */
+  int merge_to_index{-1};
+  /** Used by algorithms operating on CDT structures. */
+  int visit_index{0};
+
+  CDTVert() = default;
+  explicit CDTVert(const vec2<Arith_t> &pt);
+};
+
+template<typename Arith_t> struct CDTEdge {
+  /** List of input edge ids that this is part of. */
+  LinkNode *input_ids{nullptr};
+  /** The directed edges for this edge. */
+  SymEdge<Arith_t> symedges[2]{SymEdge<Arith_t>(), SymEdge<Arith_t>()};
+
+  CDTEdge() = default;
+};
+
+template<typename Arith_t> struct CDTFace {
+  /** A symedge in face; only used during output, so only valid then. */
+  SymEdge<Arith_t> *symedge{nullptr};
+  /** List of input face ids that this is part of. */
+  LinkNode *input_ids{nullptr};
+  /** Used by algorithms operating on CDT structures. */
+  int visit_index{0};
+  /** Marks this face no longer used. */
+  bool deleted{false};
+
+  CDTFace() = default;
+};
+
+template<typename Arith_t> struct CDTArrangement {
+  /* The arrangement owns the memory pointed to by the pointers in these vectors.
+   * They are pointers instead of actual structures because these vectors may be resized and
+   * other elements refer to the elements by pointer. */
+
+  /** The verts. Some may be merged to others (see their merge_to_index). */
+  Vector<CDTVert<Arith_t> *> verts;
+  /** The edges. Some may be deleted (SymEdge next and rot pointers are null). */
+  Vector<CDTEdge<Arith_t> *> edges;
+  /** The faces. Some may be deleted (see their delete member). */
+  Vector<CDTFace<Arith_t> *> faces;
+  /** Which CDTFace is the outer face. */
+  CDTFace<Arith_t> *outer_face{nullptr};
+
+  CDTArrangement() = default;
+  ~CDTArrangement();
+
+  /** Hint to how much space to reserve in the Vectors of the arrangement,
+   * based on these counts of input elements. */
+  void reserve(int num_verts, int num_edges, int num_faces);
+
+  /**
+   * Add a new vertex to the arrangement, with the given 2D coordinate.
+   * It will not be connected to anything yet.
+   */
+  CDTVert<Arith_t> *add_vert(const vec2<Arith_t> &pt);
+
+  /**
+   * Add an edge from v1 to v2. The edge will have a left face and a right face,
+   * specified by \a fleft and \a fright. The edge will not be connected to anything yet.
+   * If the vertices do not yet have a #SymEdge pointer,
+   * their pointer is set to the #SymEdge in this new edge.
+   */
+  CDTEdge<Arith_t> *add_edge(CDTVert<Arith_t> *v1,
+                             CDTVert<Arith_t> *v2,
+                             CDTFace<Arith_t> *fleft,
+                             CDTFace<Arith_t> *fright);
+
+  /**
+   * Add a new face. It is disconnected until an add_edge makes it the
+   * left or right face of an edge.
+   */
+  CDTFace<Arith_t> *add_face();
+
+  /** Make a new edge from v to se->vert, splicing it in. */
+  CDTEdge<Arith_t> *add_vert_to_symedge_edge(CDTVert<Arith_t> *v, SymEdge<Arith_t> *se);
+
+  /**
+   * Assuming s1 and s2 are both #SymEdge's in a face with > 3 sides and one is not the next of the
+   * other, Add an edge from `s1->v` to `s2->v`, splitting the face in two. The original face will
+   * be the one that s1 has as left face, and a new face will be added and made s2 and its
+   * next-cycle's left face.
+   */
+  CDTEdge<Arith_t> *add_diagonal(SymEdge<Arith_t> *s1, SymEdge<Arith_t> *s2);
+
+  /**
+   * Connect the verts of se1 and se2, assuming that currently those two #SymEdge's are on the
+   * outer boundary (have face == outer_face) of two components that are isolated from each other.
+   */
+  CDTEdge<Arith_t> *connect_separate_parts(SymEdge<Arith_t> *se1, SymEdge<Arith_t> *se2);
+
+  /**
+   * Split se at fraction lambda, and return the new #CDTEdge that is the new second half.
+   * Copy the edge input_ids into the new one.
+   */
+  CDTEdge<Arith_t> *split_edge(SymEdge<Arith_t> *se, Arith_t lambda);
+
+  /**
+   * Delete an edge. The new combined face on either side of the deleted edge will be the one that
+   * was e's face. There will now be an unused face, which will be marked deleted, and an unused
+   * #CDTEdge, marked by setting the next and rot pointers of its #SymEdge's to #nullptr.
+   */
+  void delete_edge(SymEdge<Arith_t> *se);
+
+  /**
+   * If the vertex with index i in the vert array has not been merge, return it.
+   * Else return the one that it has merged to.
+   */
+  CDTVert<Arith_t> *get_vert_resolve_merge(int i)
+  {
+    CDTVert<Arith_t> *v = this->verts[i];
+    if (v->merge_to_index != -1) {
+      v = this->verts[v->merge_to_index];
+    }
+    return v;
+  }
+};
+
+template<typename T> class CDT_state {
+ public:
+  CDTArrangement<T> cdt;
+  /** How many verts were in input (will be first in vert_array). */
+  int input_vert_tot;
+  /** Used for visiting things without having to initialized their visit fields. */
+  int visit_count;
+  /**
+   * Edge ids for face start with this, and each face gets this much index space
+   * to encode positions within the face.
+   */
+  int face_edge_offset;
+  /** How close before coords considered equal. */
+  T epsilon;
+
+  explicit CDT_state(int num_input_verts, int num_input_edges, int num_input_faces, T epsilon);
+  ~CDT_state()
+  {
+  }
+};
+
+template<typename T> CDTArrangement<T>::~CDTArrangement()
+{
+  for (int i : this->verts.index_range()) {
+    CDTVert<T> *v = this->verts[i];
+    BLI_linklist_free(v->input_ids, nullptr);
+    delete v;
+    this->verts[i] = nullptr;
+  }
+  for (int i : this->edges.index_range()) {
+    CDTEdge<T> *e = this->edges[i];
+    BLI_linklist_free(e->input_ids, nullptr);
+    delete e;
+    this->edges[i] = nullptr;
+  }
+  for (int i : this->faces.index_range()) {
+    CDTFace<T> *f = this->faces[i];
+    BLI_linklist_free(f->input_ids, nullptr);
+    delete f;
+    this->faces[i] = nullptr;
+  }
+}
+
+#define DEBUG_CDT
+#ifdef DEBUG_CDT
+/* Some functions to aid in debugging. */
+template<typename T> std::string vertname(const CDTVert<T> *v)
+{
+  std::stringstream ss;
+  ss << "[" << v->index << "]";
+  return ss.str();
+}
+
+/* Abbreviated pointer value is easier to read in dumps. */
+static std::string trunc_ptr(const void *p)
+{
+  constexpr int TRUNC_PTR_MASK = 0xFFFF;
+  std::stringstream ss;
+  ss << std::hex << (POINTER_AS_INT(p) & TRUNC_PTR_MASK);
+  return ss.str();
+}
+
+template<typename T> std::string sename(const SymEdge<T> *se)
+{
+  std::stringstream ss;
+  ss << "{" << trunc_ptr(se) << "}";
+  return ss.str();
+}
+
+template<typename T> std::ostream &operator<<(std::ostream &os, const SymEdge<T> &se)
+{
+  if (se.next) {
+    os << vertname(se.vert) << "(" << se.vert->co << "->" << se.next->vert->co << ")"
+       << vertname(se.next->vert);
+  }
+  else {
+    os << vertname(se.vert) << "(" << se.vert->co << "->NULL)";
+  }
+  return os;
+}
+
+template<typename T> std::ostream &operator<<(std::ostream &os, const SymEdge<T> *se)
+{
+  os << *se;
+  return os;
+}
+
+template<typename T> std::string short_se_dump(const SymEdge<T> *se)
+{
+  if (se == nullptr) {
+    return std::string("NULL");
+  }
+  return vertname(se->vert) +
+         (se->next == nullptr ? std::string("[NULL]") : vertname(se->next->vert));
+}
+
+template<typename T> std::ostream &operator<<(std::ostream &os, const CDT_state<T> &cdt_state)
+{
+  os << "\nCDT\n\nVERTS\n";
+  for (const CDTVert<T> *v : cdt_state.cdt.verts) {
+    os << vertname(v) << " " << trunc_ptr(v) << ": " << v->co
+       << " symedge=" << trunc_ptr(v->symedge);
+    if (v->merge_to_index == -1) {
+      os << "\n";
+    }
+    else {
+      os << " merge to " << vertname(cdt_state.cdt.verts[v->merge_to_index]) << "\n";
+    }
+    const SymEdge<T> *se = v->symedge;
+    int cnt = 0;
+    constexpr int print_count_limit = 25;
+    if (se) {
+      os << "  edges out:\n";
+      do {
+        if (se->next == NULL) {
+          os << "    [NULL] next/rot symedge, se=" << trunc_ptr(se) << "\n";
+          break;
+        }
+        if (se->next->next == NULL) {
+          os << "    [NULL] next-next/rot symedge, se=" << trunc_ptr(se) << "\n";
+          break;
+        }
+        const CDTVert<T> *vother = sym(se)->vert;
+        os << "    " << vertname(vother) << "(e=" << trunc_ptr(se->edge)
+           << ", se=" << trunc_ptr(se) << ")\n";
+        se = se->rot;
+        cnt++;
+      } while (se != v->symedge && cnt < print_count_limit);
+      os << "\n";
+    }
+  }
+  os << "\nEDGES\n";
+  for (const CDTEdge<T> *e : cdt_state.cdt.edges) {
+    if (e->symedges[0].next == nullptr) {
+      continue;
+    }
+    os << trunc_ptr(&e) << ":\n";
+    for (int i = 0; i < 2; ++i) {
+      const SymEdge<T> *se = &e->symedges[i];
+      os << "  se[" << i << "] @" << trunc_ptr(se) << " next=" << trunc_ptr(se->next)
+         << ", rot=" << trunc_ptr(se->rot) << ", vert=" << trunc_ptr(se->vert) << " "
+         << vertname(se->vert) << " " << se->vert->co << ", edge=" << trunc_ptr(se->edge)
+         << ", face=" << trunc_ptr(se->face) << "\n";
+    }
+  }
+  os << "\nFACES\n";
+  os << "outer_face=" << trunc_ptr(cdt_state.cdt.outer_face) << "\n";
+  /* Only after prepare_output do faces have non-null symedges. */
+  if (cdt_state.cdt.outer_face->symedge != nullptr) {
+    for (const CDTFace<T> *f : cdt_state.cdt.faces) {
+      if (!f->deleted) {
+        os << trunc_ptr(f) << " symedge=" << trunc_ptr(f->symedge) << "\n";
+      }
+    }
+  }
+  return os;
+}
+
+template<typename T> void cdt_draw(const std::string &label, const CDTArrangement<T> &cdt)
+{
+  static bool append = false; /* Will be set to true after first call. */
+
+/* Would like to use #BKE_tempdir_base() here, but that brings in dependence on kernel library.
+ * This is just for developer debugging anyway, and should never be called in production Blender.
+ */
+#  ifdef _WIN32
+  const char *drawfile = "./debug_draw.html";
+#  else
+  const char *drawfile = "/tmp/debug_draw.html";
+#  endif
+  constexpr int max_draw_width = 1800;
+  constexpr int max_draw_height = 1600;
+  constexpr double margin_expand = 0.05;
+  constexpr int thin_line = 1;
+  constexpr int thick_line = 4;
+  constexpr int vert_radius = 3;
+  constexpr bool draw_vert_labels = true;
+  constexpr bool draw_edge_labels = false;
+
+  if (cdt.verts.size() == 0) {
+    return;
+  }
+  vec2<double> vmin(DBL_MAX, DBL_MAX);
+  vec2<double> vmax(-DBL_MAX, -DBL_MAX);
+  for (const CDTVert<T> *v : cdt.verts) {
+    for (int i = 0; i < 2; ++i) {
+      double dvi = math_to_double(v->co[i]);
+      if (dvi < vmin[i]) {
+        vmin[i] = dvi;
+      }
+      if (dvi > vmax[i]) {
+        vmax[i] = dvi;
+      }
+    }
+  }
+  double draw_margin = ((vmax.x - vmin.x) + (vmax.y - vmin.y)) * margin_expand;
+  double minx = vmin.x - draw_margin;
+  double maxx = vmax.x + draw_margin;
+  double miny = vmin.y - draw_margin;
+  double maxy = vmax.y + draw_margin;
+
+  double width = maxx - minx;
+  double height = maxy - miny;
+  double aspect = height / width;
+  int view_width = max_draw_width;
+  int view_height = static_cast<int>(view_width * aspect);
+  if (view_height > max_draw_height) {
+    view_height = max_draw_height;
+    view_width = static_cast<int>(view_height / aspect);
+  }
+  double scale = view_width / width;
+
+#  define SX(x) ((math_to_double(x) - minx) * scale)
+#  define SY(y) ((maxy - math_to_double(y)) * scale)
+
+  std::ofstream f;
+  if (append) {
+    f.open(drawfile, std::ios_base::app);
+  }
+  else {
+    f.open(drawfile);
+  }
+  if (!f) {
+    std::cout << "Could not open file " << drawfile << "\n";
+    return;
+  }
+
+  f << "<div>" << label << "</div>\n<div>\n"
+    << "<svg version=\"1.1\" "
+       "xmlns=\"http://www.w3.org/2000/svg\" "
+       "xmlns:xlink=\"http://www.w3.org/1999/xlink\" "
+       "xml:space=\"preserve\"\n"
+    << "width=\"" << view_width << "\" height=\"" << view_height << "\">n";
+
+  for (const CDTEdge<T> *e : cdt.edges) {
+    if (e->symedges[0].next == nullptr) {
+      continue;
+    }
+    const CDTVert<T> *u = e->symedges[0].vert;
+    const CDTVert<T> *v = e->symedges[1].vert;
+    const vec2<T> &uco = u->co;
+    const vec2<T> &vco = v->co;
+    int strokew = e->input_ids == nullptr ? thin_line : thick_line;
+    f << "<line fill=\"none\" stroke=\"black\" stroke-width=\"" << strokew << "\" x1=\""
+      << SX(uco[0]) << "\" y1=\"" << SY(uco[1]) << "\" x2=\"" << SX(vco[0]) << "\" y2=\""
+      << SY(vco[1]) << "\">\n";
+    f << "  <title>" << vertname(u) << vertname(v) << "</title>\n";
+    f << "</line>\n";
+    if (draw_edge_labels) {
+      f << "<text x=\"" << SX((uco[0] + vco[0]) / 2) << "\" y=\"" << SY((uco[1] + vco[1]) / 2)
+        << "\" font-size=\"small\">";
+      f << vertname(u) << vertname(v) << sename(&e->symedges[0]) << sename(&e->symedges[1])
+        << "</text>\n";
+    }
+  }
+
+  int i = 0;
+  for (const CDTVert<T> *v : cdt.verts) {
+    f << "<circle fill=\"black\" cx=\"" << SX(v->co[0]) << "\" cy=\"" << SY(v->co[1]) << "\" r=\""
+      << vert_radius << "\">\n";
+    f << "  <title>[" << i << "]" << v->co << "</title>\n";
+    f << "</circle>\n";
+    if (draw_vert_labels) {
+      f << "<text x=\"" << SX(v->co[0]) + vert_radius << "\" y=\"" << SY(v->co[1]) - vert_radius
+        << "\" font-size=\"small\">[" << i << "]</text>\n";
+    }
+    ++i;
+  }
+
+  append = true;
+#  undef SX
+#  undef SY
+}
+#endif
+
+/**
+ * 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.
+ */
+template<typename T> bool in_line(const vec2<T> &a, const vec2<T> &b, const vec2<T> &c)
+{
+  vec2<T> ab = b - a;
+  vec2<T> bc = c - b;
+  vec2<T> ac = c - a;
+  if (vec2<T>::dot(ab, ac) < 0) {
+    return false;
+  }
+  return vec2<T>::dot(bc, ac) >= 0;
+}
+
+template<typename T> CDTVert<T>::CDTVert(const vec2<T> &pt)
+{
+  this->co = pt;
+  this->input_ids = nullptr;
+  this->symedge = nullptr;
+  this->index = -1;
+  this->merge_to_index = -1;
+  this->visit_index = 0;
+}
+
+template<typename T> CDTVert<T> *CDTArrangement<T>::add_vert(const vec2<T> &pt)
+{
+  CDTVert<T> *v = new CDTVert<T>(pt);
+  int index = this->verts.append_and_get_index(v);
+  v->index = index;
+  return v;
+}
+
+template<typename T>
+CDTEdge<T> *CDTArrangement<T>::add_edge(CDTVert<T> *v1,
+                                        CDTVert<T> *v2,
+                                        CDTFace<T> *fleft,
+                                        CDTFace<T> *fright)
+{
+  CDTEdge<T> *e = new CDTEdge<T>();
+  this->edges.append(e);
+  SymEdge<T> *se = &e->symedges[0];
+  SymEdge<T> *sesym = &e->symedges[1];
+  se->edge = sesym->edge = e;
+  se->face = fleft;
+  sesym->face = fright;
+  se->vert = v1;
+  if (v1->symedge == nullptr) {
+    v1->symedge = se;
+  }
+  sesym->vert = v2;
+  if (v2->symedge == nullptr) {
+    v2->symedge = sesym;
+  }
+  se->next = sesym->next = se->rot = sesym->rot = nullptr;
+  return e;
+}
+
+template<typename T> CDTFace<T> *CDTArrangement<T>::add_face()
+{
+  CDTFace<T> *f = new CDTFace<T>();
+  this->faces.append(f);
+  return f;
+}
+
+template<typename T> void CDTArrangement<T>::reserve(int num_verts, int num_edges, int num_faces)
+{
+  /* These reserves are just guesses; OK if they aren't exactly right since vectors will resize. */
+  this->verts.reserve(2 * num_verts);
+  this->edges.reserve(3 * num_verts + 2 * num_edges + 3 * 2 * num_faces);
+  this->faces.reserve(2 * num_verts + 2 * num_edges + 2 * num_faces);
+}
+
+template<typename T>
+CDT_state<T>::CDT_state(int num_input_verts, int num_input_edges, int num_input_faces, T epsilon)
+{
+  this->input_vert_tot = num_input_verts;
+  this->cdt.reserve(num_input_verts, num_input_edges, num_input_faces);
+  this->cdt.outer_face = this->cdt.add_face();
+  this->epsilon = epsilon;
+  this->visit_count = 0;
+}
+
+static bool id_in_list(const LinkNode *id_list, int id)
+{
+  const LinkNode *ln;
+
+  for (ln = id_list; ln != nullptr; ln = ln->next) {
+    if (POINTER_AS_INT(ln->link) == id) {
+      return true;
+    }
+  }
+  return false;
+}
+
+/* Is any id in (range_start, range_start+1, ... , range_end) in id_list? */
+static bool id_range_in_list(const LinkNode *id_list, int range_start, int range_end)
+{
+  const LinkNode *ln;
+  int id;
+
+  for (ln = id_list; ln != nullptr; ln = ln->next) {
+    id = POINTER_AS_INT(ln->link);
+    if (id >= range_start && id <= range_end) {
+      return true;
+    }
+  }
+  return false;
+}
+
+static void add_to_input_ids(LinkNode **dst, int input_id)
+{
+  if (!id_in_list(*dst, input_id)) {
+    BLI_linklist_prepend(dst, POINTER_FROM_INT(input_id));
+  }
+}
+
+static void add_list_to_input_ids(LinkNode **dst, const LinkNode *src)
+{
+  const LinkNode *ln;
+
+  for (ln = src; ln != nullptr; ln = ln->next) {
+    add_to_input_ids(dst, POINTER_AS_INT(ln->link));
+  }
+}
+
+template<typename T> inline bool is_border_edge(const CDTEdge<T> *e, const CDT_state<T> *cdt)
+{
+  return e->symedges[0].face == cdt->outer_face || e->symedges[1].face == cdt->outer_face;
+}
+
+template<typename T> inline bool is_constrained_edge(const CDTEdge<T> *e)
+{
+  return e->input_ids != NULL;
+}
+
+template<typename T> inline bool is_deleted_edge(const CDTEdge<T> *e)
+{
+  return e->symedges[0].next == NULL;
+}
+
+template<typename T> inline bool is_original_vert(const CDTVert<T> *v, CDT_state<T> *cdt)
+{
+  return (v->index < cdt->input_vert_tot);
+}
+
+/* Return the Symedge that goes from v1 to v2, if it exists, else return nullptr. */
+template<typename T>
+SymEdge<T> *find_symedge_between_verts(const CDTVert<T> *v1, const CDTVert<T> *v2)
+{
+  SymEdge<T> *t = v1->symedge;
+  SymEdge<T> *tstart = t;
+  do {
+    if (t->next->vert == v2) {
+      return t;
+    }
+  } while ((t = t->rot) != tstart);
+  return nullptr;
+}
+
+/**
+ * Return the SymEdge attached to v that has face f, if it exists, else return nullptr.
+ */
+template<typename T> SymEdge<T> *find_symedge_with_face(const CDTVert<T> *v, const CDTFace<T> *f)
+{
+  SymEdge<T> *t = v->symedge;
+  SymEdge<T> *tstart = t;
+  do {
+    if (t->face == f) {
+      return t;
+    }
+  } while ((t = t->rot) != tstart);
+  return nullptr;
+}
+
+/**
+ * Is there already an edge between a and b?
+ */
+template<typename T> inline bool exists_edge(const CDTVert<T> *v1, const CDTVert<T> *v2)
+{
+  return find_symedge_between_verts(v1, v2) != nullptr;
+}
+
+/**
+ * Is the vertex v incident on face f?
+ */
+template<typename T> bool vert_touches_face(const CDTVert<T> *v, const CDTFace<T> *f)
+{
+  SymEdge<T> *se = v->symedge;
+  do {
+    if (se->face == f) {
+      return true;
+    }
+  } while ((se = se->rot) != v->symedge);
+  return false;
+}
+
+/**
+ * Assume s1 and s2 are both #SymEdges in a face with > 3 sides,
+ * and one is not the next of the other.
+ * Add an edge from `s1->v` to `s2->v`, splitting the face in two.
+ * The original face will continue to be associated with the sub-face
+ * that has s1, and a new face will be made for s2's new face.
+ * Return the new diagonal's #CDTEdge pointer.
+ */
+template<typename T> CDTEdge<T> *CDTArrangement<T>::add_diagonal(SymEdge<T> *s1, SymEdge<T> *s2)
+{
+  CDTFace<T> *fold = s1->face;
+  CDTFace<T> *fnew = this->add_face();
+  SymEdge<T> *s1prev = prev(s1);
+  SymEdge<T> *s1prevsym = sym(s1prev);
+  SymEdge<T> *s2prev = prev(s2);
+  SymEdge<T> *s2prevsym = sym(s2prev);
+  CDTEdge<T> *ediag = this->add_edge(s1->vert, s2->vert, fnew, fold);
+  SymEdge<T> *sdiag = &ediag->symedges[0];
+  SymEdge<T> *sdiagsym = &ediag->symedges[1];
+  sdiag->next = s2;
+  sdiagsym->next = s1;
+  s2prev->next = sdiagsym;
+  s1prev->next = sdiag;
+  s1->rot = sdiag;
+  sdiag->rot = s1prevsym;
+  s2->rot = sdiagsym;
+  sdiagsym->rot = s2prevsym;
+  for (SymEdge<T> *se = s2; se != sdiag; se = se->next) {
+    se->face = fnew;
+  }
+  add_list_to_input_ids(&fnew->input_ids, fold->input_ids);
+  return ediag;
+}
+
+template<typename T>
+CDTEdge<T> *CDTArrangement<T>::add_vert_to_symedge_edge(CDTVert<T> *v, SymEdge<T> *se)
+{
+  SymEdge<T> *se_rot = se->rot;
+  SymEdge<T> *se_rotsym = sym(se_rot);
+  /* TODO: check: I think last arg in next should be sym(se)->face. */
+  CDTEdge<T> *e = this->add_edge(v, se->vert, se->face, se->face);
+  SymEdge<T> *new_se = &e->symedges[0];
+  SymEdge<T> *new_se_sym = &e->symedges[1];
+  new_se->next = se;
+  new_se_sym->next = new_se;
+  new_se->rot = new_se;
+  new_se_sym->rot = se_rot;
+  se->rot = new_se_sym;
+  se_rotsym->next = new_se_sym;
+  return e;
+}
+
+/**
+ * Connect the verts of se1 and se2, assuming that currently those two #SymEdge's are on
+ * the outer boundary (have face == outer_face) of two components that are isolated from
+ * each other.
+ */
+template<typename T>
+CDTEdge<T> *CDTArrangement<T>::connect_separate_parts(SymEdge<T> *se1, SymEdge<T> *se2)
+{
+  BLI_assert(se1->face == this->outer_face && se2->face == this->outer_face);
+  SymEdge<T> *se1_rot = se1->rot;
+  SymEdge<T> *se1_rotsym = sym(se1_rot);
+  SymEdge<T> *se2_rot = se2->rot;
+  SymEdge<T> *se2_rotsym = sym(se2_rot);
+  CDTEdge<T> *e = this->add_edge(se1->vert, se2->vert, this->outer_face, this->outer_face);
+  SymEdge<T> *new_se = &e->symedges[0];
+  SymEdge<T> *new_se_sym = &e->symedges[1];
+  new_se->next = se2;
+  new_se_sym->next = se1;
+  new_se->rot = se1_rot;
+  new_se_sym->rot = se2_rot;
+  se1->rot = new_se;
+  se2->rot = new_se_sym;
+  se1_rotsym->next = new_se;
+  se2_rotsym->next = new_se_sym;
+  return e;
+}
+
+/**
+ * Split se at fraction lambda,
+ * and return the new #CDTEdge that is the new second half.
+ * Copy the edge input_ids into the new one.
+ */
+template<typename T> CDTEdge<T> *CDTArrangement<T>::split_edge(SymEdge<T> *se, T lambda)
+{
+  /* Split e at lambda. */
+  const vec2<T> *a = &se->vert->co;
+  const vec2<T> *b = &se->next->vert->co;
+  SymEdge<T> *sesym = sym(se);
+  SymEdge<T> *sesymprev = prev(sesym);
+  SymEdge<T> *sesymprevsym = sym(sesymprev);
+  SymEdge<T> *senext = se->next;
+  CDTVert<T> *v = this->add_vert(vec2<T>::interpolate(*a, *b, lambda));
+  CDTEdge<T> *e = this->add_edge(v, se->next->vert, se->face, sesym->face);
+  sesym->vert = v;
+  SymEdge<T> *newse = &e->symedges[0];
+  SymEdge<T> *newsesym = &e->symedges[1];
+  se->next = newse;
+  newsesym->next = sesym;
+  newse->next = senext;
+  newse->rot = sesym;
+  sesym->rot = newse;
+  senext->rot = newsesym;
+  newsesym->rot = sesymprevsym;
+  sesymprev->next = newsesym;
+  if (newsesym->vert->symedge == sesym) {
+    newsesym->vert->symedge = newsesym;
+  }
+  add_list_to_input_ids(&e->input_ids, se->edge->input_ids);
+  return e;
+}
+
+/**
+ * Delete an edge from the structure. The new combined face on either side of
+ * the deleted edge will be the one that was e's face.
+ * There will be now an unused face, marked by setting its deleted flag,
+ * and an unused #CDTEdge, marked by setting the next and rot pointers of
+ * its #SymEdges to #nullptr.
+ * <pre>
+ *        .  v2               .
+ *       / \                 / \
+ *      /f|j\               /   \
+ *     /  |  \             /     \
+ *        |
+ *      A |  B                A
+ *    \  e|   /           \       /
+ *     \  | /              \     /
+ *      \h|i/               \   /
+ *        .  v1               .
+ * </pre>
+ * Also handle variant cases where one or both ends
+ * are attached only to e.
+ */
+template<typename T> void CDTArrangement<T>::delete_edge(SymEdge<T> *se)
+{
+  SymEdge<T> *sesym = sym(se);
+  CDTVert<T> *v1 = se->vert;
+  CDTVert<T> *v2 = sesym->vert;
+  CDTFace<T> *aface = se->face;
+  CDTFace<T> *bface = sesym->face;
+  SymEdge<T> *f = se->next;
+  SymEdge<T> *h = prev(se);
+  SymEdge<T> *i = sesym->next;
+  SymEdge<T> *j = prev(sesym);
+  SymEdge<T> *jsym = sym(j);
+  SymEdge<T> *hsym = sym(h);
+  bool v1_isolated = (i == se);
+  bool v2_isolated = (f == sesym);
+
+  if (!v1_isolated) {
+    h->next = i;
+    i->rot = hsym;
+  }
+  if (!v2_isolated) {
+    j->next = f;
+    f->rot = jsym;
+  }
+  if (!v1_isolated && !v2_isolated && aface != bface) {
+    for (SymEdge<T> *k = i; k != f; k = k->next) {
+      k->face = aface;
+    }
+  }
+
+  /* If e was representative symedge for v1 or v2, fix that. */
+  if (v1_isolated) {
+    v1->symedge = nullptr;
+  }
+  else if (v1->symedge == se) {
+    v1->symedge = i;
+  }
+  if (v2_isolated) {
+    v2->symedge = nullptr;
+  }
+  else if (v2->symedge == sesym) {
+    v2->symedge = f;
+  }
+
+  /* Mark SymEdge as deleted by setting all its pointers to NULL. */
+  se->next = se->rot = nullptr;
+  sesym->next = sesym->rot = nullptr;
+  if (!v1_isolated && !v2_isolated && aface != bface) {
+    bface->deleted = true;
+    if (this->outer_face == bface) {
+      this->outer_face = aface;
+    }
+  }
+}
+
+template<typename T> class SiteInfo {
+ public:
+  CDTVert<T> *v;
+  int orig_index;
+};
+
+/**
+ * Compare function for lexicographic sort: x, then y, then index.
+ */
+template<typename T> bool site_lexicographic_sort(const SiteInfo<T> &a, const SiteInfo<T> &b)
+{
+  const vec2<T> &co_a = a.v->co;
+  const vec2<T> &co_b = b.v->co;
+  if (co_a[0] < co_b[0]) {
+    return true;
+  }
+  if (co_a[0] > co_b[0]) {
+    return false;
+  }
+  if (co_a[1] < co_b[1]) {
+    return true;
+  }
+  if (co_a[1] > co_b[1]) {
+    return false;
+  }
+  return a.orig_index < b.orig_index;
+}
+
+/**
+ * Find series of equal vertices in the sorted sites array
+ * and use the vertices merge_to_index to indicate that
+ * all vertices after the first merge to the first.
+ */
+template<typename T> void find_site_merges(Array<SiteInfo<T>> &sites)
+{
+  int n = sites.size();
+  for (int i = 0; i < n - 1; ++i) {
+    int j = i + 1;
+    while (j < n && sites[j].v->co == sites[i].v->co) {
+      sites[j].v->merge_to_index = sites[i].orig_index;
+      ++j;
+    }
+    if (j - i > 1) {
+      i = j - 1; /* j-1 because loop head will add another 1. */
+    }
+  }
+}
+
+template<typename T> inline bool vert_left_of_symedge(CDTVert<T> *v, SymEdge<T> *se)
+{
+  return orient2d(v->co, se->vert->co, se->next->vert->co) > 0;
+}
+
+template<typename T> inline bool vert_right_of_symedge(CDTVert<T> *v, SymEdge<T> *se)
+{
+  return orient2d(v->co, se->next->vert->co, se->vert->co) > 0;
+}
+
+/* Is se above basel? */
+template<typename T>
+inline bool dc_tri_valid(SymEdge<T> *se, SymEdge<T> *basel, SymEdge<T> *basel_sym)
+{
+  return orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0;
+}
+
+/**
+ * Delaunay triangulate sites[start} to sites[end-1].
+ * Assume sites are lexicographically sorted by coordinate.
+ * Return #SymEdge of CCW convex hull at left-most point in *r_le
+ * and that of right-most point of cw convex null in *r_re.
+ */
+template<typename T>
+void dc_tri(CDTArrangement<T> *cdt,
+            Array<SiteInfo<T>> &sites,
+            int start,
+            int end,
+            SymEdge<T> **r_le,
+            SymEdge<T> **r_re)
+{
+  constexpr int dbg_level = 0;
+  if (dbg_level > 0) {
+    std::cout << "DC_TRI start=" << start << " end=" << end << "\n";
+  }
+  int n = end - start;
+  if (n <= 1) {
+    *r_le = nullptr;
+    *r_re = nullptr;
+    return;
+  }
+
+  /* Base case: if n <= 3, triangulate directly. */
+  if (n <= 3) {
+    CDTVert<T> *v1 = sites[start].v;
+    CDTVert<T> *v2 = sites[start + 1].v;
+    CDTEdge<T> *ea = cdt->add_edge(v1, v2, cdt->outer_face, cdt->outer_face);
+    ea->symedges[0].next = &ea->symedges[1];
+    ea->symedges[1].next = &ea->symedges[0];
+    ea->symedges[0].rot = &ea->symedges[0];
+    ea->symedges[1].rot = &ea->symedges[1];
+    if (n == 2) {
+      *r_le = &ea->symedges[0];
+      *r_re = &ea->symedges[1];
+      return;
+    }
+    CDTVert<T> *v3 = sites[start + 2].v;
+    CDTEdge<T> *eb = cdt->add_vert_to_symedge_edge(v3, &ea->symedges[1]);
+    int orient = orient2d(v1->co, v2->co, v3->co);
+    if (orient > 0) {
+      cdt->add_diagonal(&eb->symedges[0], &ea->symedges[0]);
+      *r_le = &ea->symedges[0];
+      *r_re = &eb->symedges[0];
+    }
+    else if (orient < 0) {
+      cdt->add_diagonal(&ea->symedges[0], &eb->symedges[0]);
+      *r_le = ea->symedges[0].rot;
+      *r_re = eb->symedges[0].rot;
+    }
+    else {
+      /* Collinear points. Just return a line. */
+      *r_le = &ea->symedges[0];
+      *r_re = &eb->symedges[0];
+    }
+    return;
+  }
+  /* Recursive case. Do left (L) and right (R) halves seperately, then join. */
+  int n2 = n / 2;
+  BLI_assert(n2 >= 2 && end - (start + n2) >= 2);
+  SymEdge<T> *ldo;
+  SymEdge<T> *ldi;
+  SymEdge<T> *rdi;
+  SymEdge<T> *rdo;
+  dc_tri(cdt, sites, start, start + n2, &ldo, &ldi);
+  dc_tri(cdt, sites, start + n2, end, &rdi, &rdo);
+  if (dbg_level > 0) {
+    std::cout << "\nDC_TRI merge step for start=" << start << ", end=" << end << "\n";
+    std::cout << "ldo " << ldo << "\n"
+              << "ldi " << ldi << "\n"
+              << "rdi " << rdi << "\n"
+              << "rdo " << rdo << "\n";
+    if (dbg_level > 1) {
+      std::string lab = "dc_tri (" + std::to_string(start) + "," + std::to_string(start + n2) +
+                        ")(" + std::to_string(start + n2) + "," + std::to_string(end) + ")";
+      cdt_draw(lab, *cdt);
+    }
+  }
+  /* Find lower common tangent of L and R. */
+  for (;;) {
+    if (vert_left_of_symedge(rdi->vert, ldi)) {
+      ldi = ldi->next;
+    }
+    else if (vert_right_of_symedge(ldi->vert, rdi)) {
+      rdi = sym(rdi)->rot; /* Previous edge to rdi with same right face. */
+    }
+    else {
+      break;
+    }
+  }
+  if (dbg_level > 0) {
+    std::cout << "common lower tangent in between\n"
+              << "rdi " << rdi << "\n"
+              << "ldi" << ldi << "\n";
+  }
+
+  CDTEdge<T> *ebasel = cdt->connect_separate_parts(sym(rdi)->next, ldi);
+  SymEdge<T> *basel = &ebasel->symedges[0];
+  SymEdge<T> *basel_sym = &ebasel->symedges[1];
+  if (dbg_level > 1) {
+    std::cout << "basel " << basel;
+    cdt_draw("after basel made", *cdt);
+  }
+  if (ldi->vert == ldo->vert) {
+    ldo = basel_sym;
+  }
+  if (rdi->vert == rdo->vert) {
+    rdo = basel;
+  }
+
+  /* Merge loop. */
+  for (;;) {
+    /* Locate the first point lcand->next->vert encountered by rising bubble,
+     * and delete L edges out of basel->next->vert that fail the circle test. */
+    SymEdge<T> *lcand = basel_sym->rot;
+    SymEdge<T> *rcand = basel_sym->next;
+    if (dbg_level > 1) {
+      std::cout << "\ntop of merge loop\n";
+      std::cout << "lcand " << lcand << "\n"
+                << "rcand " << rcand << "\n"
+                << "basel " << basel << "\n";
+    }
+    if (dc_tri_valid(lcand, basel, basel_sym)) {
+      if (dbg_level > 1) {
+        std::cout << "found valid lcand\n";
+        std::cout << "  lcand" << lcand << "\n";
+      }
+      while (incircle(basel_sym->vert->co,
+                      basel->vert->co,
+                      lcand->next->vert->co,
+                      lcand->rot->next->vert->co) > 0.0) {
+        if (dbg_level > 1) {
+          std::cout << "incircle says to remove lcand\n";
+          std::cout << "  lcand" << lcand << "\n";
+        }
+        SymEdge<T> *t = lcand->rot;
+        cdt->delete_edge(sym(lcand));
+        lcand = t;
+      }
+    }
+    /* Symmetrically, locate first R point to be hit and delete R edges. */
+    if (dc_tri_valid(rcand, basel, basel_sym)) {
+      if (dbg_level > 1) {
+        std::cout << "found valid rcand\n";
+        std::cout << "  rcand" << rcand << "\n";
+      }
+      while (incircle(basel_sym->vert->co,
+                      basel->vert->co,
+                      rcand->next->vert->co,
+                      sym(rcand)->next->next->vert->co) > 0.0) {
+        if (dbg_level > 0) {
+          std::cout << "incircle says to remove rcand\n";
+          std::cout << "  rcand" << rcand << "\n";
+        }
+        SymEdge<T> *t = sym(rcand)->next;
+        cdt->delete_edge(rcand);
+        rcand = t;
+      }
+    }
+    /* If both lcand and rcand are invalid, then basel is the common upper tangent. */
+    bool valid_lcand = dc_tri_valid(lcand, basel, basel_sym);
+    bool valid_rcand = dc_tri_valid(rcand, basel, basel_sym);
+    if (dbg_level > 0) {
+      std::cout << "after bubbling up, valid_lcand=" << valid_lcand
+                << ", valid_rand=" << valid_rcand << "\n";
+      std::cout << "lcand" << lcand << "\n"
+                << "rcand" << rcand << "\n";
+    }
+    if (!valid_lcand && !valid_rcand) {
+      break;
+    }
+    /* The next cross edge to be connected is to either `lcand->next->vert` or `rcand->next->vert`;
+     * if both are valid, choose the appropriate one using the #incircle test. */
+    if (!valid_lcand ||
+        (valid_rcand &&
+         incircle(lcand->next->vert->co, lcand->vert->co, rcand->vert->co, rcand->next->vert->co) >
+             0)) {
+      if (dbg_level > 0) {
+        std::cout << "connecting rcand\n";
+        std::cout << "  se1=basel_sym" << basel_sym << "\n";
+        std::cout << "  se2=rcand->next" << rcand->next << "\n";
+      }
+      ebasel = cdt->add_diagonal(rcand->next, basel_sym);
+    }
+    else {
+      if (dbg_level > 0) {
+        std::cout << "connecting lcand\n";
+        std::cout << "  se1=sym(lcand)" << sym(lcand) << "\n";
+        std::cout << "  se2=basel_sym->next" << basel_sym->next << "\n";
+      }
+      ebasel = cdt->add_diagonal(basel_sym->next, sym(lcand));
+    }
+    basel = &ebasel->symedges[0];
+    basel_sym = &ebasel->symedges[1];
+    BLI_assert(basel_sym->face == cdt->outer_face);
+    if (dbg_level > 2) {
+      cdt_draw("after adding new crossedge", *cdt);
+    }
+  }
+  *r_le = ldo;
+  *r_re = rdo;
+  BLI_assert(sym(ldo)->face == cdt->outer_face && rdo->face == cdt->outer_face);
+}
+
+/* Guibas-Stolfi Divide-and_Conquer algorithm. */
+template<typename T> void dc_triangulate(CDTArrangement<T> *cdt, Array<SiteInfo<T>> &sites)
+{
+  /* Compress sites in place to eliminted verts that merge to others. */
+  int i = 0;
+  int j = 0;
+  int nsites = sites.size();
+  while (j < nsites) {
+    /* Invariante: sites[0..i-1] have non-merged verts from 0..(j-1) in them. */
+    sites[i] = sites[j++];
+    if (sites[i].v->merge_to_index < 0) {
+      i++;
+    }
+  }
+  int n = i;
+  if (n == 0) {
+    return;
+  }
+  SymEdge<T> *le;
+  SymEdge<T> *re;
+  dc_tri(cdt, sites, 0, n, &le, &re);
+}
+
+/**
+ * Do a Delaunay Triangulation of the points in cdt.verts.
+ * This  is only a first step in the Constrained Delaunay triangulation,
+ * because it doesn't yet deal with the segment constraints.
+ * The algorithm used is the Divide & Conquer algorithm from the
+ * Guibas-Stolfi "Primitives for the Manipulation of General Subdivision
+ * and the Computation of Voronoi Diagrams" paper.
+ * The data structure here is similar to but not exactly the same as
+ * the quad-edge structure described in that paper.
+ * If T is not exact arithmetic, incircle and CCW tests are done using
+ * Shewchuk's exact primitives, so that this routine is robust.
+ *
+ * As a preprocessing step, we want to merge all vertices that the same.
+ * This is accomplished by lexicographically
+ * sorting the coordinates first (which is needed anyway for the D&C algorithm).
+ * The CDTVerts with merge_to_index not equal to -1 are after this regarded
+ * as having been merged into the vertex with the corresponding index.
+ */
+template<typename T> void initial_triangulation(CDTArrangement<T> *cdt)
+{
+  int n = cdt->verts.size();
+  if (n <= 1) {
+    return;
+  }
+  Array<SiteInfo<T>> sites(n);
+  for (int i = 0; i < n; ++i) {
+    sites[i].v = cdt->verts[i];
+    sites[i].orig_index = i;
+  }
+  std::sort(sites.begin(), sites.end(), site_lexicographic_sort<T>);
+  find_site_merges(sites);
+  dc_triangulate(cdt, sites);
+}
+
+/**
+ * Re-triangulates, assuring constrained delaunay condition,
+ * the pseudo-polygon that cycles from se.
+ * "pseudo" because a vertex may be repeated.
+ * See Anglada paper, "An Improved incremental algorithm
+ * for constructing restricted Delaunay triangulations".
+ */
+template<typename T> static void re_delaunay_triangulate(CDTArrangement<T> *cdt, SymEdge<T> *se)
+{
+  if (se->face == cdt->outer_face || sym(se)->face == cdt->outer_face) {
+    return;
+  }
+  /* 'se' is a diagonal just added, and it is base of area to retriangulate (face on its left) */
+  int count = 1;
+  for (SymEdge<T> *ss = se->next; ss != se; ss = ss->next) {
+    count++;
+  }
+  if (count <= 3) {
+    return;
+  }
+  /* First and last are the SymEdges whose verts are first and last off of base,
+   * continuing from 'se'. */
+  SymEdge<T> *first = se->next->next;
+  /* We want to make a triangle with 'se' as base and some other c as 3rd vertex. */
+  CDTVert<T> *a = se->vert;
+  CDTVert<T> *b = se->next->vert;
+  CDTVert<T> *c = first->vert;
+  SymEdge<T> *cse = first;
+  for (SymEdge<T> *ss = first->next; ss != se; ss = ss->next) {
+    CDTVert<T> *v = ss->vert;
+    if (incircle(a->co, b->co, c->co, v->co) > 0) {
+      c = v;
+      cse = ss;
+    }
+  }
+  /* Add diagonals necessary to make abc a triangle. */
+  CDTEdge<T> *ebc = nullptr;
+  CDTEdge<T> *eca = nullptr;
+  if (!exists_edge(b, c)) {
+    ebc = cdt->add_diagonal(se->next, cse);
+  }
+  if (!exists_edge(c, a)) {
+    eca = cdt->add_diagonal(cse, se);
+  }
+  /* Now recurse. */
+  if (ebc) {
+    re_delaunay_triangulate(cdt, &ebc->symedges[1]);
+  }
+  if (eca) {
+    re_delaunay_triangulate(cdt, &eca->symedges[1]);
+  }
+}
+
+template<typename T> inline int tri_orient(const SymEdge<T> *t)
+{
+  return orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co);
+}
+
+/**
+ * The #CrossData class 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.]
+ *
+ * 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.
+ *
+ * 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' will 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.
+ */
+template<typename T> class CrossData {
+ public:
+  T lambda = T(0);
+  CDTVert<T> *vert;
+  SymEdge<T> *in;
+  SymEdge<T> *out;
+
+  CrossData() : lambda(T(0)), vert(nullptr), in(nullptr), out(nullptr)
+  {
+  }
+  CrossData(T l, CDTVert<T> *v, SymEdge<T> *i, SymEdge<T> *o) : lambda(l), vert(v), in(i), out(o)
+  {
+  }
+  CrossData(const CrossData &other)
+      : lambda(other.lambda), vert(other.vert), in(other.in), out(other.out)
+  {
+  }
+  CrossData(CrossData &&other) noexcept
+      : lambda(std::move(other.lambda)),
+        vert(std::move(other.vert)),
+        in(std::move(other.in)),
+        out(std::move(other.out))
+  {
+  }
+  ~CrossData() = default;
+  CrossData &operator=(const CrossData &other)
+  {
+    if (this != &other) {
+      lambda = other.lambda;
+      vert = other.vert;
+      in = other.in;
+      out = other.out;
+    }
+    return *this;
+  }
+  CrossData &operator=(CrossData &&other) noexcept
+  {
+    lambda = std::move(other.lambda);
+    vert = std::move(other.vert);
+    in = std::move(other.in);
+    out = std::move(other.out);
+    return *this;
+  }
+};
+
+template<typename T>
+bool get_next_crossing_from_vert(CDT_state<T> *cdt_state,
+                                 CrossData<T> *cd,
+                                 CrossData<T> *cd_next,
+                                 const CDTVert<T> *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.
+ */
+template<typename T>
+void fill_crossdata_for_through_vert(CDTVert<T> *v,
+                                     SymEdge<T> *cd_out,
+                                     CrossData<T> *cd,
+                                     CrossData<T> *cd_next)
+{
+  SymEdge<T> *se;
+
+  cd_next->lambda = T(0);
+  cd_next->vert = v;
+  cd_next->in = NULL;
+  cd_next->out = NULL;
+  if (cd->lambda == 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;
+  }
+}
+
+/**
+ * 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 \a curco to \a 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.
+ */
+template<typename T>
+void fill_crossdata_for_intersect(const vec2<T> &curco,
+                                  const CDTVert<T> *v2,
+                                  SymEdge<T> *t,
+                                  CrossData<T> *cd,
+                                  CrossData<T> *cd_next,
+                                  const T epsilon)
+{
+  CDTVert<T> *va = t->vert;
+  CDTVert<T> *vb = t->next->vert;
+  CDTVert<T> *vc = t->next->next->vert;
+  SymEdge<T> *se_vcvb = sym(t->next);
+  SymEdge<T> *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);
+  UNUSED_VARS_NDEBUG(vc);
+  auto isect = vec2<T>::isect_seg_seg(va->co, vb->co, curco, v2->co);
+  T &lambda = isect.lambda;
+  switch (isect.kind) {
+    case vec2<T>::isect_result::LINE_LINE_CROSS: {
+#ifdef WITH_GMP
+      if (!std::is_same<T, mpq_class>::value) {
+#else
+      if (true) {
+#endif
+        T len_ab = vec2<T>::distance(va->co, vb->co);
+        if (lambda * len_ab <= epsilon) {
+          fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
+        }
+        else if ((1 - lambda) * len_ab <= epsilon) {
+          fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
+        }
+        else {
+          *cd_next = CrossData<T>(lambda, nullptr, t, nullptr);
+          if (cd->lambda == 0) {
+            cd->out = se_vcva;
+          }
+        }
+      }
+      else {
+        *cd_next = CrossData<T>(lambda, nullptr, t, nullptr);
+        if (cd->lambda == 0) {
+          cd->out = se_vcva;
+        }
+      }
+      break;
+    }
+    case vec2<T>::isect_result::LINE_LINE_EXACT: {
+      if (lambda == 0) {
+        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
+      }
+      else if (lambda == 1) {
+        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
+      }
+      else {
+        *cd_next = CrossData<T>(lambda, nullptr, t, nullptr);
+        if (cd->lambda == 0) {
+          cd->out = se_vcva;
+        }
+      }
+      break;
+    }
+    case vec2<T>::isect_result::LINE_LINE_NONE: {
+#ifdef WITH_GMP
+      if (std::is_same<T, mpq_class>::value) {
+        BLI_assert(false);
+      }
+#endif
+      /* It should be very near one end or other of segment. */
+      const T middle_lambda = 0.5;
+      if (lambda <= middle_lambda) {
+        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 vec2<T>::isect_result::LINE_LINE_COLINEAR: {
+      if (vec2<T>::distance_squared(va->co, v2->co) <= vec2<T>::distance_squared(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;
+    }
+  }
+}  // namespace blender::meshintersect
+
+/**
+ * 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.
+ */
+template<typename T>
+bool get_next_crossing_from_vert(CDT_state<T> *cdt_state,
+                                 CrossData<T> *cd,
+                                 CrossData<T> *cd_next,
+                                 const CDTVert<T> *v2)
+{
+  SymEdge<T> *tstart = cd->vert->symedge;
+  SymEdge<T> *t = tstart;
+  CDTVert<T> *vcur = cd->vert;
+  bool 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_state->cdt.outer_face && tri_orient(t) < 0) {
+      BLI_assert(false); /* Shouldn't happen. */
+    }
+    CDTVert<T> *va = t->next->vert;
+    CDTVert<T> *vb = t->next->next->vert;
+    int orient1 = orient2d(t->vert->co, va->co, v2->co);
+    if (orient1 == 0 && in_line<T>(vcur->co, va->co, v2->co)) {
+      fill_crossdata_for_through_vert(va, t, cd, cd_next);
+      ok = true;
+      break;
+    }
+    if (t->face != cdt_state->cdt.outer_face) {
+      int orient2 = orient2d(vcur->co, vb->co, v2->co);
+      /* Don't handle orient2 == 0 case here: next rotation will get it. */
+      if (orient1 > 0 && orient2 < 0) {
+        /* Segment intersection. */
+        t = t->next;
+        fill_crossdata_for_intersect(vcur->co, v2, t, cd, cd_next, cdt_state->epsilon);
+        ok = true;
+        break;
+      }
+    }
+  } while ((t = t->rot) != tstart);
+  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`.
+ */
+template<typename T>
+void get_next_crossing_from_edge(CrossData<T> *cd,
+                                 CrossData<T> *cd_next,
+                                 const CDTVert<T> *v2,
+                                 const T epsilon)
+{
+  CDTVert<T> *va = cd->in->vert;
+  CDTVert<T> *vb = cd->in->next->vert;
+  vec2<T> curco = vec2<T>::interpolate(va->co, vb->co, cd->lambda);
+  SymEdge<T> *se_ac = sym(cd->in)->next;
+  CDTVert<T> *vc = se_ac->next->vert;
+  int orient = orient2d(curco, v2->co, vc->co);
+  if (orient < 0) {
+    fill_crossdata_for_intersect<T>(curco, v2, se_ac->next, cd, cd_next, epsilon);
+  }
+  else if (orient > 0.0) {
+    fill_crossdata_for_intersect(curco, v2, se_ac, cd, cd_next, epsilon);
+  }
+  else {
+    *cd_next = CrossData<T>{0.0, vc, se_ac->next, nullptr};
+  }
+}
+
+constexpr int inline_crossings_size = 128;
+template<typename T>
+void dump_crossings(const Vector<CrossData<T>, inline_crossings_size> &crossings)
+{
+  std::cout << "CROSSINGS\n";
+  for (int i = 0; i < crossings.size(); ++i) {
+    std::cout << i << ": ";
+    const CrossData<T> &cd = crossings[i];
+    if (cd.lambda == 0) {
+      std::cout << "v" << cd.vert->index;
+    }
+    else {
+      std::cout << "lambda=" << cd.lambda;
+    }
+    if (cd.in != nullptr) {
+      std::cout << " in=" << short_se_dump(cd.in);
+      std::cout << " out=" << short_se_dump(cd.out);
+    }
+    std::cout << "\n";
+  }
+}
+
+/**
+ * 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 #blender_constrained_delaunay_get_output has not been called yet.
+ */
+template<typename T>
+void add_edge_constraint(
+    CDT_state<T> *cdt_state, CDTVert<T> *v1, CDTVert<T> *v2, int input_id, LinkNode **r_edges)
+{
+  constexpr int dbg_level = 0;
+  if (dbg_level > 0) {
+    std::cout << "\nADD EDGE CONSTRAINT\n" << vertname(v1) << " " << vertname(v2) << "\n";
+  }
+  LinkNodePair edge_list = {NULL, NULL};
+
+  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;
+  }
+  SymEdge<T> *t = find_symedge_between_verts(v1, v2);
+  if (t != nullptr) {
+    /* Segment already there. */
+    add_to_input_ids(&t->edge->input_ids, input_id);
+    if (r_edges != NULL) {
+      BLI_linklist_append(&edge_list, t->edge);
+      *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.
+   */
+  int visit = ++cdt_state->visit_count;
+  Vector<CrossData<T>, inline_crossings_size> crossings;
+  crossings.append(CrossData<T>(T(0), v1, nullptr, nullptr));
+  int n;
+  while (!((n = crossings.size()) > 0 && crossings[n - 1].vert == v2)) {
+    crossings.append(CrossData<T>());
+    CrossData<T> *cd = &crossings[n - 1];
+    CrossData<T> *cd_next = &crossings[n];
+    bool ok;
+    if (crossings[n - 1].lambda == 0) {
+      ok = get_next_crossing_from_vert(cdt_state, cd, cd_next, v2);
+    }
+    else {
+      get_next_crossing_from_edge(cd, cd_next, v2, cdt_state->epsilon);
+      ok = true;
+    }
+    constexpr int unreasonably_large_crossings = 100000;
+    if (!ok || crossings.size() == unreasonably_large_crossings) {
+      /* Shouldn't happen but if does, just bail out. */
+      BLI_assert(false);
+      return;
+    }
+    if (crossings[n].lambda == 0) {
+      if (crossings[n].vert->visit_index == visit) {
+        /* Shouldn't happen but if it does, just bail out. */
+        BLI_assert(false);
+        return;
+      }
+      crossings[n].vert->visit_index = visit;
+    }
+  }
+
+  if (dbg_level > 0) {
+    dump_crossings(crossings);
+  }
+
+  /*
+   * Post-process 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.
+   */
+  int ncrossings = crossings.size();
+  for (int i = 2; i < ncrossings; ++i) {
+    CrossData<T> *cd = &crossings[i];
+    if (cd->lambda == 0.0) {
+      CDTVert<T> *v = cd->vert;
+      int j;
+      CrossData<T> *cd_prev;
+      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;
+        }
+        cd_prev->lambda = -1.0; /* Mark cd_prev as 'deleted'. */
+      }
+      if (j < i - 1) {
+        /* Some crossings were deleted. Fix the in and out edges across gap. */
+        cd_prev = &crossings[j];
+        SymEdge<T> *se;
+        if (cd_prev->lambda == 0.0) {
+          se = find_symedge_between_verts(cd_prev->vert, v);
+          if (se == NULL) {
+            return;
+          }
+          cd_prev->out = se;
+          cd->in = NULL;
+        }
+        else {
+          se = find_symedge_with_face(v, sym(cd_prev->in)->face);
+          if (se == NULL) {
+            return;
+          }
+          cd->in = se;
+        }
+      }
+    }
+  }
+
+  /*
+   * Insert all intersection points on constrained edges.
+   */
+  for (int i = 0; i < ncrossings; ++i) {
+    CrossData<T> *cd = &crossings[i];
+    if (cd->lambda != 0.0 && cd->lambda != -1.0 && is_constrained_edge(cd->in->edge)) {
+      CDTEdge<T> *edge = cdt_state->cdt.split_edge(cd->in, cd->lambda);
+      cd->vert = edge->symedges[0].vert;
+    }
+  }
+
+  /*
+   * Remove any crossed, non-intersected edges.
+   */
+  for (int i = 0; i < ncrossings; ++i) {
+    CrossData<T> *cd = &crossings[i];
+    if (cd->lambda != 0.0 && cd->lambda != -1.0 && !is_constrained_edge(cd->in->edge)) {
+      cdt_state->cdt.delete_edge(cd->in);
+    }
+  }
+
+  /*
+   * Insert segments for v1->v2.
+   */
+  SymEdge<T> *tstart = crossings[0].out;
+  for (int i = 1; i < ncrossings; i++) {
+    CrossData<T> *cd = &crossings[i];
+    if (cd->lambda == -1.0) {
+      continue; /* This crossing was deleted. */
+    }
+    t = NULL;
+    SymEdge<T> *tnext = t;
+    CDTEdge<T> *edge;
+    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) {
+        /* Previous non-deleted crossing must also have been a vert, and segment should exist. */
+        int j;
+        CrossData<T> *cd_prev;
+        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);
+        edge = cd_prev->out->edge;
+        add_to_input_ids(&edge->input_ids, input_id);
+        if (r_edges != NULL) {
+          BLI_linklist_append(&edge_list, edge);
+        }
+      }
+    }
+    if (t != NULL) {
+      if (tstart->next->vert == t->vert) {
+        edge = tstart->edge;
+      }
+      else {
+        edge = cdt_state->cdt.add_diagonal(tstart, t);
+      }
+      add_to_input_ids(&edge->input_ids, input_id);
+      if (r_edges != NULL) {
+        BLI_linklist_append(&edge_list, edge);
+      }
+      /* Now retriangulate upper and lower gaps. */
+      re_delaunay_triangulate(&cdt_state->cdt, &edge->symedges[0]);
+      re_delaunay_triangulate(&cdt_state->cdt, &edge->symedges[1]);
+    }
+    if (i < ncrossings - 1) {
+      if (tnext != NULL) {
+        tstart = tnext;
+      }
+    }
+  }
+
+  if (r_edges) {
+    *r_edges = edge_list.list;
+  }
+}
+
+/**
+ * Incrementally add edge input edge as a constraint. This may cause the graph structure
+ * to change, in cases where the constraints intersect existing edges.
+ * The code will ensure that #CDTEdge's created will have ids that tie them back
+ * to the original edge constraint index.
+ */
+template<typename T> void add_edge_constraints(CDT_state<T> *cdt_state, const CDT_input<T> &input)
+{
+  int ne = input.edge.size();
+  int nv = input.vert.size();
+  for (int i = 0; i < ne; i++) {
+    int iv1 = input.edge[i].first;
+    int iv2 = input.edge[i].second;
+    if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
+      /* Ignore invalid indices in edges. */
+      continue;
+    }
+    CDTVert<T> *v1 = cdt_state->cdt.get_vert_resolve_merge(iv1);
+    CDTVert<T> *v2 = cdt_state->cdt.get_vert_resolve_merge(iv2);
+    add_edge_constraint(cdt_state, v1, v2, i, nullptr);
+  }
+  cdt_state->face_edge_offset = ne;
+}
+
+/**
+ * Add face_id to the input_ids lists of all #CDTFace's on the interior of the input face with that
+ * id. face_symedge is on edge of the boundary of the input face, with assumption that interior is
+ * on the left of that #SymEdge.
+ *
+ * The algorithm is: starting from the #CDTFace for face_symedge, add the face_id and then
+ * process all adjacent faces where the adjacency isn't across an edge that was a constraint added
+ * for the boundary of the input face.
+ * fedge_start..fedge_end is the inclusive range of edge input ids that are for the given face.
+ *
+ * Note: if the input face is not CCW oriented, we'll be labeling the outside, not the inside.
+ * Note 2: if the boundary has self-crossings, this method will arbitrarily pick one of the
+ * contiguous set of faces enclosed by parts of the boundary, leaving the other such un-tagged.
+ * This may be a feature instead of a bug if the first contiguous section is most of the face and
+ * the others are tiny self-crossing triangles at some parts of the boundary.
+ * On the other hand, if decide we want to handle these in full generality, then will need a more
+ * complicated algorithm (using "inside" tests and a parity rule) to decide on the interior.
+ */
+template<typename T>
+void add_face_ids(
+    CDT_state<T> *cdt_state, SymEdge<T> *face_symedge, int face_id, int fedge_start, int fedge_end)
+{
+  /* Can't loop forever since eventually would visit every face. */
+  cdt_state->visit_count++;
+  int visit = cdt_state->visit_count;
+  Vector<SymEdge<T> *> stack;
+  stack.append(face_symedge);
+  while (!stack.is_empty()) {
+    SymEdge<T> *se = stack.pop_last();
+    CDTFace<T> *face = se->face;
+    if (face->visit_index == visit) {
+      continue;
+    }
+    face->visit_index = visit;
+    add_to_input_ids(&face->input_ids, face_id);
+    SymEdge<T> *se_start = se;
+    for (se = se->next; se != se_start; se = se->next) {
+      if (!id_range_in_list(se->edge->input_ids, fedge_start, fedge_end)) {
+        SymEdge<T> *se_sym = sym(se);
+        CDTFace<T> *face_other = se_sym->face;
+        if (face_other->visit_index != visit) {
+          stack.append(se_sym);
+        }
+      }
+    }
+  }
+}
+
+/* Return a power of 10 that is greater than or equal to x. */
+static int power_of_10_greater_equal_to(int x)
+{
+  if (x <= 0) {
+    return 1;
+  }
+  int ans = 1;
+  BLI_assert(x < INT_MAX / 10);
+  while (ans < x) {
+    ans *= 10;
+  }
+  return ans;
+}
+
+/**
+   Incrementally each edge of each input face as an edge constraint.
+ * The code will ensure that the #CDTEdge's created will have ids that tie them
+ * back to the original face edge (using a numbering system for those edges
+ * that starts with cdt->face_edge_offset, and continues with the edges in
+ * order around each face in turn. And then the next face starts at
+ * cdt->face_edge_offset beyond the start for the previous face.
+ */
+template<typename T> void add_face_constraints(CDT_state<T> *cdt_state, const CDT_input<T> &input)
+{
+  int nv = input.vert.size();
+  int nf = input.face.size();
+  int fstart = 0;
+  SymEdge<T> *face_symedge0 = nullptr;
+  CDTArrangement<T> *cdt = &cdt_state->cdt;
+  int maxflen = 0;
+  for (int f = 0; f < nf; f++) {
+    maxflen = max_ii(maxflen, input.face[f].size());
+  }
+  /* For convenience in debugging, make face_edge_offset be a power of 10. */
+  cdt_state->face_edge_offset = power_of_10_greater_equal_to(
+      max_ii(maxflen, cdt_state->face_edge_offset));
+  /* The original_edge encoding scheme doesn't work if the following is false.
+   * If we really have that many faces and that large a max face length that when multiplied
+   * together the are >= INT_MAX, then the Delaunay calculation will take unreasonably long anyway.
+   */
+  BLI_assert(INT_MAX / cdt_state->face_edge_offset > nf);
+  for (int f = 0; f < nf; f++) {
+    int flen = input.face[f].size();
+    if (flen <= 2) {
+      /* Ignore faces with fewer than 3 vertices. */
+      fstart += flen;
+      continue;
+    }
+    int fedge_start = (f + 1) * cdt_state->face_edge_offset;
+    for (int i = 0; i < flen; i++) {
+      int face_edge_id = fedge_start + i;
+      int iv1 = input.face[f][i];
+      int iv2 = input.face[f][(i + 1) % flen];
+      if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
+        /* Ignore face edges with invalid vertices. */
+        continue;
+      }
+      CDTVert<T> *v1 = cdt->get_vert_resolve_merge(iv1);
+      CDTVert<T> *v2 = cdt->get_vert_resolve_merge(iv2);
+      LinkNode *edge_list;
+      add_edge_constraint(cdt_state, v1, v2, face_edge_id, &edge_list);
+      /* Set a new face_symedge0 each time since earlier ones may not
+       * survive later symedge splits. Really, just want the one when
+       * i == flen -1, but this code guards against that one somehow
+       * being null.
+       */
+      if (edge_list != nullptr) {
+        CDTEdge<T> *face_edge = static_cast<CDTEdge<T> *>(edge_list->link);
+        face_symedge0 = &face_edge->symedges[0];
+        if (face_symedge0->vert != v1) {
+          face_symedge0 = &face_edge->symedges[1];
+          BLI_assert(face_symedge0->vert == v1);
+        }
+      }
+      BLI_linklist_free(edge_list, nullptr);
+    }
+    int fedge_end = fedge_start + flen - 1;
+    if (face_symedge0 != nullptr) {
+      add_face_ids(cdt_state, face_symedge0, f, fedge_start, fedge_end);
+    }
+    fstart += flen;
+  }
+}
+
+/* Delete_edge but try not to mess up outer face.
+ * Also faces have symedges now, so make sure not
+ * to mess those up either. */
+template<typename T> void dissolve_symedge(CDT_state<T> *cdt_state, SymEdge<T> *se)
+{
+  CDTArrangement<T> *cdt = &cdt_state->cdt;
+  SymEdge<T> *symse = sym(se);
+  if (symse->face == cdt->outer_face) {
+    se = sym(se);
+    symse = sym(se);
+  }
+  if (cdt->outer_face->symedge == se || cdt->outer_face->symedge == symse) {
+    /* Advancing by 2 to get past possible 'sym(se)'. */
+    if (se->next->next == se) {
+      cdt->outer_face->symedge = NULL;
+    }
+    else {
+      cdt->outer_face->symedge = se->next->next;
+    }
+  }
+  else {
+    if (se->face->symedge == se) {
+      se->face->symedge = se->next;
+    }
+    if (symse->face->symedge == symse) {
+      symse->face->symedge = symse->next;
+    }
+  }
+  cdt->delete_edge(se);
+}
+
+/**
+ * Remove all non-constraint edges.
+ */
+template<typename T> void remove_non_constraint_edges(CDT_state<T> *cdt_state)
+{
+  for (CDTEdge<T> *e : cdt_state->cdt.edges) {
+    SymEdge<T> *se = &e->symedges[0];
+    if (!is_deleted_edge(e) && !is_constrained_edge(e)) {
+      dissolve_symedge(cdt_state, se);
+    }
+  }
+}
+
+/*
+ * Remove the non-constraint edges, but leave enough of them so that all of the
+ * faces that would be #BMesh faces (that is, the faces that have some input representative)
+ * are valid: they can't have holes, they can't have repeated vertices, and they can't have
+ * repeated edges.
+ *
+ * Not essential, but to make the result look more aesthetically nice,
+ * remove the edges in order of decreasing length, so that it is more likely that the
+ * final remaining support edges are short, and therefore likely to make a fairly
+ * direct path from an outer face to an inner hole face.
+ */
+
+/**
+ * For sorting edges by decreasing length (squared).
+ */
+template<typename T> struct EdgeToSort {
+  T len_squared = T(0);
+  CDTEdge<T> *e{nullptr};
+
+  EdgeToSort() = default;
+  EdgeToSort(const EdgeToSort &other) : len_squared(other.len_squared), e(other.e)
+  {
+  }
+  EdgeToSort(EdgeToSort &&other) noexcept : len_squared(std::move(other.len_squared)), e(other.e)
+  {
+  }
+  ~EdgeToSort() = default;
+  EdgeToSort &operator=(const EdgeToSort &other)
+  {
+    if (this != &other) {
+      len_squared = other.len_squared;
+      e = other.e;
+    }
+    return *this;
+  }
+  EdgeToSort &operator=(EdgeToSort &&other)
+  {
+    len_squared = std::move(other.len_squared);
+    e = other.e;
+    return *this;
+  }
+};
+
+template<typename T> void remove_non_constraint_edges_leave_valid_bmesh(CDT_state<T> *cdt_state)
+{
+  CDTArrangement<T> *cdt = &cdt_state->cdt;
+  size_t nedges = cdt->edges.size();
+  if (nedges == 0) {
+    return;
+  }
+  Vector<EdgeToSort<T>> dissolvable_edges;
+  dissolvable_edges.reserve(cdt->edges.size());
+  int i = 0;
+  for (CDTEdge<T> *e : cdt->edges) {
+    if (!is_deleted_edge(e) && !is_constrained_edge(e)) {
+      dissolvable_edges.append(EdgeToSort<T>());
+      dissolvable_edges[i].e = e;
+      const vec2<T> &co1 = e->symedges[0].vert->co;
+      const vec2<T> &co2 = e->symedges[1].vert->co;
+      dissolvable_edges[i].len_squared = vec2<T>::distance_squared(co1, co2);
+      i++;
+    }
+  }
+  std::sort(dissolvable_edges.begin(),
+            dissolvable_edges.end(),
+            [](const EdgeToSort<T> &a, const EdgeToSort<T> &b) -> bool {
+              return (a.len_squared < b.len_squared);
+            });
+  for (EdgeToSort<T> &ets : dissolvable_edges) {
+    CDTEdge<T> *e = ets.e;
+    SymEdge<T> *se = &e->symedges[0];
+    bool dissolve = true;
+    CDTFace<T> *fleft = se->face;
+    CDTFace<T> *fright = sym(se)->face;
+    if (fleft != cdt->outer_face && fright != cdt->outer_face &&
+        (fleft->input_ids != nullptr || fright->input_ids != nullptr)) {
+      /* Is there another #SymEdge with same left and right faces?
+       * Or is there a vertex not part of e touching the same left and right faces? */
+      for (SymEdge<T> *se2 = se->next; dissolve && se2 != se; se2 = se2->next) {
+        if (sym(se2)->face == fright ||
+            (se2->vert != se->next->vert && vert_touches_face(se2->vert, fright))) {
+          dissolve = false;
+        }
+      }
+    }
+
+    if (dissolve) {
+      dissolve_symedge(cdt_state, se);
+    }
+  }
+}
+
+template<typename T> void remove_outer_edges_until_constraints(CDT_state<T> *cdt_state)
+{
+  // LinkNode *fstack = NULL;
+  // SymEdge *se, *se_start;
+  // CDTFace *f, *fsym;
+  int visit = ++cdt_state->visit_count;
+
+  cdt_state->cdt.outer_face->visit_index = visit;
+  /* Walk around outer face, adding faces on other side of dissolvable edges to stack. */
+  Vector<CDTFace<T> *> fstack;
+  SymEdge<T> *se_start = cdt_state->cdt.outer_face->symedge;
+  SymEdge<T> *se = se_start;
+  do {
+    if (!is_constrained_edge(se->edge)) {
+      CDTFace<T> *fsym = sym(se)->face;
+      if (fsym->visit_index != visit) {
+        fstack.append(fsym);
+      }
+    }
+  } while ((se = se->next) != se_start);
+
+  while (!fstack.is_empty()) {
+    LinkNode *to_dissolve = nullptr;
+    bool dissolvable;
+    CDTFace<T> *f = fstack.pop_last();
+    if (f->visit_index == visit) {
+      continue;
+    }
+    BLI_assert(f != cdt_state->cdt.outer_face);
+    f->visit_index = visit;
+    se_start = se = f->symedge;
+    do {
+      dissolvable = !is_constrained_edge(se->edge);
+      if (dissolvable) {
+        CDTFace<T> *fsym = sym(se)->face;
+        if (fsym->visit_index != visit) {
+          fstack.append(fsym);
+        }
+        else {
+          BLI_linklist_prepend(&to_dissolve, se);
+        }
+      }
+      se = se->next;
+    } while (se != se_start);
+    while (to_dissolve != NULL) {
+      se = static_cast<SymEdge<T> *>(BLI_linklist_pop(&to_dissolve));
+      if (se->next != NULL) {
+        dissolve_symedge(cdt_state, se);
+      }
+    }
+  }
+}
+
+/**
+ * Remove edges and merge faces to get desired output, as per options.
+ * \note the cdt cannot be further changed after this.
+ */
+template<typename T>
+void prepare_cdt_for_output(CDT_state<T> *cdt_state, const CDT_output_type output_type)
+{
+  CDTArrangement<T> *cdt = &cdt_state->cdt;
+  if (cdt->edges.is_empty()) {
+    return;
+  }
+
+  /* Make sure all non-deleted faces have a symedge. */
+  for (CDTEdge<T> *e : cdt->edges) {
+    if (!is_deleted_edge(e)) {
+      if (e->symedges[0].face->symedge == nullptr) {
+        e->symedges[0].face->symedge = &e->symedges[0];
+      }
+      if (e->symedges[1].face->symedge == nullptr) {
+        e->symedges[1].face->symedge = &e->symedges[1];
+      }
+    }
+  }
+
+  if (output_type == CDT_CONSTRAINTS) {
+    remove_non_constraint_edges(cdt_state);
+  }
+  else if (output_type == CDT_CONSTRAINTS_VALID_BMESH) {
+    remove_non_constraint_edges_leave_valid_bmesh(cdt_state);
+  }
+  else if (output_type == CDT_INSIDE) {
+    remove_outer_edges_until_constraints(cdt_state);
+  }
+}
+
+template<typename T>
+CDT_result<T> get_cdt_output(CDT_state<T> *cdt_state,
+                             const CDT_input<T> UNUSED(input),
+                             CDT_output_type output_type)
+{
+  prepare_cdt_for_output(cdt_state, output_type);
+  CDT_result<T> result;
+  CDTArrangement<T> *cdt = &cdt_state->cdt;
+  result.face_edge_offset = cdt_state->face_edge_offset;
+
+  /* All verts without a merge_to_index will be output.
+   * vert_to_output_map[i] will hold the output vertex index
+   * corresponding to the vert in position i in cdt->verts.
+   * This first loop sets vert_to_output_map for un-merged verts. */
+  int verts_size = cdt->verts.size();
+  Array<int> vert_to_output_map(verts_size);
+  int nv = 0;
+  for (int i = 0; i < verts_size; ++i) {
+    CDTVert<T> *v = cdt->verts[i];
+    if (v->merge_to_index == -1) {
+      vert_to_output_map[i] = nv;
+      ++nv;
+    }
+  }
+  if (nv <= 0) {
+    return result;
+  }
+  /* Now we can set vert_to_output_map for merged verts,
+   * and also add the input indices of merged verts to the input_ids
+   * list of the merge target if they were an original input id. */
+  if (nv < verts_size) {
+    for (int i = 0; i < verts_size; ++i) {
+      CDTVert<T> *v = cdt->verts[i];
+      if (v->merge_to_index != -1) {
+        if (i < cdt_state->input_vert_tot) {
+          add_to_input_ids(&cdt->verts[v->merge_to_index]->input_ids, i);
+        }
+        vert_to_output_map[i] = vert_to_output_map[v->merge_to_index];
+      }
+    }
+  }
+  result.vert = Array<vec2<T>>(nv);
+  result.vert_orig = Array<Vector<int>>(nv);
+  int i_out = 0;
+  for (int i = 0; i < verts_size; ++i) {
+    CDTVert<T> *v = cdt->verts[i];
+    if (v->merge_to_index == -1) {
+      result.vert[i_out] = v->co;
+      if (i < cdt_state->input_vert_tot) {
+        result.vert_orig[i_out].append(i);
+      }
+      for (LinkNode *ln = v->input_ids; ln; ln = ln->next) {
+        result.vert_orig[i_out].append(POINTER_AS_INT(ln->link));
+      }
+      ++i_out;
+    }
+  }
+
+  /* All non-deleted edges will be output. */
+  int ne = std::count_if(cdt->edges.begin(), cdt->edges.end(), [](const CDTEdge<T> *e) -> bool {
+    return !is_deleted_edge(e);
+  });
+  result.edge = Array<std::pair<int, int>>(ne);
+  result.edge_orig = Array<Vector<int>>(ne);
+  int e_out = 0;
+  for (const CDTEdge<T> *e : cdt->edges) {
+    if (!is_deleted_edge(e)) {
+      int vo1 = vert_to_output_map[e->symedges[0].vert->index];
+      int vo2 = vert_to_output_map[e->symedges[1].vert->index];
+      result.edge[e_out] = std::pair<int, int>(vo1, vo2);
+      for (LinkNode *ln = e->input_ids; ln; ln = ln->next) {
+        result.edge_orig[e_out].append(POINTER_AS_INT(ln->link));
+      }
+      ++e_out;
+    }
+  }
+
+  /* All non-deleted, non-outer faces will be output. */
+  int nf = std::count_if(cdt->faces.begin(), cdt->faces.end(), [=](const CDTFace<T> *f) -> bool {
+    return !f->deleted && f != cdt->outer_face;
+  });
+  result.face = Array<Vector<int>>(nf);
+  result.face_orig = Array<Vector<int>>(nf);
+  int f_out = 0;
+  for (const CDTFace<T> *f : cdt->faces) {
+    if (!f->deleted && f != cdt->outer_face) {
+      SymEdge<T> *se = f->symedge;
+      BLI_assert(se != nullptr);
+      SymEdge<T> *se_start = se;
+      do {
+        result.face[f_out].append(vert_to_output_map[se->vert->index]);
+        se = se->next;
+      } while (se != se_start);
+      for (LinkNode *ln = f->input_ids; ln; ln = ln->next) {
+        result.face_orig[f_out].append(POINTER_AS_INT(ln->link));
+      }
+      ++f_out;
+    }
+  }
+  return result;
+}
+
+/**
+ * Add all the input verts into cdt. This will deduplicate,
+ * setting vertices merge_to_index to show merges.
+ */
+template<typename T> void add_input_verts(CDT_state<T> *cdt_state, const CDT_input<T> &input)
+{
+  for (int i = 0; i < cdt_state->input_vert_tot; ++i) {
+    cdt_state->cdt.add_vert(input.vert[i]);
+  }
+}
+
+template<typename T>
+CDT_result<T> delaunay_calc(const CDT_input<T> &input, CDT_output_type output_type)
+{
+  int nv = input.vert.size();
+  int ne = input.edge.size();
+  int nf = input.face.size();
+  CDT_state<T> cdt_state(nv, ne, nf, input.epsilon);
+  add_input_verts(&cdt_state, input);
+  initial_triangulation(&cdt_state.cdt);
+  add_edge_constraints(&cdt_state, input);
+  add_face_constraints(&cdt_state, input);
+  return get_cdt_output(&cdt_state, input, output_type);
+}
+
+blender::meshintersect::CDT_result<double> delaunay_2d_calc(const CDT_input<double> &input,
+                                                            CDT_output_type output_type)
+{
+  return delaunay_calc(input, output_type);
+}
+
+#ifdef WITH_GMP
+blender::meshintersect::CDT_result<mpq_class> delaunay_2d_calc(const CDT_input<mpq_class> &input,
+                                                               CDT_output_type output_type)
+{
+  return delaunay_calc(input, output_type);
+}
+#endif
+
+} /* namespace blender::meshintersect */
+
+/* C interface. */
+
+/**
+   This function uses the double version of #CDT::delaunay_calc.
+ * Almost all of the work here is to convert between C++ #Arrays<Vector<int>>
+ * and a C version that linearizes all the elements and uses a "start"
+ * and "len" array to say where the individual vectors start and how
+ * long they are.
+ */
+extern "C" ::CDT_result *BLI_delaunay_2d_cdt_calc(const ::CDT_input *input,
+                                                  const CDT_output_type output_type)
+{
+  blender::meshintersect::CDT_input<double> in;
+  in.vert = blender::Array<blender::meshintersect::vec2<double>>(input->verts_len);
+  in.edge = blender::Array<std::pair<int, int>>(input->edges_len);
+  in.face = blender::Array<blender::Vector<int>>(input->faces_len);
+  for (int v = 0; v < input->verts_len; ++v) {
+    double x = static_cast<double>(input->vert_coords[v][0]);
+    double y = static_cast<double>(input->vert_coords[v][1]);
+    in.vert[v] = blender::meshintersect::vec2<double>(x, y);
+  }
+  for (int e = 0; e < input->edges_len; ++e) {
+    in.edge[e] = std::pair<int, int>(input->edges[e][0], input->edges[e][1]);
+  }
+  for (int f = 0; f < input->faces_len; ++f) {
+    in.face[f] = blender::Vector<int>(input->faces_len_table[f]);
+    int fstart = input->faces_start_table[f];
+    for (int j = 0; j < input->faces_len_table[f]; ++j) {
+      in.face[f][j] = input->faces[fstart + j];
+    }
+  }
+  in.epsilon = static_cast<double>(input->epsilon);
+
+  blender::meshintersect::CDT_result<double> res = blender::meshintersect::delaunay_2d_calc(
+      in, output_type);
+
+  ::CDT_result *output = static_cast<::CDT_result *>(MEM_mallocN(sizeof(*output), __func__));
+  int nv = output->verts_len = res.vert.size();
+  int ne = output->edges_len = res.edge.size();
+  int nf = output->faces_len = res.face.size();
+  int tot_v_orig = 0;
+  int tot_e_orig = 0;
+  int tot_f_orig = 0;
+  int tot_f_lens = 0;
+  for (int v = 0; v < nv; ++v) {
+    tot_v_orig += res.vert_orig[v].size();
+  }
+  for (int e = 0; e < ne; ++e) {
+    tot_e_orig += res.edge_orig[e].size();
+  }
+  for (int f = 0; f < nf; ++f) {
+    tot_f_orig += res.face_orig[f].size();
+    tot_f_lens += res.face[f].size();
+  }
+
+  output->vert_coords = static_cast<decltype(output->vert_coords)>(
+      MEM_malloc_arrayN(nv, sizeof(output->vert_coords[0]), __func__));
+  output->verts_orig = static_cast<int *>(MEM_malloc_arrayN(tot_v_orig, sizeof(int), __func__));
+  output->verts_orig_start_table = static_cast<int *>(
+      MEM_malloc_arrayN(nv, sizeof(int), __func__));
+  output->verts_orig_len_table = static_cast<int *>(MEM_malloc_arrayN(nv, sizeof(int), __func__));
+  output->edges = static_cast<decltype(output->edges)>(
+      MEM_malloc_arrayN(ne, sizeof(output->edges[0]), __func__));
+  output->edges_orig = static_cast<int *>(MEM_malloc_arrayN(tot_e_orig, sizeof(int), __func__));
+  output->edges_orig_start_table = static_cast<int *>(
+      MEM_malloc_arrayN(ne, sizeof(int), __func__));
+  output->edges_orig_len_table = static_cast<int *>(MEM_malloc_arrayN(ne, sizeof(int), __func__));
+  output->faces = static_cast<int *>(MEM_malloc_arrayN(tot_f_lens, sizeof(int), __func__));
+  output->faces_start_table = static_cast<int *>(MEM_malloc_arrayN(nf, sizeof(int), __func__));
+  output->faces_len_table = static_cast<int *>(MEM_malloc_arrayN(nf, sizeof(int), __func__));
+  output->faces_orig = static_cast<int *>(MEM_malloc_arrayN(tot_f_orig, sizeof(int), __func__));
+  output->faces_orig_start_table = static_cast<int *>(
+      MEM_malloc_arrayN(nf, sizeof(int), __func__));
+  output->faces_orig_len_table = static_cast<int *>(MEM_malloc_arrayN(nf, sizeof(int), __func__));
+
+  int v_orig_index = 0;
+  for (int v = 0; v < nv; ++v) {
+    output->vert_coords[v][0] = static_cast<float>(res.vert[v][0]);
+    output->vert_coords[v][1] = static_cast<float>(res.vert[v][1]);
+    int this_start = v_orig_index;
+    output->verts_orig_start_table[v] = this_start;
+    for (int j : res.vert_orig[v].index_range()) {
+      output->verts_orig[v_orig_index++] = res.vert_orig[v][j];
+    }
+    output->verts_orig_len_table[v] = v_orig_index - this_start;
+  }
+  int e_orig_index = 0;
+  for (int e = 0; e < ne; ++e) {
+    output->edges[e][0] = res.edge[e].first;
+    output->edges[e][1] = res.edge[e].second;
+    int this_start = e_orig_index;
+    output->edges_orig_start_table[e] = this_start;
+    for (int j : res.edge_orig[e].index_range()) {
+      output->edges_orig[e_orig_index++] = res.edge_orig[e][j];
+    }
+    output->edges_orig_len_table[e] = e_orig_index - this_start;
+  }
+  int f_orig_index = 0;
+  int f_index = 0;
+  for (int f = 0; f < nf; ++f) {
+    output->faces_start_table[f] = f_index;
+    int flen = res.face[f].size();
+    output->faces_len_table[f] = flen;
+    for (int j = 0; j < flen; ++j) {
+      output->faces[f_index++] = res.face[f][j];
+    }
+    int this_start = f_orig_index;
+    output->faces_orig_start_table[f] = this_start;
+    for (int k : res.face_orig[f].index_range()) {
+      output->faces_orig[f_orig_index++] = res.face_orig[f][k];
+    }
+    output->faces_orig_len_table[f] = f_orig_index - this_start;
+  }
+  return output;
+}
+
+extern "C" void BLI_delaunay_2d_cdt_free(::CDT_result *result)
+{
+  MEM_freeN(result->vert_coords);
+  MEM_freeN(result->edges);
+  MEM_freeN(result->faces);
+  MEM_freeN(result->faces_start_table);
+  MEM_freeN(result->faces_len_table);
+  MEM_freeN(result->verts_orig);
+  MEM_freeN(result->verts_orig_start_table);
+  MEM_freeN(result->verts_orig_len_table);
+  MEM_freeN(result->edges_orig);
+  MEM_freeN(result->edges_orig_start_table);
+  MEM_freeN(result->edges_orig_len_table);
+  MEM_freeN(result->faces_orig);
+  MEM_freeN(result->faces_orig_start_table);
+  MEM_freeN(result->faces_orig_len_table);
+  MEM_freeN(result);
+}