delaunay_2d.cc

Go to the documentation of this file.

/*
 * 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.
 */
 
#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;
}
 
template<typename Arith_t> struct CDTVert;
template<typename Arith_t> struct CDTEdge;
template<typename Arith_t> struct CDTFace;
 
template<typename Arith_t> struct SymEdge {
  SymEdge<Arith_t> *next{nullptr};
  SymEdge<Arith_t> *rot{nullptr};
  CDTVert<Arith_t> *vert{nullptr};
  CDTEdge<Arith_t> *edge{nullptr};
  CDTFace<Arith_t> *face{nullptr};
 
  SymEdge() = default;
};
 
template<typename T> inline SymEdge<T> *sym(const SymEdge<T> *se)
{
  return se->next->rot;
}
 
template<typename T> inline SymEdge<T> *prev(const SymEdge<T> *se)
{
  return se->rot->next->rot;
}
 
template<typename T> struct FatCo {
  vec2<T> exact;
  vec2<double> approx;
  vec2<double> abs_approx;
 
  FatCo();
#ifdef WITH_GMP
  FatCo(const vec2<mpq_class> &v);
#endif
  FatCo(const vec2<double> &v);
};
 
#ifdef WITH_GMP
template<> struct FatCo<mpq_class> {
  vec2<mpq_class> exact;
  vec2<double> approx;
  vec2<double> abs_approx;
 
  FatCo()
      : exact(vec2<mpq_class>(0, 0)), approx(vec2<double>(0, 0)), abs_approx(vec2<double>(0, 0))
  {
  }
 
  FatCo(const vec2<mpq_class> &v)
  {
    exact = v;
    approx = vec2<double>(v.x.get_d(), v.y.get_d());
    abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
  }
 
  FatCo(const vec2<double> &v)
  {
    exact = vec2<mpq_class>(v.x, v.y);
    approx = v;
    abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
  }
};
#endif
 
template<> struct FatCo<double> {
  vec2<double> exact;
  vec2<double> approx;
  vec2<double> abs_approx;
 
  FatCo() : exact(vec2<double>(0, 0)), approx(vec2<double>(0, 0)), abs_approx(vec2<double>(0, 0))
  {
  }
 
#ifdef WITH_GMP
  FatCo(const vec2<mpq_class> &v)
  {
    exact = vec2<double>(v.x.get_d(), v.y.get_d());
    approx = exact;
    abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
  }
#endif
 
  FatCo(const vec2<double> &v)
  {
    exact = v;
    approx = v;
    abs_approx = vec2<double>(fabs(approx.x), fabs(approx.y));
  }
};
 
template<typename T> std::ostream &operator<<(std::ostream &stream, const FatCo<T> &co)
{
  stream << "(" << co.approx.x << ", " << co.approx.y << ")";
  return stream;
}
 
template<typename T> struct CDTVert {
  FatCo<T> co;
  SymEdge<T> *symedge{nullptr};
  LinkNode *input_ids{nullptr};
  int index{-1};
  int merge_to_index{-1};
  int visit_index{0};
 
  CDTVert() = default;
  explicit CDTVert(const vec2<T> &pt);
};
 
template<typename Arith_t> struct CDTEdge {
  LinkNode *input_ids{nullptr};
  SymEdge<Arith_t> symedges[2]{SymEdge<Arith_t>(), SymEdge<Arith_t>()};
 
  CDTEdge() = default;
};
 
template<typename Arith_t> struct CDTFace {
  SymEdge<Arith_t> *symedge{nullptr};
  LinkNode *input_ids{nullptr};
  int visit_index{0};
  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. */
 
  Vector<CDTVert<Arith_t> *> verts;
  Vector<CDTEdge<Arith_t> *> edges;
  Vector<CDTFace<Arith_t> *> faces;
  CDTFace<Arith_t> *outer_face{nullptr};
 
  CDTArrangement() = default;
  ~CDTArrangement();
 
  void reserve(int num_verts, int num_edges, int num_faces);
 
  CDTVert<Arith_t> *add_vert(const vec2<Arith_t> &pt);
 
  CDTEdge<Arith_t> *add_edge(CDTVert<Arith_t> *v1,
                             CDTVert<Arith_t> *v2,
                             CDTFace<Arith_t> *fleft,
                             CDTFace<Arith_t> *fright);
 
  CDTFace<Arith_t> *add_face();
 
  CDTEdge<Arith_t> *add_vert_to_symedge_edge(CDTVert<Arith_t> *v, SymEdge<Arith_t> *se);
 
  CDTEdge<Arith_t> *add_diagonal(SymEdge<Arith_t> *s1, SymEdge<Arith_t> *s2);
 
  CDTEdge<Arith_t> *connect_separate_parts(SymEdge<Arith_t> *se1, SymEdge<Arith_t> *se2);
 
  CDTEdge<Arith_t> *split_edge(SymEdge<Arith_t> *se, Arith_t lambda);
 
  void delete_edge(SymEdge<Arith_t> *se);
 
  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;
  int input_vert_tot;
  int visit_count;
  int face_edge_offset;
  T epsilon;
 
  explicit CDT_state(int num_input_verts, int num_input_edges, int num_input_faces, T epsilon);
};
 
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 = v->co.approx[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) (((x)-minx) * scale)
#  define SY(y) ((maxy - (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<double> &uco = u->co.approx;
    const vec2<double> &vco = v->co.approx;
    int strokew = e->input_ids == nullptr ? thin_line : thick_line;
    f << R"(<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)
        << R"(" 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 << R"(<circle fill="black" cx=")" << SX(v->co.approx[0]) << "\" cy=\"" << SY(v->co.approx[1])
      << "\" r=\"" << vert_radius << "\">\n";
    f << "  <title>[" << i << "]" << v->co.approx << "</title>\n";
    f << "</circle>\n";
    if (draw_vert_labels) {
      f << "<text x=\"" << SX(v->co.approx[0]) + vert_radius << "\" y=\""
        << SY(v->co.approx[1]) - vert_radius << R"(" font-size="small">[)" << i << "]</text>\n";
    }
    ++i;
  }
 
  append = true;
#  undef SX
#  undef SY
}
#endif
 
template<typename T>
static int filtered_orient2d(const FatCo<T> &a, const FatCo<T> &b, const FatCo<T> &c);
 
#ifdef WITH_GMP
template<>
int filtered_orient2d<mpq_class>(const FatCo<mpq_class> &a,
                                 const FatCo<mpq_class> &b,
                                 const FatCo<mpq_class> &c)
{
  double det = (a.approx[0] - c.approx[0]) * (b.approx[1] - c.approx[1]) -
               (a.approx[1] - c.approx[1]) * (b.approx[0] - c.approx[0]);
  double supremum = (a.abs_approx[0] + c.abs_approx[0]) * (b.abs_approx[1] + c.abs_approx[1]) +
                    (a.abs_approx[1] + c.abs_approx[1]) * (b.abs_approx[0] + c.abs_approx[0]);
  constexpr double index_orient2d = 6;
  double err_bound = supremum * index_orient2d * DBL_EPSILON;
  if (fabs(det) > err_bound) {
    return det > 0 ? 1 : -1;
  }
  return orient2d(a.exact, b.exact, c.exact);
}
#endif
 
template<>
int filtered_orient2d<double>(const FatCo<double> &a,
                              const FatCo<double> &b,
                              const FatCo<double> &c)
{
  return orient2d(a.approx, b.approx, c.approx);
}
 
template<typename T>
static int filtered_incircle(const FatCo<T> &a,
                             const FatCo<T> &b,
                             const FatCo<T> &c,
                             const FatCo<T> &d);
 
#ifdef WITH_GMP
template<>
int filtered_incircle<mpq_class>(const FatCo<mpq_class> &a,
                                 const FatCo<mpq_class> &b,
                                 const FatCo<mpq_class> &c,
                                 const FatCo<mpq_class> &d)
{
  double adx = a.approx[0] - d.approx[0];
  double bdx = b.approx[0] - d.approx[0];
  double cdx = c.approx[0] - d.approx[0];
  double ady = a.approx[1] - d.approx[1];
  double bdy = b.approx[1] - d.approx[1];
  double cdy = c.approx[1] - d.approx[1];
  double bdxcdy = bdx * cdy;
  double cdxbdy = cdx * bdy;
  double alift = adx * adx + ady * ady;
  double cdxady = cdx * ady;
  double adxcdy = adx * cdy;
  double blift = bdx * bdx + bdy * bdy;
  double adxbdy = adx * bdy;
  double bdxady = bdx * ady;
  double clift = cdx * cdx + cdy * cdy;
  double det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);
 
  double sup_adx = a.abs_approx[0] + d.abs_approx[0]; /* index 2. */
  double sup_bdx = b.abs_approx[0] + d.abs_approx[0];
  double sup_cdx = c.abs_approx[0] + d.abs_approx[0];
  double sup_ady = a.abs_approx[1] + d.abs_approx[1];
  double sup_bdy = b.abs_approx[1] + d.abs_approx[1];
  double sup_cdy = c.abs_approx[1] + d.abs_approx[1];
  double sup_bdxcdy = sup_bdx * sup_cdy; /* index 5. */
  double sup_cdxbdy = sup_cdx * sup_bdy;
  double sup_alift = sup_adx * sup_adx + sup_ady * sup_ady; /* index 6. */
  double sup_cdxady = sup_cdx * sup_ady;
  double sup_adxcdy = sup_adx * sup_cdy;
  double sup_blift = sup_bdx * sup_bdx + sup_bdy * sup_bdy;
  double sup_adxbdy = sup_adx * sup_bdy;
  double sup_bdxady = sup_bdx * sup_ady;
  double sup_clift = sup_cdx * sup_cdx + sup_cdy * sup_cdy;
  double sup_det = sup_alift * (sup_bdxcdy + sup_cdxbdy) + sup_blift * (sup_cdxady + sup_adxcdy) +
                   sup_clift * (sup_adxbdy + sup_bdxady);
  int index = 14;
  double err_bound = sup_det * index * DBL_EPSILON;
  if (fabs(det) > err_bound) {
    return det < 0.0 ? -1 : (det > 0.0 ? 1 : 0);
  }
  return incircle(a.exact, b.exact, c.exact, d.exact);
}
#endif
 
template<>
int filtered_incircle<double>(const FatCo<double> &a,
                              const FatCo<double> &b,
                              const FatCo<double> &c,
                              const FatCo<double> &d)
{
  return incircle(a.approx, b.approx, c.approx, d.approx);
}
 
template<typename T> static bool in_line(const FatCo<T> &a, const FatCo<T> &b, const FatCo<T> &c);
 
#ifdef WITH_GMP
template<>
bool in_line<mpq_class>(const FatCo<mpq_class> &a,
                        const FatCo<mpq_class> &b,
                        const FatCo<mpq_class> &c)
{
  vec2<double> ab = b.approx - a.approx;
  vec2<double> bc = c.approx - b.approx;
  vec2<double> ac = c.approx - a.approx;
  vec2<double> supremum_ab = b.abs_approx + a.abs_approx;
  vec2<double> supremum_bc = c.abs_approx + b.abs_approx;
  vec2<double> supremum_ac = c.abs_approx + a.abs_approx;
  double dot_ab_ac = ab.x * ac.x + ab.y * ac.y;
  double supremum_dot_ab_ac = supremum_ab.x * supremum_ac.x + supremum_ab.y * supremum_ac.y;
  constexpr double index = 6;
  double err_bound = supremum_dot_ab_ac * index * DBL_EPSILON;
  if (dot_ab_ac < -err_bound) {
    return false;
  }
  double dot_bc_ac = bc.x * ac.x + bc.y * ac.y;
  double supremum_dot_bc_ac = supremum_bc.x * supremum_ac.x + supremum_bc.y * supremum_ac.y;
  err_bound = supremum_dot_bc_ac * index * DBL_EPSILON;
  if (dot_bc_ac < -err_bound) {
    return false;
  }
  vec2<mpq_class> exact_ab = b.exact - a.exact;
  vec2<mpq_class> exact_ac = c.exact - a.exact;
  if (vec2<mpq_class>::dot(exact_ab, exact_ac) < 0) {
    return false;
  }
  vec2<mpq_class> exact_bc = c.exact - b.exact;
  return vec2<mpq_class>::dot(exact_bc, exact_ac) >= 0;
}
#endif
 
template<>
bool in_line<double>(const FatCo<double> &a, const FatCo<double> &b, const FatCo<double> &c)
{
  vec2<double> ab = b.approx - a.approx;
  vec2<double> ac = c.approx - a.approx;
  if (vec2<double>::dot(ab, ac) < 0) {
    return false;
  }
  vec2<double> bc = c.approx - b.approx;
  return vec2<double>::dot(bc, ac) >= 0;
}
 
template<> CDTVert<double>::CDTVert(const vec2<double> &pt)
{
  this->co.exact = pt;
  this->co.approx = pt;
  this->co.abs_approx = pt; /* Not used, so doesn't matter. */
  this->input_ids = nullptr;
  this->symedge = nullptr;
  this->index = -1;
  this->merge_to_index = -1;
  this->visit_index = 0;
}
 
#ifdef WITH_GMP
template<> CDTVert<mpq_class>::CDTVert(const vec2<mpq_class> &pt)
{
  this->co.exact = pt;
  this->co.approx = double2(pt.x.get_d(), pt.y.get_d());
  this->co.abs_approx = double2(fabs(this->co.approx.x), fabs(this->co.approx.y));
  this->input_ids = nullptr;
  this->symedge = nullptr;
  this->index = -1;
  this->merge_to_index = -1;
  this->visit_index = 0;
}
#endif
 
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 != nullptr;
}
 
template<typename T> inline bool is_deleted_edge(const CDTEdge<T> *e)
{
  return e->symedges[0].next == nullptr;
}
 
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;
}
 
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;
}
 
template<typename T> inline bool exists_edge(const CDTVert<T> *v1, const CDTVert<T> *v2)
{
  return find_symedge_between_verts(v1, v2) != nullptr;
}
 
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;
}
 
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;
}
 
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;
}
 
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.exact;
  const vec2<T> *b = &se->next->vert->co.exact;
  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;
}
 
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;
};
 
template<typename T> bool site_lexicographic_sort(const SiteInfo<T> &a, const SiteInfo<T> &b)
{
  const vec2<T> &co_a = a.v->co.exact;
  const vec2<T> &co_b = b.v->co.exact;
  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;
}
 
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.exact == sites[i].v->co.exact) {
      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 filtered_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 filtered_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 filtered_orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0;
}
 
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 = filtered_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 separately, 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 (filtered_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 (filtered_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 && filtered_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 eliminated verts that merge to others. */
  int i = 0;
  int j = 0;
  int nsites = sites.size();
  while (j < nsites) {
    /* Invariant: `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);
}
 
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);
}
 
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 (filtered_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 filtered_orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co);
}
 
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);
 
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 = nullptr;
  cd_next->out = nullptr;
  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;
  }
}
 
template<typename T>
void fill_crossdata_for_intersect(const FatCo<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.exact, vb->co.exact, curco.exact, v2->co.exact);
  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
        double len_ab = vec2<double>::distance(va->co.approx, vb->co.approx);
        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<double>::distance_squared(va->co.approx, v2->co.approx) <=
          vec2<double>::distance_squared(vb->co.approx, v2->co.approx)) {
        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
 
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 = filtered_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 = filtered_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;
}
 
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.exact, vb->co.exact, cd->lambda);
  FatCo<T> fat_curco(curco);
  SymEdge<T> *se_ac = sym(cd->in)->next;
  CDTVert<T> *vc = se_ac->next->vert;
  int orient = filtered_orient2d(fat_curco, v2->co, vc->co);
  if (orient < 0) {
    fill_crossdata_for_intersect<T>(fat_curco, v2, se_ac->next, cd, cd_next, epsilon);
  }
  else if (orient > 0.0) {
    fill_crossdata_for_intersect(fat_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";
  }
}
 
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 = {nullptr, nullptr};
 
  if (r_edges) {
    *r_edges = nullptr;
  }
 
  /*
   * 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 != nullptr) {
      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<T>(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 == nullptr) {
            return;
          }
          cd_prev->out = se;
          cd->in = nullptr;
        }
        else {
          se = find_symedge_with_face(v, sym(cd_prev->in)->face);
          if (se == nullptr) {
            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 = nullptr;
    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 == nullptr) {
        /* 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 != nullptr) {
          BLI_linklist_append(&edge_list, edge);
        }
      }
    }
    if (t != nullptr) {
      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 != nullptr) {
        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 != nullptr) {
        tstart = tnext;
      }
    }
  }
 
  if (r_edges) {
    *r_edges = edge_list.list;
  }
}
 
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;
}
 
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;
}
 
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 (ELEM(cdt->outer_face->symedge, se, symse)) {
    /* Advancing by 2 to get past possible 'sym(se)'. */
    if (se->next->next == se) {
      cdt->outer_face->symedge = nullptr;
    }
    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);
}
 
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.
 */
 
template<typename T> struct EdgeToSort {
  double len_squared = 0.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<double> &co1 = e->symedges[0].vert->co.approx;
      const vec2<double> &co2 = e->symedges[1].vert->co.approx;
      dissolvable_edges[i].len_squared = vec2<double>::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 != nullptr) {
      se = static_cast<SymEdge<T> *>(BLI_linklist_pop(&to_dissolve));
      if (se->next != nullptr) {
        dissolve_symedge(cdt_state, se);
      }
    }
  }
}
 
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.exact;
      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;
}
 
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. */
 
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);
}