path: root/source/blender/blenlib/intern/delaunay_2d.c

/*
 * 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
 *
 * Constrained 2d Delaunay Triangulation.
 */

#include "MEM_guardedalloc.h"

#include "BLI_array.h"
#include "BLI_bitmap.h"
#include "BLI_linklist.h"
#include "BLI_math.h"
#include "BLI_memarena.h"
#include "BLI_mempool.h"

#include "BLI_delaunay_2d.h"

/* Uncomment this define to get helpful debugging functions etc. defined. */
// #define DEBUG_CDT

struct CDTEdge;
struct CDTFace;
struct CDTVert;

typedef struct SymEdge {
  struct SymEdge *next; /* In face, doing CCW traversal of face. */
  struct SymEdge *rot;  /* CCW around vert. */
  struct CDTVert *vert; /* Vert at origin. */
  struct CDTEdge *edge; /* Undirected edge this is for. */
  struct CDTFace *face; /* Face on left side. */
} SymEdge;

typedef struct CDTVert {
  double co[2];        /* Coordinate. */
  SymEdge *symedge;    /* Some edge attached to it. */
  LinkNode *input_ids; /* List of corresponding vertex input ids. */
  int index;           /* Index into array that cdt keeps. */
  int merge_to_index;  /* Index of a CDTVert that this has merged to. -1 if no merge. */
  int visit_index;     /* Which visit epoch has this been seen. */
} CDTVert;

typedef struct CDTEdge {
  LinkNode *input_ids; /* List of input edge ids that this is part of. */
  SymEdge symedges[2]; /* The directed edges for this edge. */
  bool in_queue;       /* Used in flipping algorithm. */
} CDTEdge;

typedef struct CDTFace {
  SymEdge *symedge;    /* A symedge in face; only used during output, so only valid then. */
  LinkNode *input_ids; /* List of input face ids that this is part of. */
  int visit_index;     /* Which visit epoch has this been seen. */
  bool deleted;        /* Marks this face no longer used. */
  bool in_queue;       /* Used in remove_small_features algorithm. */
} CDTFace;

typedef struct CDT_state {
  LinkNode *edges;          /* List of CDTEdge pointer. */
  LinkNode *faces;          /* List of CDTFace pointer. */
  CDTFace *outer_face;      /* Which CDTFace is the outer face. */
  CDTVert **vert_array;     /* Array of CDTVert pointer, grows. */
  int vert_array_len;       /* Current length of vert_array. */
  int vert_array_len_alloc; /* Allocated length of vert_array. */
  int input_vert_tot;       /* How many verts were in input (will be first in vert_array). */
  double minx;              /* Used for debug drawing. */
  double miny;              /* Used for debug drawing. */
  double maxx;              /* Used for debug drawing. */
  double maxy;              /* Used for debug drawing. */
  double margin;            /* Used for debug drawing. */
  int visit_count; /* Used for visiting things without having to initialized their visit fields. */
  int face_edge_offset;  /* Input edge id where we start numbering the face edges. */
  MemArena *arena;       /* Most allocations are done from here, so can free all at once at end. */
  BLI_mempool *listpool; /* Allocations of ListNodes done from this pool. */
  double epsilon;        /* The user-specified nearness limit. */
  double epsilon_squared; /* Square of epsilon. */
  bool output_prepared;   /* Set after the mesh has been modified for output (may not be all
                             triangles now). */
} CDT_state;

#define DLNY_ARENASIZE 1 << 14

#ifdef DEBUG_CDT
#  ifdef __GNUC__
#    define ATTU __attribute__((unused))
#  else
#    define ATTU
#  endif
#  define F2(p) p[0], p[1]
#  define F3(p) p[0], p[1], p[2]
struct CrossData;
ATTU static void dump_se(const SymEdge *se, const char *lab);
ATTU static void dump_se_short(const SymEdge *se, const char *lab);
ATTU static void dump_v(const CDTVert *v, const char *lab);
ATTU static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit);
ATTU static void dump_id_list(const LinkNode *id_list, const char *lab);
ATTU static void dump_cross_data(struct CrossData *cd, const char *lab);
ATTU static void dump_cdt(const CDT_state *cdt, const char *lab);
ATTU static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab);
ATTU static void cdt_draw(CDT_state *cdt, const char *lab);
ATTU static void cdt_draw_region(
    CDT_state *cdt, const char *lab, double minx, double miny, double maxx, double maxy);

ATTU static void cdt_draw_vertex_region(CDT_state *cdt, int v, double dist, const char *lab);
ATTU static void cdt_draw_edge_region(
    CDT_state *cdt, int v1, int v2, double dist, const char *lab);
ATTU static void write_cdt_input_to_file(const CDT_input *inp);
ATTU static void validate_cdt(CDT_state *cdt,
                              bool check_all_tris,
                              bool check_delaunay,
                              bool check_visibility);
#endif

static void exactinit(void);
static double orient2d(const double *pa, const double *pb, const double *pc);
static double incircle(const double *pa, const double *pb, const double *pc, const double *pd);

/** Return other #SymEdge for same #CDTEdge as se. */
BLI_INLINE SymEdge *sym(const SymEdge *se)
{
  return se->next->rot;
}

/** Return SymEdge whose next is se. */
BLI_INLINE SymEdge *prev(const SymEdge *se)
{
  return se->rot->next->rot;
}

/**
 * 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.  */
static bool in_line(const double a[2], const double b[2], const double c[2])
{
  double ab[2], bc[2], ac[2];
  sub_v2_v2v2_db(ab, b, a);
  sub_v2_v2v2_db(bc, c, b);
  sub_v2_v2v2_db(ac, c, a);
  if (dot_v2v2_db(ab, ac) < 0.0) {
    return false;
  }
  return dot_v2v2_db(bc, ac) >= 0.0;
}

#ifndef NDEBUG
/** Is s2 reachable from s1 by next pointers with < limit hops? */
static bool reachable(SymEdge *s1, SymEdge *s2, int limit)
{
  int count = 0;
  for (SymEdge *s = s1; s && count < limit; s = s->next) {
    if (s == s2) {
      return true;
    }
    count++;
  }
  return false;
}
#endif

/** Using array to store these instead of linked list so can make a random selection from them. */
static CDTVert *add_cdtvert(CDT_state *cdt, double x, double y)
{
  CDTVert *v = BLI_memarena_alloc(cdt->arena, sizeof(*v));
  v->co[0] = x;
  v->co[1] = y;
  v->input_ids = NULL;
  v->symedge = NULL;
  if (cdt->vert_array_len == cdt->vert_array_len_alloc) {
    CDTVert **old_array = cdt->vert_array;
    cdt->vert_array_len_alloc *= 4;
    cdt->vert_array = BLI_memarena_alloc(cdt->arena,
                                         cdt->vert_array_len_alloc * sizeof(cdt->vert_array[0]));
    memmove(cdt->vert_array, old_array, cdt->vert_array_len * sizeof(cdt->vert_array[0]));
  }
  BLI_assert(cdt->vert_array_len < cdt->vert_array_len_alloc);
  v->index = cdt->vert_array_len;
  v->merge_to_index = -1;
  v->visit_index = 0;
  cdt->vert_array[cdt->vert_array_len++] = v;
  return v;
}

static CDTEdge *add_cdtedge(
    CDT_state *cdt, CDTVert *v1, CDTVert *v2, CDTFace *fleft, CDTFace *fright)
{
  CDTEdge *e = BLI_memarena_alloc(cdt->arena, sizeof(*e));
  SymEdge *se = &e->symedges[0];
  SymEdge *sesym = &e->symedges[1];
  e->input_ids = NULL;
  e->in_queue = false;
  BLI_linklist_prepend_arena(&cdt->edges, (void *)e, cdt->arena);
  se->edge = sesym->edge = e;
  se->face = fleft;
  sesym->face = fright;
  se->vert = v1;
  if (v1->symedge == NULL) {
    v1->symedge = se;
  }
  sesym->vert = v2;
  if (v2->symedge == NULL) {
    v2->symedge = sesym;
  }
  se->next = sesym->next = se->rot = sesym->rot = NULL;
  return e;
}

static CDTFace *add_cdtface(CDT_state *cdt)
{
  CDTFace *f = BLI_memarena_alloc(cdt->arena, sizeof(*f));
  f->visit_index = 0;
  f->deleted = false;
  f->symedge = NULL;
  f->input_ids = NULL;
  f->in_queue = false;
  BLI_linklist_prepend_arena(&cdt->faces, (void *)f, cdt->arena);
  return f;
}

static bool id_in_list(const LinkNode *id_list, int id)
{
  const LinkNode *ln;

  for (ln = id_list; ln; 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; 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, CDT_state *cdt)
{
  if (!id_in_list(*dst, input_id)) {
    BLI_linklist_prepend_arena(dst, POINTER_FROM_INT(input_id), cdt->arena);
  }
}

static void add_list_to_input_ids(LinkNode **dst, const LinkNode *src, CDT_state *cdt)
{
  const LinkNode *ln;

  for (ln = src; ln; ln = ln->next) {
    add_to_input_ids(dst, POINTER_AS_INT(ln->link), cdt);
  }
}

BLI_INLINE bool is_border_edge(const CDTEdge *e, const CDT_state *cdt)
{
  return e->symedges[0].face == cdt->outer_face || e->symedges[1].face == cdt->outer_face;
}

BLI_INLINE bool is_constrained_edge(const CDTEdge *e)
{
  return e->input_ids != NULL;
}

BLI_INLINE bool is_deleted_edge(const CDTEdge *e)
{
  return e->symedges[0].next == NULL;
}

BLI_INLINE bool is_original_vert(const CDTVert *v, CDT_state *cdt)
{
  return (v->index < cdt->input_vert_tot);
}

/** Return the Symedge that goes from v1 to v2, if it exists, else return NULL. */
static SymEdge *find_symedge_between_verts(const CDTVert *v1, const CDTVert *v2)
{
  SymEdge *tstart, *t;

  t = tstart = v1->symedge;
  do {
    if (t->next->vert == v2) {
      return t;
    }
  } while ((t = t->rot) != tstart);
  return NULL;
}

/** Return the SymEdge attached to v that has face f, if it exists, else return NULL. */
static SymEdge *find_symedge_with_face(const CDTVert *v, const CDTFace *f)
{
  SymEdge *tstart, *t;

  t = tstart = v->symedge;
  do {
    if (t->face == f) {
      return t;
    }
  } while ((t = t->rot) != tstart);
  return NULL;
}

/** Is there already an edge between a and b? */
static inline bool exists_edge(const CDTVert *v1, const CDTVert *v2)
{
  return find_symedge_between_verts(v1, v2) != NULL;
}

/** Is the vertex v incident on face f? */
static bool vert_touches_face(const CDTVert *v, const CDTFace *f)
{
  SymEdge *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 subface
 * that has s1, and a new face will be made for s2's new face.
 * Return the new diagonal's CDTEdge *.
 */
static CDTEdge *add_diagonal(CDT_state *cdt, SymEdge *s1, SymEdge *s2)
{
  CDTEdge *ediag;
  CDTFace *fold, *fnew;
  SymEdge *sdiag, *sdiagsym;
  SymEdge *s1prev, *s1prevsym, *s2prev, *s2prevsym, *se;
  BLI_assert(reachable(s1, s2, 20000));
  BLI_assert(reachable(s2, s1, 20000));
  fold = s1->face;
  fnew = add_cdtface(cdt);
  s1prev = prev(s1);
  s1prevsym = sym(s1prev);
  s2prev = prev(s2);
  s2prevsym = sym(s2prev);
  ediag = add_cdtedge(cdt, s1->vert, s2->vert, fnew, fold);
  sdiag = &ediag->symedges[0];
  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;
#ifdef DEBUG_CDT
  BLI_assert(reachable(s2, sdiag, 2000));
#endif
  for (se = s2; se != sdiag; se = se->next) {
    se->face = fnew;
  }
  add_list_to_input_ids(&fnew->input_ids, fold->input_ids, cdt);
  return ediag;
}

/**
 * Add a dangling edge from an isolated v to the vert at se in the same face as se->face.
 */
static CDTEdge *add_vert_to_symedge_edge(CDT_state *cdt, CDTVert *v, SymEdge *se)
{
  CDTEdge *e;
  SymEdge *se_rot, *se_rotsym, *new_se, *new_se_sym;

  se_rot = se->rot;
  se_rotsym = sym(se_rot);
  e = add_cdtedge(cdt, v, se->vert, se->face, se->face);
  new_se = &e->symedges[0];
  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 #SymEdges are on
 * the outer boundary (have face == outer_face) of two components that are isolated from
 * each other.
 */
static CDTEdge *connect_separate_parts(CDT_state *cdt, SymEdge *se1, SymEdge *se2)
{
  CDTEdge *e;
  SymEdge *se1_rot, *se1_rotsym, *se2_rot, *se2_rotsym, *new_se, *new_se_sym;

  BLI_assert(se1->face == cdt->outer_face && se2->face == cdt->outer_face);
  se1_rot = se1->rot;
  se1_rotsym = sym(se1_rot);
  se2_rot = se2->rot;
  se2_rotsym = sym(se2_rot);
  e = add_cdtedge(cdt, se1->vert, se2->vert, cdt->outer_face, cdt->outer_face);
  new_se = &e->symedges[0];
  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 \a se at fraction \a lambda,
 * and return the new #CDTEdge that is the new second half.
 * Copy the edge input_ids into the new one.
 */
static CDTEdge *split_edge(CDT_state *cdt, SymEdge *se, double lambda)
{
  const double *a, *b;
  double p[2];
  CDTVert *v;
  CDTEdge *e;
  SymEdge *sesym, *newse, *newsesym, *senext, *sesymprev, *sesymprevsym;
  /* Split e at lambda. */
  a = se->vert->co;
  b = se->next->vert->co;
  sesym = sym(se);
  sesymprev = prev(sesym);
  sesymprevsym = sym(sesymprev);
  senext = se->next;
  p[0] = (1.0 - lambda) * a[0] + lambda * b[0];
  p[1] = (1.0 - lambda) * a[1] + lambda * b[1];
  v = add_cdtvert(cdt, p[0], p[1]);
  e = add_cdtedge(cdt, v, se->next->vert, se->face, sesym->face);
  sesym->vert = v;
  newse = &e->symedges[0];
  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, cdt);
  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 #SymEdge(s) to NULL.
 * <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.
 */
static void delete_edge(CDT_state *cdt, SymEdge *e)
{
  SymEdge *esym, *f, *h, *i, *j, *k, *jsym, *hsym;
  CDTFace *aface, *bface;
  CDTVert *v1, *v2;
  bool v1_isolated, v2_isolated;

  esym = sym(e);
  v1 = e->vert;
  v2 = esym->vert;
  aface = e->face;
  bface = esym->face;
  f = e->next;
  h = prev(e);
  i = esym->next;
  j = prev(esym);
  jsym = sym(j);
  hsym = sym(h);
  v1_isolated = (i == e);
  v2_isolated = (f == esym);

  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 (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 = NULL;
  }
  else if (v1->symedge == e) {
    v1->symedge = i;
  }
  if (v2_isolated) {
    v2->symedge = NULL;
  }
  else if (v2->symedge == esym) {
    v2->symedge = f;
  }

  /* Mark SymEdge as deleted by setting all its pointers to NULL. */
  e->next = e->rot = NULL;
  esym->next = esym->rot = NULL;
  if (!v1_isolated && !v2_isolated && aface != bface) {
    bface->deleted = true;
    if (cdt->outer_face == bface) {
      cdt->outer_face = aface;
    }
  }
}

static CDT_state *cdt_init(const CDT_input *in)
{
  int i;
  MemArena *arena = BLI_memarena_new(DLNY_ARENASIZE, __func__);
  CDT_state *cdt = BLI_memarena_calloc(arena, sizeof(CDT_state));

  cdt->epsilon = (double)in->epsilon;
  cdt->epsilon_squared = cdt->epsilon * cdt->epsilon;
  cdt->arena = arena;
  cdt->input_vert_tot = in->verts_len;
  cdt->vert_array_len_alloc = 2 * in->verts_len;
  cdt->vert_array = BLI_memarena_alloc(arena,
                                       cdt->vert_array_len_alloc * sizeof(*cdt->vert_array));
  cdt->listpool = BLI_mempool_create(
      sizeof(LinkNode), 128 + 4 * in->verts_len, 128 + in->verts_len, 0);

  for (i = 0; i < in->verts_len; i++) {
    add_cdtvert(cdt, (double)(in->vert_coords[i][0]), (double)(in->vert_coords[i][1]));
  }
  cdt->outer_face = add_cdtface(cdt);
  return cdt;
}

static void new_cdt_free(CDT_state *cdt)
{
  BLI_mempool_destroy(cdt->listpool);
  BLI_memarena_free(cdt->arena);
}

typedef struct SiteInfo {
  CDTVert *v;
  int orig_index;
} SiteInfo;

static int site_lexicographic_cmp(const void *a, const void *b)
{
  const SiteInfo *s1 = a;
  const SiteInfo *s2 = b;
  const double *co1 = s1->v->co;
  const double *co2 = s2->v->co;

  if (co1[0] < co2[0]) {
    return -1;
  }
  if (co1[0] > co2[0]) {
    return 1;
  }
  if (co1[1] < co2[1]) {
    return -1;
  }
  if (co1[1] > co2[1]) {
    return 1;
  }
  if (s1->orig_index < s2->orig_index) {
    return -1;
  }
  if (s1->orig_index > s2->orig_index) {
    return 1;
  }
  return 0;
}

BLI_INLINE bool vert_left_of_symedge(CDTVert *v, SymEdge *se)
{
  return orient2d(v->co, se->vert->co, se->next->vert->co) > 0.0;
}

BLI_INLINE bool vert_right_of_symedge(CDTVert *v, SymEdge *se)
{
  return orient2d(v->co, se->next->vert->co, se->vert->co) > 0.0;
}

/* Is se above basel? */
BLI_INLINE bool dc_tri_valid(SymEdge *se, SymEdge *basel, SymEdge *basel_sym)
{
  return orient2d(se->next->vert->co, basel_sym->vert->co, basel->vert->co) > 0.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.
 */
static void dc_tri(
    CDT_state *cdt, SiteInfo *sites, int start, int end, SymEdge **r_le, SymEdge **r_re)
{
  int n = end - start;
  int n2;
  CDTVert *v1, *v2, *v3;
  CDTEdge *ea, *eb, *ebasel;
  SymEdge *ldo, *ldi, *rdi, *rdo, *basel, *basel_sym, *lcand, *rcand, *t;
  double orient;
  bool valid_lcand, valid_rcand;
#ifdef DEBUG_CDT
  char label_buf[100];
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "DC_TRI start=%d end=%d\n", start, end);
  }
#endif

  BLI_assert(r_le != NULL && r_re != NULL);
  if (n <= 1) {
    *r_le = NULL;
    *r_re = NULL;
    return;
  }
  if (n <= 3) {
    v1 = sites[start].v;
    v2 = sites[start + 1].v;
    ea = add_cdtedge(cdt, 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;
    }
    v3 = sites[start + 2].v;
    eb = add_vert_to_symedge_edge(cdt, v3, &ea->symedges[1]);
    orient = orient2d(v1->co, v2->co, v3->co);
    if (orient > 0.0) {
      add_diagonal(cdt, &eb->symedges[0], &ea->symedges[0]);
      *r_le = &ea->symedges[0];
      *r_re = &eb->symedges[0];
    }
    else if (orient < 0.0) {
      add_diagonal(cdt, &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;
  }
  /* Here: n >= 4.  Divide and conquer. */
  n2 = n / 2;
  BLI_assert(n2 >= 2 && end - (start + n2) >= 2);

  /* Delaunay triangulate two halves, L and R. */
  dc_tri(cdt, sites, start, start + n2, &ldo, &ldi);
  dc_tri(cdt, sites, start + n2, end, &rdi, &rdo);
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "\nDC_TRI merge step for start=%d, end=%d\n", start, end);
    dump_se(ldo, "ldo");
    dump_se(ldi, "ldi");
    dump_se(rdi, "rdi");
    dump_se(rdo, "rdo");
    if (dbg_level > 1) {
      sprintf(label_buf, "dc_tri(%d,%d)(%d,%d)", start, start + n2, start + n2, end);
      /* dump_cdt(cdt, label_buf); */
      cdt_draw(cdt, label_buf);
    }
  }
#endif

  /* 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;
    }
  }
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "common lower tangent is between\n");
    dump_se(rdi, "rdi");
    dump_se(ldi, "ldi");
  }
#endif
  ebasel = connect_separate_parts(cdt, sym(rdi)->next, ldi);
  basel = &ebasel->symedges[0];
  basel_sym = &ebasel->symedges[1];
#ifdef DEBUG_CDT
  if (dbg_level > 1) {
    dump_se(basel, "basel");
    cdt_draw(cdt, "after basel made");
  }
#endif
  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. */
    lcand = basel_sym->rot;
    rcand = basel_sym->next;
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "\ntop of merge loop\n");
      dump_se(lcand, "lcand");
      dump_se(rcand, "rcand");
      dump_se(basel, "basel");
    }
#endif
    if (dc_tri_valid(lcand, basel, basel_sym)) {
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "found valid lcand\n");
        dump_se(lcand, "  lcand");
      }
#endif
      while (incircle(basel_sym->vert->co,
                      basel->vert->co,
                      lcand->next->vert->co,
                      lcand->rot->next->vert->co) > 0.0) {
#ifdef DEBUG_CDT
        if (dbg_level > 1) {
          fprintf(stderr, "incircle says to remove lcand\n");
          dump_se(lcand, "  lcand");
        }
#endif
        t = lcand->rot;
        delete_edge(cdt, sym(lcand));
        lcand = t;
      }
    }
    /* Symmetrically, locate first R point to be hit and delete R edges. */
    if (dc_tri_valid(rcand, basel, basel_sym)) {
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "found valid rcand\n");
        dump_se(rcand, "  rcand");
      }
#endif
      while (incircle(basel_sym->vert->co,
                      basel->vert->co,
                      rcand->next->vert->co,
                      sym(rcand)->next->next->vert->co) > 0.0) {
#ifdef DEBUG_CDT
        if (dbg_level > 0) {
          fprintf(stderr, "incircle says to remove rcand\n");
          dump_se(lcand, "  rcand");
        }
#endif
        t = sym(rcand)->next;
        delete_edge(cdt, rcand);
        rcand = t;
      }
    }
    /* If both lcand and rcand are invalid, then basel is the common upper tangent. */
    valid_lcand = dc_tri_valid(lcand, basel, basel_sym);
    valid_rcand = dc_tri_valid(rcand, basel, basel_sym);
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(
          stderr, "after bubbling up, valid_lcand=%d, valid_rcand=%d\n", valid_lcand, valid_rcand);
      dump_se(lcand, "lcand");
      dump_se(rcand, "rcand");
    }
#endif
    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.0)) {
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr, "connecting rcand\n");
        dump_se(basel_sym, "  se1=basel_sym");
        dump_se(rcand->next, "  se2=rcand->next");
      }
#endif
      ebasel = add_diagonal(cdt, rcand->next, basel_sym);
    }
    else {
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr, "connecting lcand\n");
        dump_se(sym(lcand), "  se1=sym(lcand)");
        dump_se(basel_sym->next, "  se2=basel_sym->next");
      }
#endif
      ebasel = add_diagonal(cdt, basel_sym->next, sym(lcand));
    }
    basel = &ebasel->symedges[0];
    basel_sym = &ebasel->symedges[1];
    BLI_assert(basel_sym->face == cdt->outer_face);
#ifdef DEBUG_CDT
    if (dbg_level > 2) {
      cdt_draw(cdt, "after adding new crossedge");
      // dump_cdt(cdt, "after adding new crossedge");
    }
#endif
  }
  *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. */
static void dc_triangulate(CDT_state *cdt, SiteInfo *sites, int nsites)
{
  int i, j, n;
  SymEdge *le, *re;

  /* Compress sites in place to eliminated verts that merge to others. */
  i = 0;
  j = 0;
  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++;
    }
  }
  n = i;
  if (n == 0) {
    return;
  }
  dc_tri(cdt, sites, 0, n, &le, &re);
}

/**
 * Do a Delaunay Triangulation of the points in cdt->vert_array.
 * 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.
 * The incircle and ccw tests are done using Shewchuk's exact
 * primitives (see below), so that this routine is robust.
 *
 * As a preprocessing step, we want to merge all vertices that are
 * within cdt->epsilon of each other. 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.
 */
static void initial_triangulation(CDT_state *cdt)
{
  int i, j, n;
  SiteInfo *sites;
  double *ico, *jco;
  double xend, yend, xcur;
  double epsilon = cdt->epsilon;
  double epsilon_squared = cdt->epsilon_squared;
#ifdef SJF_WAY
  CDTEdge *e;
  CDTVert *va, *vb;
#endif
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nINITIAL TRIANGULATION\n\n");
  }
#endif

  /* First sort the vertices by lexicographic order of their
   * coordinates, breaking ties by putting earlier original-index
   * vertices first.
   */
  n = cdt->vert_array_len;
  if (n <= 1) {
    return;
  }
  sites = MEM_malloc_arrayN(n, sizeof(SiteInfo), __func__);
  for (i = 0; i < n; i++) {
    sites[i].v = cdt->vert_array[i];
    sites[i].orig_index = i;
  }
  qsort(sites, n, sizeof(SiteInfo), site_lexicographic_cmp);
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "after sorting\n");
    for (i = 0; i < n; i++) {
      fprintf(stderr, "%d: orig index: %d, (%f,%f)\n", i, sites[i].orig_index, F2(sites[i].v->co));
    }
  }
#endif

  /* Now de-duplicate according to user-defined epsilon.
   * We will merge a vertex into an earlier-indexed vertex
   * that is within epsilon (Euclidean distance).
   * Merges may cascade. So we may end up merging two things
   * that are farther than epsilon by transitive merging. Oh well.
   * Assume that merges are rare, so use simple searches in the
   * lexicographic ordering - likely we will soon hit y's with
   * the same x that are farther away than epsilon, and then
   * skipping ahead to the next biggest x, are likely to soon
   * find one of those farther away than epsilon.
   */
  for (i = 0; i < n - 1; i++) {
    ico = sites[i].v->co;
    /* Start j at next place that has both x and y coords within epsilon. */
    xend = ico[0] + epsilon;
    yend = ico[1] + epsilon;
    j = i + 1;
    while (j < n) {
      jco = sites[j].v->co;
      if (jco[0] > xend) {
        break; /* No more j's to process. */
      }
      if (jco[1] > yend) {
        /* Get past any string of v's with the same x and too-big y. */
        xcur = jco[0];
        while (++j < n) {
          if (sites[j].v->co[0] > xcur) {
            break;
          }
        }
        BLI_assert(j == n || sites[j].v->co[0] > xcur);
        if (j == n) {
          break;
        }
        jco = sites[j].v->co;
        if (jco[0] > xend || jco[1] > yend) {
          break;
        }
      }
      /* When here, vertex i and j are within epsilon by box test.
       * The Euclidean distance test is stricter, so need to do it too, now.
       */
      BLI_assert(j < n && jco[0] <= xend && jco[1] <= yend);
      if (len_squared_v2v2_db(ico, jco) <= epsilon_squared) {
        sites[j].v->merge_to_index = (sites[i].v->merge_to_index == -1) ?
                                         sites[i].orig_index :
                                         sites[i].v->merge_to_index;
#ifdef DEBUG_CDT
        if (dbg_level > 0) {
          fprintf(stderr,
                  "merged orig vert %d to %d\n",
                  sites[j].orig_index,
                  sites[j].v->merge_to_index);
        }
#endif
      }
      j++;
    }
  }

  /* Now add non-dup vertices into triangulation in lexicographic order. */

  dc_triangulate(cdt, sites, n);
  MEM_freeN(sites);
}

/**
 * Use #LinkNode linked list as stack of #SymEdges, allocating from `cdt->listpool` .
 */
typedef LinkNode *Stack;

BLI_INLINE void push(Stack *stack, SymEdge *se, CDT_state *cdt)
{
  BLI_linklist_prepend_pool(stack, se, cdt->listpool);
}

BLI_INLINE SymEdge *pop(Stack *stack, CDT_state *cdt)
{
  return (SymEdge *)BLI_linklist_pop_pool(stack, cdt->listpool);
}

BLI_INLINE bool is_empty(Stack *stack)
{
  return *stack == NULL;
}

/**
 * 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".
 */
static void re_delaunay_triangulate(CDT_state *cdt, SymEdge *se)
{
  SymEdge *ss, *first, *cse;
  CDTVert *a, *b, *c, *v;
  CDTEdge *ebc, *eca;
  int count;
#ifdef DEBUG_CDT
  SymEdge *last;
  const int dbg_level = 0;
#endif

  if (se->face == cdt->outer_face || sym(se)->face == cdt->outer_face) {
    return;
  }
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "retriangulate");
    dump_se_cycle(se, "poly ", 1000);
  }
#endif
  /* 'se' is a diagonal just added, and it is base of area to retriangulate (face on its left) */
  count = 1;
  for (ss = se->next; ss != se; ss = ss->next) {
    count++;
  }
  if (count <= 3) {
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "nothing to do\n");
    }
#endif
    return;
  }
  /* First and last are the SymEdges whose verts are first and last off of base,
   * continuing from 'se'. */
  first = se->next->next;
  /* We want to make a triangle with 'se' as base and some other c as 3rd vertex. */
  a = se->vert;
  b = se->next->vert;
  c = first->vert;
  cse = first;
#ifdef DEBUG_CDT
  last = prev(se);
  if (dbg_level > 1) {
    dump_se(first, "first");
    dump_se(last, "last");
    dump_v(a, "a");
    dump_v(b, "b");
    dump_v(c, "c");
  }
#endif
  for (ss = first->next; ss != se; ss = ss->next) {
    v = ss->vert;
    if (incircle(a->co, b->co, c->co, v->co) > 0.0) {
      c = v;
      cse = ss;
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        dump_v(c, "new c ");
      }
#endif
    }
  }
  /* Add diagonals necessary to make abc a triangle. */
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "make triangle abc exist where\n");
    dump_v(a, "  a");
    dump_v(b, "  b");
    dump_v(c, "  c");
  }
#endif
  ebc = NULL;
  eca = NULL;
  if (!exists_edge(b, c)) {
    ebc = add_diagonal(cdt, se->next, cse);
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "added edge ebc\n");
      dump_se(&ebc->symedges[0], "  ebc");
    }
#endif
  }
  if (!exists_edge(c, a)) {
    eca = add_diagonal(cdt, cse, se);
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "added edge eca\n");
      dump_se(&eca->symedges[0], "  eca");
    }
#endif
  }
  /* Now recurse. */
  if (ebc) {
    re_delaunay_triangulate(cdt, &ebc->symedges[1]);
  }
  if (eca) {
    re_delaunay_triangulate(cdt, &eca->symedges[1]);
  }
}

static double tri_orient(const SymEdge *t)
{
  return orient2d(t->vert->co, t->next->vert->co, t->next->next->vert->co);
}

/**
 * The CrossData struct gives defines either an endpoint or an intermediate point
 * in the path we will take to insert an edge constraint.
 * Each such point will either be
 * (a) a vertex or
 * (b) a fraction lambda (0 < lambda < 1) along some #SymEdge.]
 *
 * 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' wlll be NULL for the last one.
 *
 * For case (b), vert will be NULL at first, and later filled in with the created split vertex,
 * and 'in' will be the #SymEdge that we go through, and lambda will be between 0 and 1,
 * the fraction from in's vert to in->next's vert to put the split vertex.
 * 'out' is not needed in this case, since the attachment point will be the sym of the first
 * half of the split edge.
 */
typedef struct CrossData {
  double lambda;
  CDTVert *vert;
  SymEdge *in;
  SymEdge *out;
} CrossData;

static bool get_next_crossing_from_vert(CDT_state *cdt,
                                        CrossData *cd,
                                        CrossData *cd_next,
                                        const CDTVert *v2);

/**
 * As part of finding crossings, we found a case where the next crossing goes through vert v.
 * If it came from a previous vert in cd, then cd_out is the edge that leads from that to v.
 * Else cd_out can be NULL, because it won't be used.
 * Set *cd_next to indicate this. We can set 'in' but not 'out'.  We can set the 'out' of the
 * current cd.
 */
static void fill_crossdata_for_through_vert(CDTVert *v,
                                            SymEdge *cd_out,
                                            CrossData *cd,
                                            CrossData *cd_next)
{
  SymEdge *se;
#ifdef DEBUG_CDT
  int dbg_level = 0;
#endif

  cd_next->lambda = 0.0;
  cd_next->vert = v;
  cd_next->in = NULL;
  cd_next->out = NULL;
  if (cd->lambda == 0.0) {
    cd->out = cd_out;
  }
  else {
    /* One of the edges in the triangle with edge sym(cd->in) contains v. */
    se = sym(cd->in);
    if (se->vert != v) {
      se = se->next;
      if (se->vert != v) {
        se = se->next;
      }
    }
    BLI_assert(se->vert == v);
    cd_next->in = se;
  }
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    dump_cross_data(cd, "cd through vert, cd");
    dump_cross_data(cd_next, "cd_next through vert, cd");
  }
#endif
}

/**
 * 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.
 */
static void fill_crossdata_for_intersect(CDT_state *cdt,
                                         const double *curco,
                                         const CDTVert *v2,
                                         SymEdge *t,
                                         CrossData *cd,
                                         CrossData *cd_next)
{
  CDTVert *va, *vb, *vc;
  double lambda, mu, len_ab;
  SymEdge *se_vcva, *se_vcvb;
  int isect;
#ifdef DEBUG_CDT
  int dbg_level = 0;
#endif

  va = t->vert;
  vb = t->next->vert;
  vc = t->next->next->vert;
  se_vcvb = sym(t->next);
  se_vcva = t->next->next;
  BLI_assert(se_vcva->vert == vc && se_vcva->next->vert == va);
  BLI_assert(se_vcvb->vert == vc && se_vcvb->next->vert == vb);
  UNUSED_VARS_NDEBUG(vc);
  isect = isect_seg_seg_v2_lambda_mu_db(va->co, vb->co, curco, v2->co, &lambda, &mu);
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    double co[2];
    fprintf(stderr, "crossdata for intersect gets lambda=%.17g, mu=%.17g\n", lambda, mu);
    fprintf(stderr,
            "isect=%s\n",
            isect == 2 ? "cross" : (isect == 1 ? "exact" : (isect == 0 ? "none" : "colinear")));
    fprintf(stderr,
            "va=v%d=(%g,%g), vb=v%d=(%g,%g), vc=v%d, curco=(%g,%g), v2=(%g,%g)\n",
            va->index,
            F2(va->co),
            vb->index,
            F2(vb->co),
            vc->index,
            F2(curco),
            F2(v2->co));
    dump_se_short(se_vcva, "vcva=");
    dump_se_short(se_vcvb, " vcvb=");
    interp_v2_v2v2_db(co, va->co, vb->co, lambda);
    fprintf(stderr, "\nco=(%.17g,%.17g)\n", F2(co));
  }
#endif
  switch (isect) {
    case ISECT_LINE_LINE_CROSS:
      len_ab = len_v2v2_db(va->co, vb->co);
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr,
                "len_ab=%g, near a=%g, near b=%g\n",
                len_ab,
                lambda * len_ab,
                (1.0 - lambda) * len_ab);
      }
#endif
      if (lambda * len_ab <= cdt->epsilon) {
        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
      }
      else if ((1.0 - lambda) * len_ab <= cdt->epsilon) {
        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
      }
      else {
        *cd_next = (CrossData){lambda, NULL, t, NULL};
        if (cd->lambda == 0.0) {
          cd->out = se_vcva;
        }
      }
      break;
    case ISECT_LINE_LINE_EXACT:
      if (lambda == 0.0) {
        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
      }
      else if (lambda == 1.0) {
        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
      }
      else {
        *cd_next = (CrossData){lambda, NULL, t, NULL};
        if (cd->lambda == 0.0) {
          cd->out = se_vcva;
        }
      }
      break;
    case ISECT_LINE_LINE_NONE:
      /* It should be very near one end or other of segment. */
      if (lambda <= 0.5) {
        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
      }
      else {
        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
      }
      break;
    case ISECT_LINE_LINE_COLINEAR:
      if (len_squared_v2v2_db(va->co, v2->co) <= len_squared_v2v2_db(vb->co, v2->co)) {
        fill_crossdata_for_through_vert(va, se_vcva, cd, cd_next);
      }
      else {
        fill_crossdata_for_through_vert(vb, se_vcvb, cd, cd_next);
      }
      break;
  }
}

/**
 * As part of finding the crossings of a ray to v2, find the next crossing after 'cd', assuming
 * 'cd' represents a crossing that goes through a vertex.
 *
 * We do a rotational scan around cd's vertex, looking for the triangle where the ray from cd->vert
 * to v2 goes between the two arms from cd->vert, or where it goes along one of the edges.
 */
static bool get_next_crossing_from_vert(CDT_state *cdt,
                                        CrossData *cd,
                                        CrossData *cd_next,
                                        const CDTVert *v2)
{
  SymEdge *t, *tstart;
  CDTVert *vcur, *va, *vb;
  double orient1, orient2;
  bool ok;
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nget_next_crossing_from_vert\n");
    dump_v(cd->vert, "  cd->vert");
  }
#endif

  t = tstart = cd->vert->symedge;
  vcur = cd->vert;
  ok = false;
  do {
    /*
     * The ray from vcur to v2 has to go either between two successive
     * edges around vcur or exactly along them. This time through the
     * loop, check to see if the ray goes along vcur-va
     * or between vcur-va and vcur-vb, where va is the end of t
     * and vb is the next vertex (on the next rot edge around vcur, but
     * should also be the next vert of triangle starting with vcur-va.
     */
    if (t->face != cdt->outer_face && tri_orient(t) < 0.0) {
      fprintf(stderr, "BAD TRI orientation %g\n", tri_orient(t)); /* TODO: handle this */
#ifdef DEBUG_CDT
      dump_se_cycle(t, "top of ccw scan loop", 4);
#endif
    }
    va = t->next->vert;
    vb = t->next->next->vert;
    orient1 = orient2d(t->vert->co, va->co, v2->co);
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "non-final through vert case\n");
      dump_v(va, " va");
      dump_v(vb, " vb");
      fprintf(stderr, "orient1=%g\n", orient1);
    }
#endif
    if (orient1 == 0.0 && in_line(vcur->co, va->co, v2->co)) {
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "ray goes through va\n");
      }
#endif
      fill_crossdata_for_through_vert(va, t, cd, cd_next);
      ok = true;
      break;
    }
    if (t->face != cdt->outer_face) {
      orient2 = orient2d(vcur->co, vb->co, v2->co);
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "orient2=%g\n", orient2);
      }
#endif
      /* Don't handle orient2 == 0.0 case here: next rotation will get it. */
      if (orient1 > 0.0 && orient2 < 0.0) {
#ifdef DEBUG_CDT
        if (dbg_level > 1) {
          fprintf(stderr, "segment intersects\n");
        }
#endif
        t = t->next;
        fill_crossdata_for_intersect(cdt, vcur->co, v2, t, cd, cd_next);
        ok = true;
        break;
      }
    }
  } while ((t = t->rot) != tstart);
#ifdef DEBUG_CDT
  if (!ok) {
    /* Didn't find the exit! Shouldn't happen. */
    fprintf(stderr, "shouldn't happen: can't find next crossing from vert\n");
  }
#endif
  return ok;
}

/**
 * As part of finding the crossings of a ray to 'v2', find the next crossing after 'cd', assuming
 * 'cd' represents a crossing that goes through a an edge, not at either end of that edge.
 *
 * We have the triangle 'vb-va-vc', where va and vb are the split edge and 'vc' is the third vertex
 * on that new side of the edge (should be closer to v2). The next crossing should be through 'vc'
 * or intersecting 'vb-vc' or 'va-vc'.
 */
static void get_next_crossing_from_edge(CDT_state *cdt,
                                        CrossData *cd,
                                        CrossData *cd_next,
                                        const CDTVert *v2)
{
  double curco[2];
  double orient;
  CDTVert *va, *vb, *vc;
  SymEdge *se_ac;
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nget_next_crossing_from_edge\n");
    fprintf(stderr, "  lambda=%.17g\n", cd->lambda);
    dump_se_short(cd->in, "  cd->in");
  }
#endif

  va = cd->in->vert;
  vb = cd->in->next->vert;
  interp_v2_v2v2_db(curco, va->co, vb->co, cd->lambda);
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "  curco=(%.17g,%.17g)\n", F2(curco));
  }
#endif
  se_ac = sym(cd->in)->next;
  vc = se_ac->next->vert;
  orient = orient2d(curco, v2->co, vc->co);
#ifdef DEBUG_CDT
  if (dbg_level > 1) {
    fprintf(stderr, "now searching with third vertex ");
    dump_v(vc, "vc");
    fprintf(stderr, "orient2d(cur, v2, vc) = %g\n", orient);
  }
#endif
  if (orient < 0.0) {
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "curco--v2 intersects vb--vc\n");
    }
#endif
    fill_crossdata_for_intersect(cdt, curco, v2, se_ac->next, cd, cd_next);
  }
  else if (orient > 0.0) {
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "curco--v2 intersects va--vc\n");
    }
#endif
    fill_crossdata_for_intersect(cdt, curco, v2, se_ac, cd, cd_next);
  }
  else {
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "orient==0 case, so going through or to vc\n");
    }
#endif
    *cd_next = (CrossData){0.0, vc, se_ac->next, NULL};
  }
}

/**
 * Add a constrained edge between v1 and v2 to cdt structure.
 * This may result in a number of #CDTEdges created, due to intersections
 * and partial overlaps with existing cdt vertices and edges.
 * Each created #CDTEdge will have input_id added to its input_ids list.
 *
 * If \a r_edges is not NULL, the #CDTEdges generated or found that go from
 * v1 to v2 are put into that linked list, in order.
 *
 * Assumes that #BLI_constrained_delaunay_get_output has not been called yet.
 */
static void add_edge_constraint(
    CDT_state *cdt, CDTVert *v1, CDTVert *v2, int input_id, LinkNode **r_edges)
{
  SymEdge *t, *se, *tstart, *tnext;
  int i, j, n, visit;
  bool ok;
  CrossData *crossings = NULL;
  CrossData *cd, *cd_prev, *cd_next;
  CDTVert *v;
  CDTEdge *edge;
  LinkNodePair edge_list = {NULL, NULL};
  BLI_array_staticdeclare(crossings, 128);
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nADD_EDGE_CONSTRAINT\n");
    dump_v(v1, "  1");
    dump_v(v2, "  2");
  }
#endif

  if (r_edges) {
    *r_edges = NULL;
  }

  /*
   * Handle two special cases first:
   * 1) The two end vertices are the same (can happen because of merging).
   * 2) There is already an edge between v1 and v2.
   */
  if (v1 == v2) {
    return;
  }
  if ((t = find_symedge_between_verts(v1, v2)) != NULL) {
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "segment already there\n");
    }
#endif
    add_to_input_ids(&t->edge->input_ids, input_id, cdt);
    if (r_edges != NULL) {
      BLI_linklist_append_pool(&edge_list, t->edge, cdt->listpool);
      *r_edges = edge_list.list;
    }
    return;
  }

  /*
   * Fill crossings array with CrossData points for intersection path from v1 to v2.
   *
   * At every point, the crossings array has the path so far, except that
   * the .out field of the last element of it may not be known yet -- if that
   * last element is a vertex, then we won't know the output edge until we
   * find the next crossing.
   *
   * To protect against infinite loops, we keep track of which vertices
   * we have visited by setting their visit_index to a new visit epoch.
   *
   * We check a special case first: where the segment is already there in
   * one hop. Saves a bunch of orient2d tests in that common case.
   */
  visit = ++cdt->visit_count;
  BLI_array_grow_one(crossings);
  crossings[0] = (CrossData){0.0, v1, NULL, NULL};
  while (!((n = BLI_array_len(crossings)) > 0 && crossings[n - 1].vert == v2)) {
    BLI_array_grow_one(crossings);
    cd = &crossings[n - 1];
    cd_next = &crossings[n];
    if (crossings[n - 1].lambda == 0.0) {
      ok = get_next_crossing_from_vert(cdt, cd, cd_next, v2);
    }
    else {
      get_next_crossing_from_edge(cdt, cd, cd_next, v2);
    }
    if (!ok || BLI_array_len(crossings) == 100000) {
      /* Shouldn't happen but if does, just bail out. */
#ifdef DEBUG_CDT
      fprintf(stderr, "FAILURE adding segment, bailing out\n");
#endif
      BLI_array_free(crossings);
      return;
    }
#ifdef DEBUG_CDT
    if (dbg_level > 1) {
      fprintf(stderr, "crossings[%d]: ", n);
      dump_cross_data(&crossings[n], "");
    }
#endif
    if (crossings[n].lambda == 0.0) {
      if (crossings[n].vert->visit_index == visit) {
        fprintf(stderr, "WHOOPS, REVISIT. Bailing out.\n"); /*TODO: desperation search here. */
        BLI_array_free(crossings);
        return;
      }
      crossings[n].vert->visit_index = visit;
    }
  }

#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "\ncrossings found\n");
    for (i = 0; i < BLI_array_len(crossings); i++) {
      fprintf(stderr, "%d: ", i);
      dump_cross_data(&crossings[i], "");
      if (crossings[i].lambda == 0.0) {
        if (i == 0 || crossings[i - 1].lambda == 0.0) {
          BLI_assert(crossings[i].in == NULL);
        }
        else {
          BLI_assert(crossings[i].in != NULL && crossings[i].in->vert == crossings[i].vert);
          BLI_assert(crossings[i].in->face == sym(crossings[i - 1].in)->face);
        }
        if (i == BLI_array_len(crossings) - 1) {
          BLI_assert(crossings[i].vert == v2);
          BLI_assert(crossings[i].out == NULL);
        }
        else {
          BLI_assert(crossings[i].out->vert == crossings[i].vert);
          if (crossings[i + 1].lambda == 0.0) {
            BLI_assert(crossings[i].out->next->vert == crossings[i + 1].vert);
          }
          else {
            BLI_assert(crossings[i].out->face == crossings[i + 1].in->face);
          }
        }
      }
      else {
        if (i > 0 && crossings[i - 1].lambda == 0.0) {
          BLI_assert(crossings[i].in->face == crossings[i - 1].out->face);
        }
        BLI_assert(crossings[i].out == NULL);
      }
    }
  }
#endif

  /*
   * Postprocess crossings.
   * Some crossings may have an intersection crossing followed
   * by a vertex crossing that is on the same edge that was just
   * intersected. We prefer to go directly from the previous
   * crossing directly to the vertex. This may chain backwards.
   *
   * This loop marks certain crossings as "deleted", by setting
   * their lambdas to -1.0.
   */
  for (i = 2; i < BLI_array_len(crossings); i++) {
    cd = &crossings[i];
    if (cd->lambda == 0.0) {
      v = cd->vert;
      for (j = i - 1; j > 0; j--) {
        cd_prev = &crossings[j];
        if ((cd_prev->lambda == 0.0 && cd_prev->vert != v) ||
            (cd_prev->lambda != 0.0 && cd_prev->in->vert != v && cd_prev->in->next->vert != v)) {
          break;
        }
        cd_prev->lambda = -1.0; /* Mark cd_prev as 'deleted'. */
#ifdef DEBUG_CDT
        if (dbg_level > 0) {
          fprintf(stderr, "deleted crossing %d\n", j);
        }
#endif
      }
      if (j < i - 1) {
        /* Some crossings were deleted. Fix the in and out edges across gap. */
        cd_prev = &crossings[j];
        if (cd_prev->lambda == 0.0) {
          se = find_symedge_between_verts(cd_prev->vert, v);
          if (se == NULL) {
#ifdef DEBUG_CDT
            fprintf(stderr, "FAILURE(a) in delete crossings, bailing out.\n");
#endif
            BLI_array_free(crossings);
            return;
          }
          cd_prev->out = se;
          cd->in = NULL;
        }
        else {
          se = find_symedge_with_face(v, sym(cd_prev->in)->face);
          if (se == NULL) {
#ifdef DEBUG_CDT
            fprintf(stderr, "FAILURE(b) in delete crossings, bailing out.\n");
#endif
            BLI_array_free(crossings);
            return;
          }
          cd->in = se;
        }
#ifdef DEBUG_CDT
        if (dbg_level > 0) {
          fprintf(stderr, "after deleting crossings:\n");
          fprintf(stderr, "cross[%d]: ", j);
          dump_cross_data(cd_prev, "");
          fprintf(stderr, "cross[%d]: ", i);
          dump_cross_data(cd, "");
        }
#endif
      }
    }
  }

  /*
   * Insert all intersection points on constrained edges.
   */
  for (i = 0; i < BLI_array_len(crossings); i++) {
    cd = &crossings[i];
    if (cd->lambda != 0.0 && cd->lambda != -1.0 && is_constrained_edge(cd->in->edge)) {
      edge = split_edge(cdt, cd->in, cd->lambda);
      cd->vert = edge->symedges[0].vert;
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "insert vert for crossing %d: ", i);
        dump_v(cd->vert, "inserted");
      }
#endif
    }
  }

  /*
   * Remove any crossed, non-intersected edges.
   */
  for (i = 0; i < BLI_array_len(crossings); i++) {
    cd = &crossings[i];
    if (cd->lambda != 0.0 && cd->lambda != -1.0 && !is_constrained_edge(cd->in->edge)) {
      delete_edge(cdt, cd->in);
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "delete edge for crossing %d\n", i);
      }
#endif
    }
  }

  /*
   * Insert segments for v1->v2.
   */
  tstart = crossings[0].out;
  for (i = 1; i < BLI_array_len(crossings); i++) {
    cd = &crossings[i];
    if (cd->lambda == -1.0) {
      continue; /* This crossing was deleted. */
    }
    t = tnext = NULL;
    if (cd->lambda != 0.0) {
      if (is_constrained_edge(cd->in->edge)) {
        t = cd->vert->symedge;
        tnext = sym(t)->next;
      }
    }
    else if (cd->lambda == 0.0) {
      t = cd->in;
      tnext = cd->out;
      if (t == NULL) {
        /* Previous non-deleted crossing must also have been a vert, and segment should exist. */
        for (j = i - 1; j >= 0; j--) {
          cd_prev = &crossings[j];
          if (cd_prev->lambda != -1.0) {
            break;
          }
        }
        BLI_assert(cd_prev->lambda == 0.0);
        BLI_assert(cd_prev->out->next->vert == cd->vert);
#ifdef DEBUG_CDT
        if (dbg_level > 1) {
          fprintf(stderr, "edge to crossing %d already there\n", i);
        }
#endif
        edge = cd_prev->out->edge;
        add_to_input_ids(&edge->input_ids, input_id, cdt);
      }
    }
    if (t != NULL) {
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "edge to crossing %d: insert diagonal between\n", i);
        dump_se(tstart, "  ");
        dump_se(t, " ");
        dump_se_cycle(tstart, "tstart", 100);
        dump_se_cycle(t, "t", 100);
      }
#endif
      if (tstart->next->vert == t->vert) {
        edge = tstart->edge;
#ifdef DEBUG_CDT
        if (dbg_level > 1) {
          fprintf(stderr, "already there (b)\n");
        }
#endif
      }
      else {
        edge = add_diagonal(cdt, tstart, t);
      }
      add_to_input_ids(&edge->input_ids, input_id, cdt);
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        fprintf(stderr, "added\n");
      }
#endif
      if (r_edges != NULL) {
        BLI_linklist_append_pool(&edge_list, edge, cdt->listpool);
      }
      /* Now retriangulate upper and lower gaps. */
      re_delaunay_triangulate(cdt, &edge->symedges[0]);
      re_delaunay_triangulate(cdt, &edge->symedges[1]);
    }
    if (i < BLI_array_len(crossings) - 1) {
      if (tnext != NULL) {
        tstart = tnext;
#ifdef DEBUG_CDT
        if (dbg_level > 1) {
          fprintf(stderr, "now tstart = ");
          dump_se(tstart, "");
        }
#endif
      }
    }
  }

  if (r_edges) {
    *r_edges = edge_list.list;
  }
  BLI_array_free(crossings);
}

/**
 * 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 untagged. 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.
 */
static void add_face_ids(
    CDT_state *cdt, SymEdge *face_symedge, int face_id, int fedge_start, int fedge_end)
{
  Stack stack;
  SymEdge *se, *se_start, *se_sym;
  CDTFace *face, *face_other;
  int visit;

  /* Can't loop forever since eventually would visit every face. */
  cdt->visit_count++;
  visit = cdt->visit_count;
  stack = NULL;
  push(&stack, face_symedge, cdt);
  while (!is_empty(&stack)) {
    se = pop(&stack, cdt);
    face = se->face;
    if (face->visit_index == visit) {
      continue;
    }
    face->visit_index = visit;
    add_to_input_ids(&face->input_ids, face_id, cdt);
    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)) {
        se_sym = sym(se);
        face_other = se_sym->face;
        if (face_other->visit_index != visit) {
          push(&stack, se_sym, cdt);
        }
      }
    }
  }
}

/* Delete_edge but try not to mess up outer face.
 * Also faces have symedges now, so make sure not
 * to mess those up either. */
static void dissolve_symedge(CDT_state *cdt, SymEdge *se)
{
  SymEdge *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;
    }
  }
  delete_edge(cdt, se);
}

/* Slow way to get face and start vertex index within face for edge id e. */
static bool get_face_edge_id_indices(const CDT_input *in, int e, int *r_f, int *r_fi)
{
  int f;
  int id;

  id = in->edges_len;
  if (e < id) {
    return false;
  }
  for (f = 0; f < in->faces_len; f++) {
    if (e < id + in->faces_len_table[f]) {
      *r_f = f;
      *r_fi = e - id;
      return true;
    }
    id += in->faces_len_table[f];
  }
  return false;
}

/* Is pt_co when snapped to segment seg1 seg2 all of:
 * a) strictly within that segment
 * b) within epsilon from original pt_co
 * c) pt_co is not within epsilon of either seg1 or seg2.
 * Return true if so, and return in *r_lambda the fraction of the way from seg1 to seg2 of the
 * snapped point.
 */
static bool check_vert_near_segment(const double *pt_co,
                                    const double *seg1,
                                    const double *seg2,
                                    double epsilon_squared,
                                    double *r_lambda)
{
  double lambda, snap_co[2];

  lambda = closest_to_line_v2_db(snap_co, pt_co, seg1, seg2);
  *r_lambda = lambda;
  if (lambda <= 0.0 || lambda >= 1.0) {
    return false;
  }
  if (len_squared_v2v2_db(pt_co, snap_co) > epsilon_squared) {
    return false;
  }
  if (len_squared_v2v2_db(pt_co, seg1) <= epsilon_squared ||
      len_squared_v2v2_db(pt_co, seg2) <= epsilon_squared) {
    return false;
  }
  return true;
}

typedef struct EdgeVertLambda {
  int e_id;
  int v_id;
  double lambda;
} EdgeVertLambda;

/* For sorting first by edge id, then by lambda, then by vert id. */
static int evl_cmp(const void *a, const void *b)
{
  const EdgeVertLambda *area = a;
  const EdgeVertLambda *sb = b;

  if (area->e_id < sb->e_id) {
    return -1;
  }
  if (area->e_id > sb->e_id) {
    return 1;
  }
  if (area->lambda < sb->lambda) {
    return -1;
  }
  if (area->lambda > sb->lambda) {
    return 1;
  }
  if (area->v_id < sb->v_id) {
    return -1;
  }
  if (area->v_id > sb->v_id) {
    return 1;
  }
  return 0;
}

/**
 * If epsilon > 0, and input doesn't have skip_modify_input == true,
 * check input to see if any constraint edge ends (including face edges) come
 * within epsilon of another edge.
 * For all such cases, we want to split the constraint edge at the point nearest to near vertex
 * and move the vertex coordinates to be on that edge.
 * But exclude cases where they come within epsilon of either end because those will be handled
 * by vertex merging in the main triangulation algorithm.
 *
 * If any such splits are found, make a new CDT_input reflecting this change, and provide an
 * edge map to map from edge ids in the new input space to edge ids in the old input space.
 *
 * TODO: replace naive O(n^2) algorithm with kdopbvh-based one.
 */
static const CDT_input *modify_input_for_near_edge_ends(const CDT_input *input, int **r_edge_map)
{
  CDT_input *new_input = NULL;
  int e, eprev, e1, e2, f, fi, flen, start, i, j;
  int i_new, i_old, i_evl;
  int v11, v12, v21, v22;
  double co11[2], co12[2], co21[2], co22[2];
  double lambda;
  double eps = (double)input->epsilon;
  double eps_sq = eps * eps;
  int tot_edge_constraints, edges_len, tot_face_edges;
  int new_tot_face_edges, new_tot_con_edges;
  int delta_con_edges, delta_face_edges, cur_e_cnt;
  int *edge_map;
  int evl_len;
  EdgeVertLambda *edge_vert_lambda = NULL;
  BLI_array_staticdeclare(edge_vert_lambda, 128);
#ifdef DEBUG_CDT
  EdgeVertLambda *evl;
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nMODIFY INPUT\n\n");
  }
#endif

  *r_edge_map = NULL;
  if (input->epsilon == 0.0 || input->skip_input_modify ||
      (input->edges_len == 0 && input->faces_len == 0)) {
    return input;
  }

  /* Edge constraints are union of the explicitly provided edges and the implicit edges
   * that are part of the provided faces. We index constraints by have the first input->edges_len
   * ints standing for the explicit edge with the same index, and the rest being face edges in
   * the order that the faces appear and then edges within those faces, with indices offset by
   * input->edges_len.
   * Calculate tot_edge_constraints to be the sum of the two kinds of edges.
   * We first have to count the number of face edges.
   * That is the same as the number of vertices in the faces table, which
   * we can find by adding the last length to the last start.
   */
  edges_len = input->edges_len;
  tot_edge_constraints = edges_len;
  if (input->faces_len > 0) {
    tot_face_edges = input->faces_start_table[input->faces_len - 1] +
                     input->faces_len_table[input->faces_len - 1];
  }
  else {
    tot_face_edges = 0;
  }
  tot_edge_constraints = edges_len + tot_face_edges;

  for (e1 = 0; e1 < tot_edge_constraints - 1; e1++) {
    if (e1 < edges_len) {
      v11 = input->edges[e1][0];
      v12 = input->edges[e1][1];
    }
    else {
      if (!get_face_edge_id_indices(input, e1, &f, &fi)) {
        /* Must be bad input. Will be caught later so don't need to signal here. */
        continue;
      }
      start = input->faces_start_table[f];
      flen = input->faces_len_table[f];
      v11 = input->faces[start + fi];
      v12 = input->faces[(fi == flen - 1) ? start : start + fi + 1];
    }
    for (e2 = e1 + 1; e2 < tot_edge_constraints; e2++) {
      if (e2 < edges_len) {
        v21 = input->edges[e2][0];
        v22 = input->edges[e2][1];
      }
      else {
        if (!get_face_edge_id_indices(input, e2, &f, &fi)) {
          continue;
        }
        start = input->faces_start_table[f];
        flen = input->faces_len_table[f];
        v21 = input->faces[start + fi];
        v22 = input->faces[(fi == flen - 1) ? start : start + fi + 1];
      }
      copy_v2db_v2fl(co11, input->vert_coords[v11]);
      copy_v2db_v2fl(co12, input->vert_coords[v12]);
      copy_v2db_v2fl(co21, input->vert_coords[v21]);
      copy_v2db_v2fl(co22, input->vert_coords[v22]);
      if (check_vert_near_segment(co11, co21, co22, eps_sq, &lambda)) {

        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e2, v11, lambda}));
      }
      if (check_vert_near_segment(co12, co21, co22, eps_sq, &lambda)) {
        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e2, v12, lambda}));
      }
      if (check_vert_near_segment(co21, co11, co12, eps_sq, &lambda)) {
        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e1, v21, lambda}));
      }
      if (check_vert_near_segment(co22, co11, co12, eps_sq, &lambda)) {
        BLI_array_append(edge_vert_lambda, ((EdgeVertLambda){e1, v22, lambda}));
      }
    }
  }

  evl_len = BLI_array_len(edge_vert_lambda);
  if (evl_len > 0) {
    /* Sort to bring splits for each edge together,
     * and for each edge, to be in order of lambda. */
    qsort(edge_vert_lambda, evl_len, sizeof(EdgeVertLambda), evl_cmp);
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "\nafter sorting\n");
      for (i = 0; i < evl_len; i++) {
        evl = &edge_vert_lambda[i];
        fprintf(stderr, "e%d, v%d, %g\n", evl->e_id, evl->v_id, evl->lambda);
      }
    }
#endif

    /* Remove dups in edge_vert_lambda, where dup means that the edge is the
     * same, and the verts are either the same or will be merged by epsilon-nearness.
     */
    i = 0;
    j = 0;
    /* In loop, copy from position j to position i. */
    for (j = 0; j < evl_len;) {
      int k;
      if (i != j) {
        memmove(&edge_vert_lambda[i], &edge_vert_lambda[j], sizeof(EdgeVertLambda));
      }
      for (k = j + 1; k < evl_len; k++) {
        int vj = edge_vert_lambda[j].v_id;
        int vk = edge_vert_lambda[k].v_id;
        if (vj != vk) {
          if (len_squared_v2v2(input->vert_coords[vj], input->vert_coords[vk]) > (float)eps_sq) {
            break;
          }
        }
      }
      j = k;
      i++;
    }

    if (i != evl_len) {
      evl_len = i;
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr, "\nduplicates eliminated\n");
        for (i = 0; i < evl_len; i++) {
          evl = &edge_vert_lambda[i];
          fprintf(stderr, "e%d, v%d, %g\n", evl->e_id, evl->v_id, evl->lambda);
        }
      }
#endif
    }
    /* Find delta in number of constraint edges and face edges.
     * This may be overestimates of true number, due to duplicates. */
    delta_con_edges = 0;
    delta_face_edges = 0;
    cur_e_cnt = 0;
    eprev = -1;
    for (i = 0; i < evl_len; i++) {
      e = edge_vert_lambda[i].e_id;
      if (i > 0 && e > eprev) {
        /* New edge group. Previous group had cur_e_cnt split vertices.
         * That is the delta in the number of edges needed in input since
         * there will be cur_e_cnt + 1 edges replacing one edge.
         */
        if (eprev < edges_len) {
          delta_con_edges += cur_e_cnt;
        }
        else {
          delta_face_edges += cur_e_cnt;
        }
        cur_e_cnt = 1;
        ;
      }
      else {
        cur_e_cnt++;
      }
      eprev = e;
    }
    if (eprev < edges_len) {
      delta_con_edges += cur_e_cnt;
    }
    else {
      delta_face_edges += cur_e_cnt;
    }
    new_tot_con_edges = input->edges_len + delta_con_edges;
    if (input->faces_len > 0) {
      new_tot_face_edges = input->faces_start_table[input->faces_len - 1] +
                           input->faces_len_table[input->faces_len - 1] + delta_face_edges;
    }
    else {
      new_tot_face_edges = 0;
    }

    /* Allocate new CDT_input, now we know sizes needed (perhaps overestimated a bit).
     * Caller will be responsible for freeing it and its arrays.
     */
    new_input = MEM_callocN(sizeof(CDT_input), __func__);
    new_input->epsilon = input->epsilon;
    new_input->verts_len = input->verts_len;
    new_input->vert_coords = (float(*)[2])MEM_malloc_arrayN(
        new_input->verts_len, 2 * sizeof(float), __func__);
    /* We don't do it now, but may decide to change coords of snapped verts. */
    memmove(new_input->vert_coords,
            input->vert_coords,
            (size_t)new_input->verts_len * sizeof(float) * 2);

    if (edges_len > 0) {
      new_input->edges_len = new_tot_con_edges;
      new_input->edges = (int(*)[2])MEM_malloc_arrayN(
          new_tot_con_edges, 2 * sizeof(int), __func__);
    }

    if (input->faces_len > 0) {
      new_input->faces_len = input->faces_len;
      new_input->faces_start_table = (int *)MEM_malloc_arrayN(
          new_input->faces_len, sizeof(int), __func__);
      new_input->faces_len_table = (int *)MEM_malloc_arrayN(
          new_input->faces_len, sizeof(int), __func__);
      new_input->faces = (int *)MEM_malloc_arrayN(new_tot_face_edges, sizeof(int), __func__);
    }

    edge_map = (int *)MEM_malloc_arrayN(
        new_tot_con_edges + new_tot_face_edges, sizeof(int), __func__);
    *r_edge_map = edge_map;

    i_new = i_old = i_evl = 0;
    e = edge_vert_lambda[0].e_id;
    /* First do new constraint edges. */
    for (i_old = 0; i_old < edges_len; i_old++) {
      if (i_old < e) {
        /* Edge for i_old not split; copy it into new_input. */
        new_input->edges[i_new][0] = input->edges[i_old][0];
        new_input->edges[i_new][1] = input->edges[i_old][1];
        edge_map[i_new] = i_old;
        i_new++;
      }
      else {
        /* Edge for i_old is split. */
        BLI_assert(i_old == e);
        new_input->edges[i_new][0] = input->edges[i_old][0];
        new_input->edges[i_new][1] = edge_vert_lambda[i_evl].v_id;
        edge_map[i_new] = i_old;
        i_new++;
        i_evl++;
        while (i_evl < evl_len && e == edge_vert_lambda[i_evl].e_id) {
          new_input->edges[i_new][0] = new_input->edges[i_new - 1][1];
          new_input->edges[i_new][1] = edge_vert_lambda[i_evl].v_id;
          edge_map[i_new] = i_old;
          i_new++;
          i_evl++;
        }
        new_input->edges[i_new][0] = new_input->edges[i_new - 1][1];
        new_input->edges[i_new][1] = input->edges[i_old][1];
        edge_map[i_new] = i_old;
        i_new++;
        if (i_evl < evl_len) {
          e = edge_vert_lambda[i_evl].e_id;
        }
        else {
          e = INT_MAX;
        }
      }
    }
    BLI_assert(i_new <= new_tot_con_edges);
    new_input->edges_len = i_new;

    /* Now do face constraints. */
    if (input->faces_len > 0) {
      f = 0;
      i_new = 0; /* Now will index cur place in new_input->faces. */
      while (i_old < tot_edge_constraints) {
        flen = input->faces_len_table[f];
        BLI_assert(i_old - edges_len == input->faces_start_table[f]);
        new_input->faces_start_table[f] = i_new;
        if (i_old + flen - 1 < e) {
          /* Face f is not split. */
          for (j = 0; j < flen; j++) {
            new_input->faces[i_new] = input->faces[i_old - edges_len + j];
            edge_map[i_new + new_input->edges_len] = i_old + j;
            i_new++;
          }
          i_old += flen;
          new_input->faces_len_table[f] = flen;
          f++;
        }
        else {
          /* Face f has at least one split edge. */
          int i_new_start = i_new;
          for (j = 0; j < flen; j++) {
            if (i_old + j < e) {
              /* jth edge of f is not split. */
              new_input->faces[i_new] = input->faces[i_old - edges_len + j];
              edge_map[i_new + new_input->edges_len] = i_old + j;
              i_new++;
            }
            else {
              /* jth edge of f is split. */
              BLI_assert(i_old + j == e);
              new_input->faces[i_new] = input->faces[i_old - edges_len + j];
              edge_map[i_new + new_input->edges_len] = i_old + j;
              i_new++;
              while (i_evl < evl_len && e == edge_vert_lambda[i_evl].e_id) {
                new_input->faces[i_new] = edge_vert_lambda[i_evl].v_id;
                edge_map[i_new + new_input->edges_len] = i_old + j;
                i_new++;
                i_evl++;
              }
              if (i_evl < evl_len) {
                e = edge_vert_lambda[i_evl].e_id;
              }
              else {
                e = INT_MAX;
              }
            }
          }
          new_input->faces_len_table[f] = i_new - i_new_start;
          i_old += flen;
          f++;
        }
      }
    }

#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "\nnew constraint edges\n");
      for (i = 0; i < new_input->edges_len; i++) {
        fprintf(stderr, "  e%d: (%d,%d)\n", i, new_input->edges[i][0], new_input->edges[i][1]);
      }
      fprintf(stderr, "\nnew faces\n");
      for (f = 0; f < new_input->faces_len; f++) {
        flen = new_input->faces_len_table[f];
        start = new_input->faces_start_table[f];
        fprintf(stderr, "  f%d: start=%d, len=%d\n    ", f, start, flen);
        for (i = start; i < start + flen; i++) {
          fprintf(stderr, "%d ", new_input->faces[i]);
        }
        fprintf(stderr, "\n");
      }
      fprintf(stderr, "\nedge map (new->old)\n");
      for (i = 0; i < new_tot_con_edges + new_tot_face_edges; i++) {
        fprintf(stderr, "  %d->%d\n", i, edge_map[i]);
      }
    }
#endif
  }

  BLI_array_free(edge_vert_lambda);
  if (new_input != NULL) {
    return (const CDT_input *)new_input;
  }
  return input;
}

static void free_modified_input(CDT_input *input)
{
  MEM_freeN(input->vert_coords);
  if (input->edges != NULL) {
    MEM_freeN(input->edges);
  }
  if (input->faces != NULL) {
    MEM_freeN(input->faces);
    MEM_freeN(input->faces_len_table);
    MEM_freeN(input->faces_start_table);
  }
  MEM_freeN(input);
}

/* Return true if we can merge se's vert into se->next's vert
 * without making the area of any new triangle formed by doing
 * that into a zero or negative area triangle.*/
static bool can_collapse(const SymEdge *se)
{
  SymEdge *loop_se;
  const double *co = se->next->vert->co;

  for (loop_se = se->rot; loop_se != se && loop_se->rot != se; loop_se = loop_se->rot) {
    if (orient2d(co, loop_se->next->vert->co, loop_se->rot->next->vert->co) <= 0.0) {
      return false;
    }
  }
  return true;
}

/*
 * Merge one end of e onto the other, fixing up surrounding faces.
 *
 * General situation looks something like:
 *
 *             c-----e
 *           /  \   / \
 *          /    \ /   \
 *         a------b-----f
 *          \    / \   /
 *           \  /   \ /
 *             d-----g
 *
 * where ab is the tiny edge. We want to merge a and b and delete edge ab.
 * We don't want to change the coordinates of input vertices [We could revisit this
 * in the future, as API def doesn't prohibit this, but callers will appreciate if they
 * don't change.]
 * Sometimes the collapse shouldn't happen because the triangles formed by the changed
 * edges may end up with zero or negative area (see can_collapse, above).
 * So don't choose a collapse direction that is not allowed or one that has an original vertex
 * as origin and a non-original vertex as destination.
 * If both collapse directions are allowed by that rule, pick the one with the lower original
 * index.
 *
 * After merging, the faces abc and adb disappear (if they are not the outer face).
 * Suppose we merge b onto a.
 * Then edges cb and db are deleted. Face cbe becomes cae and face bdg becomes adg.
 * Any other faces attached to b now have a in their place.
 * We can do this by rotating edges round b, replacing their vert references with a.
 * Similar statements can be made about what happens when a merges into b;
 * in code below we'll swap a and b to make above lettering work for a b->a merge.
 * Return the vert at the collapsed edge, if a collapse happens.
 */
static CDTVert *collapse_tiny_edge(CDT_state *cdt, CDTEdge *e)
{
  CDTVert *va, *vb;
  SymEdge *ab_se, *ba_se, *bd_se, *bc_se, *ad_se, *ac_se;
  SymEdge *bg_se, *be_se, *se, *gb_se, *ca_se;
  bool can_collapse_a_to_b, can_collapse_b_to_a;
#ifdef DEBUG_CDT
  int dbg_level = 0;
#endif

  ab_se = &e->symedges[0];
  ba_se = &e->symedges[1];
  va = ab_se->vert;
  vb = ba_se->vert;
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "\ncollapse_tiny_edge\n");
    dump_se(&e->symedges[0], "tiny edge");
    fprintf(stderr, "a = [%d], b = [%d]\n", va->index, vb->index);
    validate_cdt(cdt, true, false, true);
  }
#endif
  can_collapse_a_to_b = can_collapse(ab_se);
  can_collapse_b_to_a = can_collapse(ba_se);
  /* Now swap a and b if necessary and possible, so that from this point on we are collapsing b to
   * a. */
  if (va->index > vb->index || !can_collapse_b_to_a) {
    if (can_collapse_a_to_b && !(is_original_vert(va, cdt) && !is_original_vert(vb, cdt))) {
      SWAP(CDTVert *, va, vb);
      ab_se = &e->symedges[1];
      ba_se = &e->symedges[0];
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr, "swapped a and b\n");
      }
#endif
    }
    else {
      /* Neither collapse direction is OK. */
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr, "neither collapse direction ok\n");
      }
#endif
      return NULL;
    }
  }
  bc_se = ab_se->next;
  bd_se = ba_se->rot;
  if (bd_se == ba_se) {
    /* Shouldn't happen. Wire edge in outer face. */
    fprintf(stderr, "unexpected wire edge\n");
    return NULL;
  }
  vb->merge_to_index = va->merge_to_index == -1 ? va->index : va->merge_to_index;
  vb->symedge = NULL;
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr,
            "vb = v[%d] merges to va = v[%d], vb->merge_to_index=%d\n",
            vb->index,
            va->index,
            vb->merge_to_index);
  }
#endif
  /* First fix the vertex of intermediate triangles, like bgf. */
  for (se = bd_se->rot; se != bc_se; se = se->rot) {
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      dump_se(se, "intermediate tri edge, setting vert to va");
    }
#endif
    se->vert = va;
  }
  ad_se = sym(sym(bd_se)->rot);
  ca_se = bc_se->next;
  ac_se = sym(ca_se);
  if (bd_se->rot != bc_se) {
    bg_se = bd_se->rot;
    be_se = sym(bc_se)->next;
    gb_se = sym(bg_se);
  }
  else {
    bg_se = NULL;
    be_se = NULL;
  }
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr, "delete bd, inputs to ad\n");
    dump_se(bd_se, "  bd");
    dump_se(ad_se, "  ad");
    fprintf(stderr, "delete bc, inputs to ac\n");
    dump_se(bc_se, "  bc");
    dump_se(ac_se, "  ac");
    fprintf(stderr, "delete ab\n");
    dump_se(ab_se, "  ab");
    if (bg_se != NULL) {
      fprintf(stderr, "fix up bg, be\n");
      dump_se(bg_se, "  bg");
      dump_se(be_se, "  be");
    }
  }
#endif
  add_list_to_input_ids(&ad_se->edge->input_ids, bd_se->edge->input_ids, cdt);
  delete_edge(cdt, bd_se);
  add_list_to_input_ids(&ac_se->edge->input_ids, bc_se->edge->input_ids, cdt);
  delete_edge(cdt, sym(bc_se));
  /* At this point we have this:
   *
   *             c-----e
   *           /      / \
   *          /      /   \
   *         a------b-----f
   *          \      \   /
   *           \      \ /
   *             d-----g
   *
   * Or, if there is not bg_se and be_se, like this:
   *
   *             c
   *           /
   *          /
   *         a------b
   *          \
   *           \
   *             d
   *
   * (But we've already changed the vert field for bg, bf, ..., be to be va.)
   */
  if (bg_se != NULL) {
    gb_se->next = ad_se;
    ad_se->rot = bg_se;
    ca_se->next = be_se;
    be_se->rot = ac_se;
    bg_se->vert = va;
    be_se->vert = va;
  }
  else {
    ca_se->next = ad_se;
    ad_se->rot = ac_se;
  }
  /* Don't use delete_edge as it changes too much. */
  ab_se->next = ab_se->rot = NULL;
  ba_se->next = ba_se->rot = NULL;
  if (va->symedge == ab_se) {
    va->symedge = ac_se;
  }
  return va;
}

/*
 * Check to see if e is tiny (length <= epsilon) and queue it if so.
 */
static void maybe_enqueue_small_feature(CDT_state *cdt, CDTEdge *e, LinkNodePair *tiny_edge_queue)
{
  SymEdge *se, *sesym;
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nmaybe_enqueue_small_features\n");
    dump_se(&e->symedges[0], "  se0");
  }
#endif

  if (is_deleted_edge(e) || e->in_queue) {
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "returning because of e conditions\n");
    }
#endif
    return;
  }
  se = &e->symedges[0];
  sesym = &e->symedges[1];
  if (len_squared_v2v2_db(se->vert->co, sesym->vert->co) <= cdt->epsilon_squared) {
    BLI_linklist_append_pool(tiny_edge_queue, e, cdt->listpool);
    e->in_queue = true;
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "Queue tiny edge\n");
    }
#endif
  }
}

/* Consider all edges in rot ring around v for possible enqueing as small features .*/
static void maybe_enqueue_small_features(CDT_state *cdt, CDTVert *v, LinkNodePair *tiny_edge_queue)
{
  SymEdge *se, *se_start;

  se = se_start = v->symedge;
  if (!se_start) {
    return;
  }
  do {
    maybe_enqueue_small_feature(cdt, se->edge, tiny_edge_queue);
  } while ((se = se->rot) != se_start);
}

/* Collapse small edges (length <= epsilon) until no more exist.
 */
static void remove_small_features(CDT_state *cdt)
{
  double epsilon = cdt->epsilon;
  LinkNodePair tiny_edge_queue = {NULL, NULL};
  BLI_mempool *pool = cdt->listpool;
  LinkNode *ln;
  CDTEdge *e;
  CDTVert *v_change;
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nREMOVE_SMALL_FEATURES, epsilon=%g\n", epsilon);
  }
#endif

  if (epsilon == 0.0) {
    return;
  }

  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    maybe_enqueue_small_feature(cdt, e, &tiny_edge_queue);
  }

  while (tiny_edge_queue.list != NULL) {
    e = (CDTEdge *)BLI_linklist_pop_pool(&tiny_edge_queue.list, pool);
    if (tiny_edge_queue.list == NULL) {
      tiny_edge_queue.last_node = NULL;
    }
    e->in_queue = false;
    if (is_deleted_edge(e)) {
      continue;
    }
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "collapse tiny edge\n");
      dump_se(&e->symedges[0], "");
    }
#endif
    v_change = collapse_tiny_edge(cdt, e);
    if (v_change) {
      maybe_enqueue_small_features(cdt, v_change, &tiny_edge_queue);
    }
  }
}

/* Remove all non-constraint edges. */
static void remove_non_constraint_edges(CDT_state *cdt)
{
  LinkNode *ln;
  CDTEdge *e;
  SymEdge *se;

  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    se = &e->symedges[0];
    if (!is_deleted_edge(e) && !is_constrained_edge(e)) {
      dissolve_symedge(cdt, 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). */
struct EdgeToSort {
  double len_squared;
  CDTEdge *e;
};

static int edge_to_sort_cmp(const void *a, const void *b)
{
  const struct EdgeToSort *e1 = a;
  const struct EdgeToSort *e2 = b;

  if (e1->len_squared > e2->len_squared) {
    return -1;
  }
  if (e1->len_squared < e2->len_squared) {
    return 1;
  }
  return 0;
}

static void remove_non_constraint_edges_leave_valid_bmesh(CDT_state *cdt)
{
  LinkNode *ln;
  CDTEdge *e;
  SymEdge *se, *se2;
  CDTFace *fleft, *fright;
  bool dissolve;
  size_t nedges;
  int i, ndissolvable;
  const double *co1, *co2;
  struct EdgeToSort *sorted_edges;

  nedges = 0;
  for (ln = cdt->edges; ln; ln = ln->next) {
    nedges++;
  }
  if (nedges == 0) {
    return;
  }
  sorted_edges = BLI_memarena_alloc(cdt->arena, nedges * sizeof(*sorted_edges));
  i = 0;
  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    if (!is_deleted_edge(e) && !is_constrained_edge(e)) {
      sorted_edges[i].e = e;
      co1 = e->symedges[0].vert->co;
      co2 = e->symedges[1].vert->co;
      sorted_edges[i].len_squared = len_squared_v2v2_db(co1, co2);
      i++;
    }
  }
  ndissolvable = i;
  qsort(sorted_edges, ndissolvable, sizeof(*sorted_edges), edge_to_sort_cmp);
  for (i = 0; i < ndissolvable; i++) {
    e = sorted_edges[i].e;
    se = &e->symedges[0];
    dissolve = true;
    if (true /*!edge_touches_frame(e)*/) {
      fleft = se->face;
      fright = sym(se)->face;
      if (fleft != cdt->outer_face && fright != cdt->outer_face &&
          (fleft->input_ids != NULL || fright->input_ids != NULL)) {
        /* 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 (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, se);
    }
  }
}

static void remove_outer_edges_until_constraints(CDT_state *cdt)
{
  LinkNode *fstack = NULL;
  SymEdge *se, *se_start;
  CDTFace *f, *fsym;
  int visit = ++cdt->visit_count;
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "remove_outer_edges_until_constraints\n");
  }
#endif

  cdt->outer_face->visit_index = visit;
  /* Walk around outer face, adding faces on other side of dissolvable edges to stack. */
  se_start = se = cdt->outer_face->symedge;
  do {
    if (!is_constrained_edge(se->edge)) {
      fsym = sym(se)->face;
      if (fsym->visit_index != visit) {
#ifdef DEBUG_CDT
        if (dbg_level > 0) {
          fprintf(stderr, "pushing f=%p from symedge ", fsym);
          dump_se(se, "an outer edge");
        }
#endif
        BLI_linklist_prepend_pool(&fstack, fsym, cdt->listpool);
      }
    }
  } while ((se = se->next) != se_start);

  while (fstack != NULL) {
    LinkNode *to_dissolve = NULL;
    bool dissolvable;
    f = (CDTFace *)BLI_linklist_pop_pool(&fstack, cdt->listpool);
    if (f->visit_index == visit) {
#ifdef DEBUG_CDT
      if (dbg_level > 0) {
        fprintf(stderr, "skipping f=%p, already visited\n", f);
      }
#endif
      continue;
    }
    BLI_assert(f != cdt->outer_face);
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "top of loop, f=%p\n", f);
      dump_se_cycle(f->symedge, "visit", 10000);
      if (dbg_level > 1) {
        dump_cdt(cdt, "cdt at top of loop");
        cdt_draw(cdt, "top of dissolve loop");
      }
    }
#endif
    f->visit_index = visit;
    se_start = se = f->symedge;
    do {
      dissolvable = !is_constrained_edge(se->edge);
#ifdef DEBUG_CDT
      if (dbg_level > 1) {
        dump_se(se, "edge in f");
        fprintf(stderr, "  dissolvable=%d\n", dissolvable);
      }
#endif
      if (dissolvable) {
        fsym = sym(se)->face;
#ifdef DEBUG_CDT
        if (dbg_level > 1) {
          dump_se_cycle(fsym->symedge, "fsym", 10000);
          fprintf(stderr, "  visited=%d\n", fsym->visit_index == visit);
        }
#endif
        if (fsym->visit_index != visit) {
#ifdef DEBUG_CDT
          if (dbg_level > 0) {
            fprintf(stderr, "pushing face %p\n", fsym);
            dump_se_cycle(fsym->symedge, "pushed", 10000);
          }
#endif
          BLI_linklist_prepend_pool(&fstack, fsym, cdt->listpool);
        }
        else {
          BLI_linklist_prepend_pool(&to_dissolve, se, cdt->listpool);
        }
      }
      se = se->next;
    } while (se != se_start);
    while (to_dissolve != NULL) {
      se = (SymEdge *)BLI_linklist_pop_pool(&to_dissolve, cdt->listpool);
      if (se->next != NULL) {
        dissolve_symedge(cdt, se);
      }
    }
  }
}

/**
 * Remove edges and merge faces to get desired output, as per options.
 * \note the cdt cannot be further changed after this.
 */
static void prepare_cdt_for_output(CDT_state *cdt, const CDT_output_type output_type)
{
  CDTFace *f;
  CDTEdge *e;
  LinkNode *ln;

  cdt->output_prepared = true;
  if (cdt->edges == NULL) {
    return;
  }

  /* Make sure all non-deleted faces have a symedge. */
  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    if (!is_deleted_edge(e)) {
      if (e->symedges[0].face->symedge == NULL) {
        e->symedges[0].face->symedge = &e->symedges[0];
      }
      if (e->symedges[1].face->symedge == NULL) {
        e->symedges[1].face->symedge = &e->symedges[1];
      }
    }
  }
#ifdef DEBUG_CDT
  /* All non-deleted faces should have a symedge now. */
  for (ln = cdt->faces; ln; ln = ln->next) {
    f = (CDTFace *)ln->link;
    if (!f->deleted) {
      BLI_assert(f->symedge != NULL);
    }
  }
#else
  UNUSED_VARS(f);
#endif

  if (output_type == CDT_CONSTRAINTS) {
    remove_non_constraint_edges(cdt);
  }
  else if (output_type == CDT_CONSTRAINTS_VALID_BMESH) {
    remove_non_constraint_edges_leave_valid_bmesh(cdt);
  }
  else if (output_type == CDT_INSIDE) {
    remove_outer_edges_until_constraints(cdt);
  }
}

static CDT_result *cdt_get_output(CDT_state *cdt,
                                  const CDT_input *input,
                                  const CDT_output_type output_type)
{
  int i, j, nv, ne, nf, faces_len_total;
  int orig_map_size, orig_map_index;
  int *vert_to_output_map;
  CDT_result *result;
  CDTVert *v;
  LinkNode *lne, *lnf, *ln;
  SymEdge *se, *se_start;
  CDTEdge *e;
  CDTFace *f;
#ifdef DEBUG_CDT
  int dbg_level = 0;

  if (dbg_level > 0) {
    fprintf(stderr, "\nCDT_GET_OUTPUT\n\n");
  }
#endif

  prepare_cdt_for_output(cdt, output_type);

  result = (CDT_result *)MEM_callocN(sizeof(*result), __func__);
  if (cdt->vert_array_len == 0) {
    return result;
  }

#ifdef DEBUG_CDT
  if (dbg_level > 1) {
    dump_cdt(cdt, "cdt to output");
  }
#endif

  /* 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->vert_array.
   * Since merging picked the leftmost-lowermost representative,
   * that is not necessarily the same as the vertex with the lowest original
   * index (i.e., index in cdt->vert_array), so we need two passes:
   * one to get the non-merged-to vertices in vert_to_output_map,
   * and a second to put in the merge targets for merged-to vertices.
   */
  vert_to_output_map = BLI_memarena_alloc(cdt->arena, (size_t)cdt->vert_array_len * sizeof(int *));
  nv = 0;
  for (i = 0; i < cdt->vert_array_len; i++) {
    v = cdt->vert_array[i];
    if (v->merge_to_index == -1) {
      vert_to_output_map[i] = nv;
      nv++;
    }
  }
  if (nv <= 0) {
    return result;
  }
  if (nv < cdt->vert_array_len) {
    for (i = 0; i < input->verts_len; i++) {
      v = cdt->vert_array[i];
      if (v->merge_to_index != -1) {
        add_to_input_ids(&cdt->vert_array[v->merge_to_index]->input_ids, i, cdt);
        vert_to_output_map[i] = vert_to_output_map[v->merge_to_index];
      }
    }
  }

  result->verts_len = nv;
  result->vert_coords = MEM_malloc_arrayN(nv, sizeof(result->vert_coords[0]), __func__);

  /* Make the vertex "orig" map arrays, mapping output verts to lists of input ones. */
  orig_map_size = 0;
  for (i = 0; i < cdt->vert_array_len; i++) {
    if (cdt->vert_array[i]->merge_to_index == -1) {
      orig_map_size += 1 + BLI_linklist_count(cdt->vert_array[i]->input_ids);
    }
  }
  result->verts_orig_len_table = MEM_malloc_arrayN(nv, sizeof(int), __func__);
  result->verts_orig_start_table = MEM_malloc_arrayN(nv, sizeof(int), __func__);
  result->verts_orig = MEM_malloc_arrayN(orig_map_size, sizeof(int), __func__);

  orig_map_index = 0;
  i = 0;
  for (j = 0; j < cdt->vert_array_len; j++) {
    v = cdt->vert_array[j];
    if (v->merge_to_index == -1) {
      result->vert_coords[i][0] = (float)v->co[0];
      result->vert_coords[i][1] = (float)v->co[1];
      result->verts_orig_start_table[i] = orig_map_index;
      if (j < input->verts_len) {
        result->verts_orig[orig_map_index++] = j;
      }
      for (ln = v->input_ids; ln; ln = ln->next) {
        result->verts_orig[orig_map_index++] = POINTER_AS_INT(ln->link);
      }
      result->verts_orig_len_table[i] = orig_map_index - result->verts_orig_start_table[i];
      i++;
    }
  }

  ne = 0;
  orig_map_size = 0;
  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    if (!is_deleted_edge(e)) {
      ne++;
      if (e->input_ids) {
        orig_map_size += BLI_linklist_count(e->input_ids);
      }
    }
  }
  if (ne != 0) {
    result->edges_len = ne;
    result->face_edge_offset = cdt->face_edge_offset;
    result->edges = MEM_malloc_arrayN(ne, sizeof(result->edges[0]), __func__);
    result->edges_orig_len_table = MEM_malloc_arrayN(ne, sizeof(int), __func__);
    result->edges_orig_start_table = MEM_malloc_arrayN(ne, sizeof(int), __func__);
    if (orig_map_size > 0) {
      result->edges_orig = MEM_malloc_arrayN(orig_map_size, sizeof(int), __func__);
    }
    orig_map_index = 0;
    i = 0;
    for (lne = cdt->edges; lne; lne = lne->next) {
      e = (CDTEdge *)lne->link;
      if (!is_deleted_edge(e)) {
        result->edges[i][0] = vert_to_output_map[e->symedges[0].vert->index];
        result->edges[i][1] = vert_to_output_map[e->symedges[1].vert->index];
        result->edges_orig_start_table[i] = orig_map_index;
        for (ln = e->input_ids; ln; ln = ln->next) {
          result->edges_orig[orig_map_index++] = POINTER_AS_INT(ln->link);
        }
        result->edges_orig_len_table[i] = orig_map_index - result->edges_orig_start_table[i];
        i++;
      }
    }
  }

  nf = 0;
  faces_len_total = 0;
  orig_map_size = 0;
  for (ln = cdt->faces; ln; ln = ln->next) {
    f = (CDTFace *)ln->link;
    if (!f->deleted && f != cdt->outer_face) {
      nf++;
      se = se_start = f->symedge;
      BLI_assert(se != NULL);
      do {
        faces_len_total++;
        se = se->next;
      } while (se != se_start);
      if (f->input_ids) {
        orig_map_size += BLI_linklist_count(f->input_ids);
      }
    }
  }

  if (nf != 0) {
    result->faces_len = nf;
    result->faces_len_table = MEM_malloc_arrayN(nf, sizeof(int), __func__);
    result->faces_start_table = MEM_malloc_arrayN(nf, sizeof(int), __func__);
    result->faces = MEM_malloc_arrayN(faces_len_total, sizeof(int), __func__);
    result->faces_orig_len_table = MEM_malloc_arrayN(nf, sizeof(int), __func__);
    result->faces_orig_start_table = MEM_malloc_arrayN(nf, sizeof(int), __func__);
    if (orig_map_size > 0) {
      result->faces_orig = MEM_malloc_arrayN(orig_map_size, sizeof(int), __func__);
    }
    orig_map_index = 0;
    i = 0;
    j = 0;
    for (lnf = cdt->faces; lnf; lnf = lnf->next) {
      f = (CDTFace *)lnf->link;
      if (!f->deleted && f != cdt->outer_face) {
        result->faces_start_table[i] = j;
        se = se_start = f->symedge;
        do {
          result->faces[j++] = vert_to_output_map[se->vert->index];
          se = se->next;
        } while (se != se_start);
        result->faces_len_table[i] = j - result->faces_start_table[i];
        result->faces_orig_start_table[i] = orig_map_index;
        for (ln = f->input_ids; ln; ln = ln->next) {
          result->faces_orig[orig_map_index++] = POINTER_AS_INT(ln->link);
        }
        result->faces_orig_len_table[i] = orig_map_index - result->faces_orig_start_table[i];
        i++;
      }
    }
  }
  return result;
}

/**
 * Calculate the Constrained Delaunay Triangulation of the 2d elements given in \a input.
 *
 * A Delaunay triangulation of a set of vertices is a triangulation where no triangle in the
 * triangulation has a circumcircle that strictly contains another vertex. Delaunay triangulations
 * are avoid long skinny triangles: they maximize the minimum angle of all triangles in the
 * triangulation.
 *
 * A Constrained Delaunay Triangulation adds the requirement that user-provided line segments must
 * appear as edges in the output (perhaps divided into several sub-segments). It is not required
 * that the input edges be non-intersecting: this routine will calculate the intersections. This
 * means that besides triangulating, this routine is also useful for general and robust 2d edge and
 * face intersection.
 *
 * This routine also takes an epsilon parameter in the \a input. Input vertices closer than epsilon
 * will be merged, and we collapse tiny edges (less than epsilon length).
 *
 * The current initial Deluanay triangulation algorithm is the Guibas-Stolfi Divide and Conquer
 * algorithm (see "Primitives for the Manipulation of General Subdivisions and the Computation of
 * Voronoi Diagrams"). and uses Shewchuk's exact predicates to issues where numeric errors cause
 * inconsistent geometric judgments. This is followed by inserting edge constraints (including the
 * edges implied by faces) using the algorithms discussed in "Fully Dynamic Constrained Delaunay
 * Triangulations" by Kallmann, Bieri, and Thalmann.
 *
 * \param input: points to a CDT_input struct which contains the vertices, edges, and faces to be
 * triangulated. \param output_type: specifies which edges to remove after doing the triangulation.
 * \return A pointer to an allocated CDT_result struct, which describes the triangulation in terms
 * of vertices, edges, and faces, and also has tables to map output elements back to input
 * elements. The caller must use BLI_delaunay_2d_cdt_free() on the result when done with it.
 *
 * See the header file BLI_delaunay_2d.h for details of the CDT_input and CDT_result structs and
 * the CDT_output_type enum.
 */
CDT_result *BLI_delaunay_2d_cdt_calc(const CDT_input *input, const CDT_output_type output_type)
{
  int nv = input->verts_len;
  int ne = input->edges_len;
  int nf = input->faces_len;
  int i, iv1, iv2, f, fedge_start, fedge_end, ei;
  CDT_state *cdt;
  CDTVert *v1, *v2;
  CDTEdge *face_edge;
  SymEdge *face_symedge;
  LinkNode *edge_list;
  CDT_result *result;
  const CDT_input *input_orig;
  int *new_edge_map;
  static bool called_exactinit = false;
#ifdef DEBUG_CDT
  int dbg_level = 0;
#endif

  /* The exact orientation and incircle primitives need a one-time initialization of certain
   * constants. */
  if (!called_exactinit) {
    exactinit();
    called_exactinit = true;
  }
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    fprintf(stderr,
            "\n\nCDT CALC, nv=%d, ne=%d, nf=%d, eps=%g\n",
            input->verts_len,
            input->edges_len,
            input->faces_len,
            input->epsilon);
  }
  if (dbg_level == -1) {
    write_cdt_input_to_file(input);
  }
#endif

  if ((nv > 0 && input->vert_coords == NULL) || (ne > 0 && input->edges == NULL) ||
      (nf > 0 && (input->faces == NULL || input->faces_start_table == NULL ||
                  input->faces_len_table == NULL))) {
#ifdef DEBUG_CDT
    fprintf(stderr, "invalid input: unexpected NULL array(s)\n");
#endif
    return NULL;
  }

  input_orig = input;
  input = modify_input_for_near_edge_ends(input, &new_edge_map);
  if (input != input_orig) {
    nv = input->verts_len;
    ne = input->edges_len;
    nf = input->faces_len;
#ifdef DEBUG_CDT
    if (dbg_level > 0) {
      fprintf(stderr, "input modified for near ends; now ne=%d\n", ne);
    }
#endif
  }
  cdt = cdt_init(input);
  initial_triangulation(cdt);
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    validate_cdt(cdt, true, false, false);
    if (dbg_level > 1) {
      cdt_draw(cdt, "after initial triangulation");
    }
  }
#endif

  for (i = 0; i < ne; i++) {
    iv1 = input->edges[i][0];
    iv2 = input->edges[i][1];
    if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
#ifdef DEBUG_CDT
      fprintf(stderr, "edge indices for e%d not valid: v1=%d, v2=%d, nv=%d\n", i, iv1, iv2, nv);
#endif
      continue;
    }
    v1 = cdt->vert_array[iv1];
    v2 = cdt->vert_array[iv2];
    if (v1->merge_to_index != -1) {
      v1 = cdt->vert_array[v1->merge_to_index];
    }
    if (v2->merge_to_index != -1) {
      v2 = cdt->vert_array[v2->merge_to_index];
    }
    if (new_edge_map) {
      ei = new_edge_map[i];
    }
    else {
      ei = i;
    }
    add_edge_constraint(cdt, v1, v2, ei, NULL);
#ifdef DEBUG_CDT
    if (dbg_level > 3) {
      char namebuf[60];
      sprintf(namebuf, "after edge constraint %d = (%d,%d)\n", i, iv1, iv2);
      cdt_draw(cdt, namebuf);
      // dump_cdt(cdt, namebuf);
      validate_cdt(cdt, true, true, false);
    }
#endif
  }

  cdt->face_edge_offset = ne;
  for (f = 0; f < nf; f++) {
    int flen = input->faces_len_table[f];
    int fstart = input->faces_start_table[f];
    if (flen <= 2) {
#ifdef DEBUG_CDT
      fprintf(stderr, "face %d has length %d; ignored\n", f, flen);
#endif
      continue;
    }
    for (i = 0; i < flen; i++) {
      int face_edge_id = cdt->face_edge_offset + fstart + i;
      if (new_edge_map) {
        face_edge_id = new_edge_map[face_edge_id];
      }
      iv1 = input->faces[fstart + i];
      iv2 = input->faces[fstart + ((i + 1) % flen)];
      if (iv1 < 0 || iv1 >= nv || iv2 < 0 || iv2 >= nv) {
#ifdef DEBUG_CDT
        fprintf(stderr, "face indices not valid: f=%d, iv1=%d, iv2=%d, nv=%d\n", f, iv1, iv2, nv);
#endif
        continue;
      }
      v1 = cdt->vert_array[iv1];
      v2 = cdt->vert_array[iv2];
      if (v1->merge_to_index != -1) {
        v1 = cdt->vert_array[v1->merge_to_index];
      }
      if (v2->merge_to_index != -1) {
        v2 = cdt->vert_array[v2->merge_to_index];
      }
      add_edge_constraint(cdt, v1, v2, face_edge_id, &edge_list);
#ifdef DEBUG_CDT
      if (dbg_level > 2) {
        fprintf(stderr, "edges for edge %d:\n", i);
        for (LinkNode *ln = edge_list; ln; ln = ln->next) {
          CDTEdge *cdt_e = (CDTEdge *)ln->link;
          fprintf(stderr,
                  "  (%.2f,%.2f)->(%.2f,%.2f)\n",
                  F2(cdt_e->symedges[0].vert->co),
                  F2(cdt_e->symedges[1].vert->co));
        }
      }
      if (dbg_level > 2) {
        cdt_draw(cdt, "after a face edge");
        if (dbg_level > 3) {
          dump_cdt(cdt, "after a face edge");
        }
        validate_cdt(cdt, true, true, false);
      }
#endif
      if (i == 0) {
        face_edge = (CDTEdge *)edge_list->link;
        face_symedge = &face_edge->symedges[0];
        if (face_symedge->vert != v1) {
          face_symedge = &face_edge->symedges[1];
          BLI_assert(face_symedge->vert == v1);
        }
      }
      BLI_linklist_free_pool(edge_list, NULL, cdt->listpool);
    }
    fedge_start = cdt->face_edge_offset + fstart;
    fedge_end = fedge_start + flen - 1;
    add_face_ids(cdt, face_symedge, f, fedge_start, fedge_end);
  }
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    validate_cdt(cdt, true, true, false);
  }
  if (dbg_level > 1) {
    cdt_draw(cdt, "after adding edges and faces");
    if (dbg_level > 2) {
      dump_cdt(cdt, "after adding edges and faces");
    }
  }
#endif

  if (cdt->epsilon > 0.0) {
    remove_small_features(cdt);
#ifdef DEBUG_CDT
    if (dbg_level > 2) {
      cdt_draw(cdt, "after remove small features\n");
      if (dbg_level > 3) {
        dump_cdt(cdt, "after remove small features\n");
      }
    }
#endif
  }

  result = cdt_get_output(cdt, input, output_type);
#ifdef DEBUG_CDT
  if (dbg_level > 0) {
    cdt_draw(cdt, "final");
  }
#endif

  if (input != input_orig) {
    free_modified_input((CDT_input *)input);
  }
  new_cdt_free(cdt);
  return result;
}

void BLI_delaunay_2d_cdt_free(CDT_result *result)
{
  if (result == NULL) {
    return;
  }
  if (result->vert_coords) {
    MEM_freeN(result->vert_coords);
  }
  if (result->edges) {
    MEM_freeN(result->edges);
  }
  if (result->faces) {
    MEM_freeN(result->faces);
  }
  if (result->faces_start_table) {
    MEM_freeN(result->faces_start_table);
  }
  if (result->faces_len_table) {
    MEM_freeN(result->faces_len_table);
  }
  if (result->verts_orig) {
    MEM_freeN(result->verts_orig);
  }
  if (result->verts_orig_start_table) {
    MEM_freeN(result->verts_orig_start_table);
  }
  if (result->verts_orig_len_table) {
    MEM_freeN(result->verts_orig_len_table);
  }
  if (result->edges_orig) {
    MEM_freeN(result->edges_orig);
  }
  if (result->edges_orig_start_table) {
    MEM_freeN(result->edges_orig_start_table);
  }
  if (result->edges_orig_len_table) {
    MEM_freeN(result->edges_orig_len_table);
  }
  if (result->faces_orig) {
    MEM_freeN(result->faces_orig);
  }
  if (result->faces_orig_start_table) {
    MEM_freeN(result->faces_orig_start_table);
  }
  if (result->faces_orig_len_table) {
    MEM_freeN(result->faces_orig_len_table);
  }
  MEM_freeN(result);
}

#ifdef DEBUG_CDT

ATTU static const char *vertname(const CDTVert *v)
{
  static char vertnamebuf[20];

  sprintf(vertnamebuf, "[%d]", v->index);
  return vertnamebuf;
}

ATTU static const char *sename(const SymEdge *se)
{
  static char senamebuf[20];

  sprintf(senamebuf, "{%x}", (POINTER_AS_UINT(se)) & 0xFFFF);
  return senamebuf;
}

static void dump_v(const CDTVert *v, const char *lab)
{
  fprintf(stderr, "%s%s(%.10f,%.10f)\n", lab, vertname(v), F2(v->co));
}

static void dump_se(const SymEdge *se, const char *lab)
{
  if (se->next) {
    fprintf(stderr,
            "%s%s((%.10f,%.10f)->(%.10f,%.10f))",
            lab,
            vertname(se->vert),
            F2(se->vert->co),
            F2(se->next->vert->co));
    fprintf(stderr, "%s\n", vertname(se->next->vert));
  }
  else {
    fprintf(stderr, "%s%s((%.10f,%.10f)->NULL)\n", lab, vertname(se->vert), F2(se->vert->co));
  }
}

static void dump_se_short(const SymEdge *se, const char *lab)
{
  if (se == NULL) {
    fprintf(stderr, "%sNULL", lab);
  }
  else {
    fprintf(stderr, "%s%s", lab, vertname(se->vert));
    fprintf(stderr, "%s", se->next == NULL ? "[NULL]" : vertname(se->next->vert));
  }
}

static void dump_se_cycle(const SymEdge *se, const char *lab, const int limit)
{
  int count = 0;
  const SymEdge *s = se;
  fprintf(stderr, "%s:\n", lab);
  do {
    dump_se(s, "  ");
    s = s->next;
    count++;
  } while (s != se && count < limit);
  if (count == limit) {
    fprintf(stderr, "  limit hit without cycle!\n");
  }
}

static void dump_id_list(const LinkNode *id_list, const char *lab)
{
  const LinkNode *ln;
  if (!id_list) {
    return;
  }
  fprintf(stderr, "%s", lab);
  for (ln = id_list; ln; ln = ln->next) {
    fprintf(stderr, "%d%c", POINTER_AS_INT(ln->link), ln->next ? ' ' : '\n');
  }
}

static void dump_cross_data(struct CrossData *cd, const char *lab)
{
  fprintf(stderr, "%s", lab);
  if (cd->lambda == 0.0) {
    fprintf(stderr, "v%d", cd->vert->index);
  }
  else {
    fprintf(stderr, "lambda=%.17g", cd->lambda);
  }
  dump_se_short(cd->in, " in=");
  dump_se_short(cd->out, " out=");
  fprintf(stderr, "\n");
}

/* If filter_fn != NULL, only dump vert v its edges when filter_fn(cdt, v, filter_data) is true. */
#  define PL(p) (POINTER_AS_UINT(p) & 0xFFFF)
static void dump_cdt_filtered(const CDT_state *cdt,
                              bool (*filter_fn)(const CDT_state *, int, void *),
                              void *filter_data,
                              const char *lab)
{
  LinkNode *ln;
  CDTVert *v, *vother;
  CDTEdge *e;
  CDTFace *f;
  SymEdge *se;
  int i, cnt;

  fprintf(stderr, "\nCDT %s\n", lab);
  fprintf(stderr, "\nVERTS\n");
  for (i = 0; i < cdt->vert_array_len; i++) {
    if (filter_fn && !filter_fn(cdt, i, filter_data)) {
      continue;
    }
    v = cdt->vert_array[i];
    fprintf(stderr, "%s %x: (%f,%f) symedge=%x", vertname(v), PL(v), F2(v->co), PL(v->symedge));
    if (v->merge_to_index == -1) {
      fprintf(stderr, "\n");
    }
    else {
      fprintf(stderr, " merge to %s\n", vertname(cdt->vert_array[v->merge_to_index]));
      continue;
    }
    dump_id_list(v->input_ids, "  ");
    se = v->symedge;
    cnt = 0;
    if (se) {
      fprintf(stderr, "  edges out:\n");
      do {
        if (se->next == NULL) {
          fprintf(stderr, "    [NULL next/rot symedge, se=%x\n", PL(se));
          break;
        }
        if (se->next->next == NULL) {
          fprintf(stderr, "    [NULL next-next/rot symedge, se=%x\n", PL(se));
          break;
        }
        vother = sym(se)->vert;
        fprintf(stderr, "    %s (e=%x, se=%x)\n", vertname(vother), PL(se->edge), PL(se));
        se = se->rot;
        cnt++;
      } while (se != v->symedge && cnt < 25);
      fprintf(stderr, "\n");
    }
  }
  if (filter_fn) {
    return;
  }
  fprintf(stderr, "\nEDGES\n");
  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    if (e->symedges[0].next == NULL) {
      continue;
    }
    fprintf(stderr, "%x:\n", PL(e));
    for (i = 0; i < 2; i++) {
      se = &e->symedges[i];
      fprintf(stderr,
              "  se[%d] @%x: next=%x, rot=%x, vert=%x [%s] (%.2f,%.2f), edge=%x, face=%x\n",
              i,
              PL(se),
              PL(se->next),
              PL(se->rot),
              PL(se->vert),
              vertname(se->vert),
              F2(se->vert->co),
              PL(se->edge),
              PL(se->face));
    }
    dump_id_list(e->input_ids, "  ");
  }
  fprintf(stderr, "\nFACES\n");
  for (ln = cdt->faces; ln; ln = ln->next) {
    f = (CDTFace *)ln->link;
    if (f->deleted) {
      continue;
    }
    if (f == cdt->outer_face) {
      fprintf(stderr, "%x: outer", PL(f));
    }
    fprintf(stderr, " symedge=%x\n", PL(f->symedge));
    dump_id_list(f->input_ids, "  ");
  }
  fprintf(stderr, "\nOTHER\n");
  fprintf(stderr, "outer_face=%x\n", PL(cdt->outer_face));
  fprintf(
      stderr, "minx=%f, maxx=%f, miny=%f, maxy=%f\n", cdt->minx, cdt->maxx, cdt->miny, cdt->maxy);
  fprintf(stderr, "margin=%f\n", cdt->margin);
}
#  undef PL

static void dump_cdt(const CDT_state *cdt, const char *lab)
{
  dump_cdt_filtered(cdt, NULL, NULL, lab);
}

typedef struct ReachableFilterData {
  int vstart_index;
  int maxdist;
} ReachableFilterData;

/* Stupid algorithm will repeat itself. Don't use for large n. */
static bool reachable_filter(const CDT_state *cdt, int v_index, void *filter_data)
{
  CDTVert *v;
  SymEdge *se;
  ReachableFilterData *rfd_in = (ReachableFilterData *)filter_data;
  ReachableFilterData rfd_next;

  if (v_index == rfd_in->vstart_index) {
    return true;
  }
  if (rfd_in->maxdist <= 0 || v_index < 0 || v_index >= cdt->vert_array_len) {
    return false;
  }
  else {
    v = cdt->vert_array[v_index];
    se = v->symedge;
    if (se != NULL) {
      rfd_next.vstart_index = rfd_in->vstart_index;
      rfd_next.maxdist = rfd_in->maxdist - 1;
      do {
        if (reachable_filter(cdt, se->next->vert->index, &rfd_next)) {
          return true;
        }
        se = se->rot;
      } while (se != v->symedge);
    }
  }
  return false;
}

static void set_min_max(CDT_state *cdt)
{
  int i;
  double minx, maxx, miny, maxy;
  double *co;

  minx = miny = DBL_MAX;
  maxx = maxy = -DBL_MAX;
  for (i = 0; i < cdt->vert_array_len; i++) {
    co = cdt->vert_array[i]->co;
    if (co[0] < minx) {
      minx = co[0];
    }
    if (co[0] > maxx) {
      maxx = co[0];
    }
    if (co[1] < miny) {
      miny = co[1];
    }
    if (co[1] > maxy) {
      maxy = co[1];
    }
  }
  if (minx != DBL_MAX) {
    cdt->minx = minx;
    cdt->miny = miny;
    cdt->maxx = maxx;
    cdt->maxy = maxy;
  }
}

static void dump_cdt_vert_neighborhood(CDT_state *cdt, int v, int maxdist, const char *lab)
{
  ReachableFilterData rfd;
  rfd.vstart_index = v;
  rfd.maxdist = maxdist;
  dump_cdt_filtered(cdt, reachable_filter, &rfd, lab);
}

/*
 * Make an html file with svg in it to display the argument cdt.
 * Mouse-overs will reveal the coordinates of vertices and edges.
 * Constraint edges are drawn thicker than non-constraint edges.
 * The first call creates DRAWFILE; subsequent calls append to it.
 */
#  define DRAWFILE "/tmp/debug_draw.html"
#  define MAX_DRAW_WIDTH 2000
#  define MAX_DRAW_HEIGHT 1400
#  define THIN_LINE 1
#  define THICK_LINE 4
#  define VERT_RADIUS 3
#  define DRAW_VERT_LABELS 1
#  define DRAW_EDGE_LABELS 0

static void cdt_draw_region(
    CDT_state *cdt, const char *lab, double minx, double miny, double maxx, double maxy)
{
  static bool append = false;
  FILE *f = fopen(DRAWFILE, append ? "a" : "w");
  int view_width, view_height;
  double width, height, aspect, scale;
  LinkNode *ln;
  CDTVert *v, *u;
  CDTEdge *e;
  int i, strokew;

  width = maxx - minx;
  height = maxy - miny;
  aspect = height / width;
  view_width = MAX_DRAW_WIDTH;
  view_height = (int)(view_width * aspect);
  if (view_height > MAX_DRAW_HEIGHT) {
    view_height = MAX_DRAW_HEIGHT;
    view_width = (int)(view_height / aspect);
  }
  scale = view_width / width;

#  define SX(x) ((x - minx) * scale)
#  define SY(y) ((maxy - y) * scale)

  if (!f) {
    printf("couldn't open file %s\n", DRAWFILE);
    return;
  }
  fprintf(f, "<div>%s</div>\n<div>\n", lab);
  fprintf(f,
          "<svg version=\"1.1\" "
          "xmlns=\"http://www.w3.org/2000/svg\" "
          "xmlns:xlink=\"http://www.w3.org/1999/xlink\" "
          "xml:space=\"preserve\"\n");
  fprintf(f, "width=\"%d\" height=\"%d\">/n", view_width, view_height);

  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    if (is_deleted_edge(e)) {
      continue;
    }
    u = e->symedges[0].vert;
    v = e->symedges[1].vert;
    strokew = is_constrained_edge(e) ? THICK_LINE : THIN_LINE;
    fprintf(f,
            "<line fill=\"none\" stroke=\"black\" stroke-width=\"%d\" "
            "x1=\"%f\" y1=\"%f\" x2=\"%f\" y2=\"%f\">\n",
            strokew,
            SX(u->co[0]),
            SY(u->co[1]),
            SX(v->co[0]),
            SY(v->co[1]));
    fprintf(f, "  <title>%s", vertname(u));
    fprintf(f, "%s</title>\n", vertname(v));
    fprintf(f, "</line>\n");
#  if DRAW_EDGE_LABELS
    fprintf(f,
            "<text x=\"%f\" y=\"%f\" font-size=\"small\">",
            SX(0.5 * (u->co[0] + v->co[0])),
            SY(0.5 * (u->co[1] + v->co[1])));
    fprintf(f, "%s", vertname(u));
    fprintf(f, "%s", vertname(v));
    fprintf(f, "%s", sename(&e->symedges[0]));
    fprintf(f, "%s</text>\n", sename(&e->symedges[1]));
#  endif
  }
  i = 0;
  for (; i < cdt->vert_array_len; i++) {
    v = cdt->vert_array[i];
    if (v->merge_to_index != -1) {
      continue;
    }
    fprintf(f,
            "<circle fill=\"black\" cx=\"%f\" cy=\"%f\" r=\"%d\">\n",
            SX(v->co[0]),
            SY(v->co[1]),
            VERT_RADIUS);
    fprintf(f, "  <title>%s(%.10f,%.10f)</title>\n", vertname(v), v->co[0], v->co[1]);
    fprintf(f, "</circle>\n");
#  if DRAW_VERT_LABELS
    fprintf(f,
            "<text x=\"%f\" y=\"%f\" font-size=\"small\">%s</text>\n",
            SX(v->co[0]) + VERT_RADIUS,
            SY(v->co[1]) - VERT_RADIUS,
            vertname(v));
#  endif
  }

  fprintf(f, "</svg>\n</div>\n");
  fclose(f);
  append = true;
#  undef SX
#  undef SY
}

static void cdt_draw(CDT_state *cdt, const char *lab)
{
  double draw_margin, minx, maxx, miny, maxy;

  set_min_max(cdt);
  draw_margin = (cdt->maxx - cdt->minx + cdt->maxy - cdt->miny + 1) * 0.05;
  minx = cdt->minx - draw_margin;
  maxx = cdt->maxx + draw_margin;
  miny = cdt->miny - draw_margin;
  maxy = cdt->maxy + draw_margin;
  cdt_draw_region(cdt, lab, minx, miny, maxx, maxy);
}

static void cdt_draw_vertex_region(CDT_state *cdt, int v, double dist, const char *lab)
{
  const double *co = cdt->vert_array[v]->co;
  cdt_draw_region(cdt, lab, co[0] - dist, co[1] - dist, co[0] + dist, co[1] + dist);
}

static void cdt_draw_edge_region(CDT_state *cdt, int v1, int v2, double dist, const char *lab)
{
  const double *co1 = cdt->vert_array[v1]->co;
  const double *co2 = cdt->vert_array[v2]->co;
  double minx, miny, maxx, maxy;

  minx = min_dd(co1[0], co2[0]);
  miny = min_dd(co1[1], co2[1]);
  maxx = max_dd(co1[0], co2[0]);
  maxy = max_dd(co1[1], co2[1]);
  cdt_draw_region(cdt, lab, minx - dist, miny - dist, maxx + dist, maxy + dist);
}

#  define CDTFILE "/tmp/cdtinput.txt"
static void write_cdt_input_to_file(const CDT_input *inp)
{
  int i, j;
  FILE *f = fopen(CDTFILE, "w");

  fprintf(f, "%d %d %d\n", inp->verts_len, inp->edges_len, inp->faces_len);
  for (i = 0; i < inp->verts_len; i++) {
    fprintf(f, "%.17f %.17f\n", inp->vert_coords[i][0], inp->vert_coords[i][1]);
  }
  for (i = 0; i < inp->edges_len; i++) {
    fprintf(f, "%d %d\n", inp->edges[i][0], inp->edges[i][1]);
  }
  for (i = 0; i < inp->faces_len; i++) {
    for (j = 0; j < inp->faces_len_table[i]; j++) {
      fprintf(f, "%d ", inp->faces[j + inp->faces_start_table[i]]);
    }
    fprintf(f, "\n");
  }
  fclose(f);
}

#  ifndef NDEBUG /* Only used in assert. */
/*
 * Is a visible from b: i.e., ab crosses no edge of cdt?
 * If constrained is true, consider only constrained edges as possible crossed edges.
 * In any case, don't count an edge ab itself.
 * Note: this is an expensive test if there are a lot of edges.
 */
static bool is_visible(const CDTVert *a, const CDTVert *b, bool constrained, const CDT_state *cdt)
{
  const LinkNode *ln;
  const CDTEdge *e;
  const SymEdge *se, *senext;
  double lambda, mu;
  int ikind;

  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (const CDTEdge *)ln->link;
    if (is_deleted_edge(e) || is_border_edge(e, cdt)) {
      continue;
    }
    if (constrained && !is_constrained_edge(e)) {
      continue;
    }
    se = (const SymEdge *)&e->symedges[0];
    senext = se->next;
    if ((a == se->vert || a == senext->vert) || b == se->vert || b == se->next->vert) {
      continue;
    }
    ikind = isect_seg_seg_v2_lambda_mu_db(
        a->co, b->co, se->vert->co, senext->vert->co, &lambda, &mu);
    if (ikind != ISECT_LINE_LINE_NONE) {
      if (ikind == ISECT_LINE_LINE_COLINEAR) {
        /* TODO: special test here for overlap. */
        continue;
      }
      /* Allow an intersection very near or at ends, to allow for numerical error. */
      if (lambda > FLT_EPSILON && (1.0 - lambda) > FLT_EPSILON && mu > FLT_EPSILON &&
          (1.0 - mu) > FLT_EPSILON) {
        return false;
      }
    }
  }
  return true;
}
#  endif

#  ifndef NDEBUG /* Only used in assert. */
/*
 * Check that edge ab satisfies constrained delaunay condition:
 * That is, for all non-constraint, non-border edges ab,
 * (1) ab is visible in the constraint graph; and
 * (2) there is a circle through a and b such that any vertex v connected by an edge to a or b
 *     is not inside that circle.
 * The argument 'se' specifies ab by: a is se's vert and b is se->next's vert.
 * Return true if check is OK.
 */
static bool is_delaunay_edge(const SymEdge *se)
{
  int i;
  CDTVert *a, *b, *c;
  const SymEdge *sesym, *curse, *ss;
  bool ok[2];

  if (!is_constrained_edge(se->edge)) {
    return true;
  }
  sesym = sym(se);
  a = se->vert;
  b = se->next->vert;
  /* Try both the triangles adjacent to se's edge for circle. */
  for (i = 0; i < 2; i++) {
    ok[i] = true;
    curse = (i == 0) ? se : sesym;
    a = curse->vert;
    b = curse->next->vert;
    c = curse->next->next->vert;
    for (ss = curse->rot; ss != curse; ss = ss->rot) {
      ok[i] |= incircle(a->co, b->co, c->co, ss->next->vert->co) <= 0.0;
    }
  }
  return ok[0] || ok[1];
}
#  endif

#  ifndef NDEBUG
static bool plausible_non_null_ptr(void *p)
{
  return p > (void *)0x1000;
}
#  endif

static void validate_cdt(CDT_state *cdt,
                         bool check_all_tris,
                         bool check_delaunay,
                         bool check_visibility)
{
  LinkNode *ln;
  int totedges, totfaces, totverts;
  CDTEdge *e;
  SymEdge *se, *sesym, *s;
  CDTVert *v, *v1, *v2, *v3;
  CDTFace *f;
  int i, limit;
  bool isborder;

  if (cdt->output_prepared) {
    return;
  }
  if (cdt->edges == NULL || cdt->edges->next == NULL) {
    return;
  }

  BLI_assert(cdt != NULL);
  totedges = 0;
  for (ln = cdt->edges; ln; ln = ln->next) {
    e = (CDTEdge *)ln->link;
    se = &e->symedges[0];
    sesym = &e->symedges[1];
    if (is_deleted_edge(e)) {
      BLI_assert(se->rot == NULL && sesym->next == NULL && sesym->rot == NULL);
      continue;
    }
    totedges++;
    isborder = is_border_edge(e, cdt);
    BLI_assert(se->vert != sesym->vert);
    BLI_assert(se->edge == sesym->edge && se->edge == e);
    BLI_assert(sym(se) == sesym && sym(sesym) == se);
    for (i = 0; i < 2; i++) {
      se = &e->symedges[i];
      v = se->vert;
      f = se->face;
      BLI_assert(plausible_non_null_ptr(v));
      if (f != NULL) {
        BLI_assert(plausible_non_null_ptr(f));
      }
      BLI_assert(plausible_non_null_ptr(se->next));
      BLI_assert(plausible_non_null_ptr(se->rot));
      if (check_all_tris && se->face != cdt->outer_face) {
        limit = 3;
      }
      else {
        limit = 10000;
      }
      BLI_assert(reachable(se->next, se, limit));
      if (limit == 3) {
        v1 = se->vert;
        v2 = se->next->vert;
        v3 = se->next->next->vert;
        /* The triangle should be positively oriented, but because
         * the insertion of intersection vertices doesn't use exact
         * arithmetic, this may not be true, so allow a little slop. */
        BLI_assert(orient2d(v1->co, v2->co, v3->co) >= -FLT_EPSILON);
        BLI_assert(orient2d(v2->co, v3->co, v1->co) >= -FLT_EPSILON);
        BLI_assert(orient2d(v3->co, v1->co, v2->co) >= -FLT_EPSILON);
      }
      UNUSED_VARS_NDEBUG(limit);
      BLI_assert(se->next->next != se);
      s = se;
      do {
        BLI_assert(prev(s)->next == s);
        BLI_assert(s->rot == sym(prev(s)));
        s = s->next;
      } while (s != se);
    }
    if (check_visibility) {
      BLI_assert(isborder || is_visible(se->vert, se->next->vert, false, cdt));
    }
    if (!isborder && check_delaunay) {
      BLI_assert(is_delaunay_edge(se));
    }
  }
  totverts = 0;
  for (i = 0; i < cdt->vert_array_len; i++) {
    v = cdt->vert_array[i];
    BLI_assert(plausible_non_null_ptr(v));
    if (v->merge_to_index != -1) {
      BLI_assert(v->merge_to_index >= 0 && v->merge_to_index < cdt->vert_array_len);
      continue;
    }
    totverts++;
    BLI_assert(cdt->vert_array_len <= 1 || v->symedge->vert == v);
  }
  totfaces = 0;
  for (ln = cdt->faces; ln; ln = ln->next) {
    f = (CDTFace *)ln->link;
    BLI_assert(plausible_non_null_ptr(f));
    if (f->deleted) {
      continue;
    }
    totfaces++;
    if (f == cdt->outer_face) {
      continue;
    }
  }
  /* Euler's formula for planar graphs. */
  if (check_all_tris && totfaces > 1) {
    BLI_assert(totverts - totedges + totfaces == 2);
  }
}
#endif

/* Jonathan Shewchuk's adaptive predicates, trimmed to those needed here.
 * Permission obtained by private communication from Jonathan to include this code in Blender.
 */

/*
 *  Routines for Arbitrary Precision Floating-point Arithmetic
 *  and Fast Robust Geometric Predicates
 *  (predicates.c)
 *
 *  May 18, 1996
 *
 *  Placed in the public domain by
 *  Jonathan Richard Shewchuk
 *  School of Computer Science
 *  Carnegie Mellon University
 *  5000 Forbes Avenue
 *  Pittsburgh, Pennsylvania  15213-3891
 *  jrs@cs.cmu.edu
 *
 *  This file contains C implementation of algorithms for exact addition
 *    and multiplication of floating-point numbers, and predicates for
 *    robustly performing the orientation and incircle tests used in
 *    computational geometry.  The algorithms and underlying theory are
 *    described in Jonathan Richard Shewchuk.  "Adaptive Precision Floating-
 *    Point Arithmetic and Fast Robust Geometric Predicates."  Technical
 *    Report CMU-CS-96-140, School of Computer Science, Carnegie Mellon
 *    University, Pittsburgh, Pennsylvania, May 1996.  (Submitted to
 *    Discrete & Computational Geometry.)
 *
 *  This file, the paper listed above, and other information are available
 *    from the Web page http://www.cs.cmu.edu/~quake/robust.html .
 *
 *  Using this code:
 *
 *  First, read the short or long version of the paper (from the Web page
 *    above).
 *
 *  Be sure to call exactinit() once, before calling any of the arithmetic
 *    functions or geometric predicates.  Also be sure to turn on the
 *    optimizer when compiling this file.
 *
 * On some machines, the exact arithmetic routines might be defeated by the
 *   use of internal extended precision floating-point registers.  Sometimes
 *   this problem can be fixed by defining certain values to be volatile,
 *   thus forcing them to be stored to memory and rounded off.  This isn't
 *   a great solution, though, as it slows the arithmetic down.
 *
 * To try this out, write "#define INEXACT volatile" below.  Normally,
 *   however, INEXACT should be defined to be nothing.  ("#define INEXACT".)
 */

#define INEXACT /* Nothing */
/* #define INEXACT volatile */

/* Which of the following two methods of finding the absolute values is
 *   fastest is compiler-dependent.  A few compilers can inline and optimize
 *   the fabs() call; but most will incur the overhead of a function call,
 *   which is disastrously slow.  A faster way on IEEE machines might be to
 *   mask the appropriate bit, but that's difficult to do in C.
 */

#define Absolute(a) ((a) >= 0.0 ? (a) : -(a))
/* #define Absolute(a)  fabs(a) */

/* Many of the operations are broken up into two pieces, a main part that
 *   performs an approximate operation, and a "tail" that computes the
 *   roundoff error of that operation.
 *
 * The operations Fast_Two_Sum(), Fast_Two_Diff(), Two_Sum(), Two_Diff(),
 *   Split(), and Two_Product() are all implemented as described in the
 *   reference.  Each of these macros requires certain variables to be
 *   defined in the calling routine.  The variables `bvirt', `c', `abig',
 *   `_i', `_j', `_k', `_l', `_m', and `_n' are declared `INEXACT' because
 *   they store the result of an operation that may incur roundoff error.
 *   The input parameter `x' (or the highest numbered `x_' parameter) must
 *   also be declared `INEXACT'.
 */

#define Fast_Two_Sum_Tail(a, b, x, y) \
  bvirt = x - a; \
  y = b - bvirt

#define Fast_Two_Sum(a, b, x, y) \
  x = (double)(a + b); \
  Fast_Two_Sum_Tail(a, b, x, y)

#define Fast_Two_Diff_Tail(a, b, x, y) \
  bvirt = a - x; \
  y = bvirt - b

#define Fast_Two_Diff(a, b, x, y) \
  x = (double)(a - b); \
  Fast_Two_Diff_Tail(a, b, x, y)

#define Two_Sum_Tail(a, b, x, y) \
  bvirt = (double)(x - a); \
  avirt = x - bvirt; \
  bround = b - bvirt; \
  around = a - avirt; \
  y = around + bround

#define Two_Sum(a, b, x, y) \
  x = (double)(a + b); \
  Two_Sum_Tail(a, b, x, y)

#define Two_Diff_Tail(a, b, x, y) \
  bvirt = (double)(a - x); \
  avirt = x + bvirt; \
  bround = bvirt - b; \
  around = a - avirt; \
  y = around + bround

#define Two_Diff(a, b, x, y) \
  x = (double)(a - b); \
  Two_Diff_Tail(a, b, x, y)

#define Split(a, ahi, alo) \
  c = (double)(splitter * a); \
  abig = (double)(c - a); \
  ahi = c - abig; \
  alo = a - ahi

#define Two_Product_Tail(a, b, x, y) \
  Split(a, ahi, alo); \
  Split(b, bhi, blo); \
  err1 = x - (ahi * bhi); \
  err2 = err1 - (alo * bhi); \
  err3 = err2 - (ahi * blo); \
  y = (alo * blo) - err3

#define Two_Product(a, b, x, y) \
  x = (double)(a * b); \
  Two_Product_Tail(a, b, x, y)

#define Two_Product_Presplit(a, b, bhi, blo, x, y) \
  x = (double)(a * b); \
  Split(a, ahi, alo); \
  err1 = x - (ahi * bhi); \
  err2 = err1 - (alo * bhi); \
  err3 = err2 - (ahi * blo); \
  y = (alo * blo) - err3

#define Square_Tail(a, x, y) \
  Split(a, ahi, alo); \
  err1 = x - (ahi * ahi); \
  err3 = err1 - ((ahi + ahi) * alo); \
  y = (alo * alo) - err3

#define Square(a, x, y) \
  x = (double)(a * a); \
  Square_Tail(a, x, y)

#define Two_One_Sum(a1, a0, b, x2, x1, x0) \
  Two_Sum(a0, b, _i, x0); \
  Two_Sum(a1, _i, x2, x1)

#define Two_One_Diff(a1, a0, b, x2, x1, x0) \
  Two_Diff(a0, b, _i, x0); \
  Two_Sum(a1, _i, x2, x1)

#define Two_Two_Sum(a1, a0, b1, b0, x3, x2, x1, x0) \
  Two_One_Sum(a1, a0, b0, _j, _0, x0); \
  Two_One_Sum(_j, _0, b1, x3, x2, x1)

#define Two_Two_Diff(a1, a0, b1, b0, x3, x2, x1, x0) \
  Two_One_Diff(a1, a0, b0, _j, _0, x0); \
  Two_One_Diff(_j, _0, b1, x3, x2, x1)

static double splitter;  /* = 2^ceiling(p / 2) + 1.  Used to split floats in half. */
static double m_epsilon; /* = 2^(-p).  Used to estimate roundoff errors. */
/* A set of coefficients used to calculate maximum roundoff errors. */
static double resulterrbound;
static double ccwerrboundA, ccwerrboundB, ccwerrboundC;
static double o3derrboundA, o3derrboundB, o3derrboundC;
static double iccerrboundA, iccerrboundB, iccerrboundC;
static double isperrboundA, isperrboundB, isperrboundC;

/*  exactinit()   Initialize the variables used for exact arithmetic.
 *
 *  `epsilon' is the largest power of two such that 1.0 + epsilon = 1.0 in
 *  floating-point arithmetic.  `epsilon' bounds the relative roundoff
 *  error.  It is used for floating-point error analysis.
 *
 *  `splitter' is used to split floating-point numbers into two
 *  half-length significands for exact multiplication.
 *
 *  I imagine that a highly optimizing compiler might be too smart for its
 *  own good, and somehow cause this routine to fail, if it pretends that
 *  floating-point arithmetic is too much like real arithmetic.
 *
 *  Don't change this routine unless you fully understand it.
 */

static void exactinit(void)
{
  double half;
  double check, lastcheck;
  int every_other;

  every_other = 1;
  half = 0.5;
  m_epsilon = 1.0;
  splitter = 1.0;
  check = 1.0;
  /* Repeatedly divide `epsilon' by two until it is too small to add to
   *   one without causing roundoff.  (Also check if the sum is equal to
   *   the previous sum, for machines that round up instead of using exact
   *   rounding.  Not that this library will work on such machines anyway.
   */
  do {
    lastcheck = check;
    m_epsilon *= half;
    if (every_other) {
      splitter *= 2.0;
    }
    every_other = !every_other;
    check = 1.0 + m_epsilon;
  } while ((check != 1.0) && (check != lastcheck));
  splitter += 1.0;

  /* Error bounds for orientation and incircle tests. */
  resulterrbound = (3.0 + 8.0 * m_epsilon) * m_epsilon;
  ccwerrboundA = (3.0 + 16.0 * m_epsilon) * m_epsilon;
  ccwerrboundB = (2.0 + 12.0 * m_epsilon) * m_epsilon;
  ccwerrboundC = (9.0 + 64.0 * m_epsilon) * m_epsilon * m_epsilon;
  o3derrboundA = (7.0 + 56.0 * m_epsilon) * m_epsilon;
  o3derrboundB = (3.0 + 28.0 * m_epsilon) * m_epsilon;
  o3derrboundC = (26.0 + 288.0 * m_epsilon) * m_epsilon * m_epsilon;
  iccerrboundA = (10.0 + 96.0 * m_epsilon) * m_epsilon;
  iccerrboundB = (4.0 + 48.0 * m_epsilon) * m_epsilon;
  iccerrboundC = (44.0 + 576.0 * m_epsilon) * m_epsilon * m_epsilon;
  isperrboundA = (16.0 + 224.0 * m_epsilon) * m_epsilon;
  isperrboundB = (5.0 + 72.0 * m_epsilon) * m_epsilon;
  isperrboundC = (71.0 + 1408.0 * m_epsilon) * m_epsilon * m_epsilon;
}

/*  fast_expansion_sum_zeroelim()   Sum two expansions, eliminating zero
 *                                  components from the output expansion.
 *
 *  Sets h = e + f.  See the long version of my paper for details.
 *
 *  If round-to-even is used (as with IEEE 754), maintains the strongly
 *  non-overlapping property.  (That is, if e is strongly non-overlapping, h
 *  will be also.)  Does NOT maintain the non-overlapping or non-adjacent
 *  properties.
 */

static int fast_expansion_sum_zeroelim(
    int elen, const double *e, int flen, const double *f, double *h) /* h cannot be e or f. */
{
  double Q;
  INEXACT double Qnew;
  INEXACT double hh;
  INEXACT double bvirt;
  double avirt, bround, around;
  int eindex, findex, hindex;
  double enow, fnow;

  enow = e[0];
  fnow = f[0];
  eindex = findex = 0;
  if ((fnow > enow) == (fnow > -enow)) {
    Q = enow;
    enow = e[++eindex];
  }
  else {
    Q = fnow;
    fnow = f[++findex];
  }
  hindex = 0;
  if ((eindex < elen) && (findex < flen)) {
    if ((fnow > enow) == (fnow > -enow)) {
      Fast_Two_Sum(enow, Q, Qnew, hh);
      enow = e[++eindex];
    }
    else {
      Fast_Two_Sum(fnow, Q, Qnew, hh);
      fnow = f[++findex];
    }
    Q = Qnew;
    if (hh != 0.0) {
      h[hindex++] = hh;
    }
    while ((eindex < elen) && (findex < flen)) {
      if ((fnow > enow) == (fnow > -enow)) {
        Two_Sum(Q, enow, Qnew, hh);
        enow = e[++eindex];
      }
      else {
        Two_Sum(Q, fnow, Qnew, hh);
        fnow = f[++findex];
      }
      Q = Qnew;
      if (hh != 0.0) {
        h[hindex++] = hh;
      }
    }
  }
  while (eindex < elen) {
    Two_Sum(Q, enow, Qnew, hh);
    enow = e[++eindex];
    Q = Qnew;
    if (hh != 0.0) {
      h[hindex++] = hh;
    }
  }
  while (findex < flen) {
    Two_Sum(Q, fnow, Qnew, hh);
    fnow = f[++findex];
    Q = Qnew;
    if (hh != 0.0) {
      h[hindex++] = hh;
    }
  }
  if ((Q != 0.0) || (hindex == 0)) {
    h[hindex++] = Q;
  }
  return hindex;
}

/*  scale_expansion_zeroelim()   Multiply an expansion by a scalar,
 *                               eliminating zero components from the
 *                               output expansion.
 *
 *  Sets h = be.  See either version of my paper for details.
 *
 *  Maintains the nonoverlapping property.  If round-to-even is used (as
 *  with IEEE 754), maintains the strongly nonoverlapping and nonadjacent
 *  properties as well.  (That is, if e has one of these properties, so
 *  will h.)
 */

static int scale_expansion_zeroelim(int elen,
                                    const double *e,
                                    double b,
                                    double *h) /* e and h cannot be the same. */
{
  INEXACT double Q, sum;
  double hh;
  INEXACT double product1;
  double product0;
  int eindex, hindex;
  double enow;
  INEXACT double bvirt;
  double avirt, bround, around;
  INEXACT double c;
  INEXACT double abig;
  double ahi, alo, bhi, blo;
  double err1, err2, err3;

  Split(b, bhi, blo);
  Two_Product_Presplit(e[0], b, bhi, blo, Q, hh);
  hindex = 0;
  if (hh != 0) {
    h[hindex++] = hh;
  }
  for (eindex = 1; eindex < elen; eindex++) {
    enow = e[eindex];
    Two_Product_Presplit(enow, b, bhi, blo, product1, product0);
    Two_Sum(Q, product0, sum, hh);
    if (hh != 0) {
      h[hindex++] = hh;
    }
    Fast_Two_Sum(product1, sum, Q, hh);
    if (hh != 0) {
      h[hindex++] = hh;
    }
  }
  if ((Q != 0.0) || (hindex == 0)) {
    h[hindex++] = Q;
  }
  return hindex;
}

/*  estimate()   Produce a one-word estimate of an expansion's value.
 *
 *  See either version of my paper for details.
 */

static double estimate(int elen, const double *e)
{
  double Q;
  int eindex;

  Q = e[0];
  for (eindex = 1; eindex < elen; eindex++) {
    Q += e[eindex];
  }
  return Q;
}

/*  orient2d()   Adaptive exact 2D orientation test.  Robust.
 *
 *               Return a positive value if the points pa, pb, and pc occur
 *               in counterclockwise order; a negative value if they occur
 *               in clockwise order; and zero if they are collinear.  The
 *               result is also a rough approximation of twice the signed
 *               area of the triangle defined by the three points.
 *
 *  This uses exact arithmetic to ensure a correct answer.  The
 *  result returned is the determinant of a matrix.
 *  This determinant is computed adaptively, in the sense that exact
 *  arithmetic is used only to the degree it is needed to ensure that the
 *  returned value has the correct sign.  Hence, orient2d() is usually quite
 *  fast, but will run more slowly when the input points are collinear or
 *  nearly so.
 */

static double orient2dadapt(const double *pa, const double *pb, const double *pc, double detsum)
{
  INEXACT double acx, acy, bcx, bcy;
  double acxtail, acytail, bcxtail, bcytail;
  INEXACT double detleft, detright;
  double detlefttail, detrighttail;
  double det, errbound;
  double B[4], C1[8], C2[12], D[16];
  INEXACT double B3;
  int C1length, C2length, Dlength;
  double u[4];
  INEXACT double u3;
  INEXACT double s1, t1;
  double s0, t0;

  INEXACT double bvirt;
  double avirt, bround, around;
  INEXACT double c;
  INEXACT double abig;
  double ahi, alo, bhi, blo;
  double err1, err2, err3;
  INEXACT double _i, _j;
  double _0;

  acx = (double)(pa[0] - pc[0]);
  bcx = (double)(pb[0] - pc[0]);
  acy = (double)(pa[1] - pc[1]);
  bcy = (double)(pb[1] - pc[1]);

  Two_Product(acx, bcy, detleft, detlefttail);
  Two_Product(acy, bcx, detright, detrighttail);

  Two_Two_Diff(detleft, detlefttail, detright, detrighttail, B3, B[2], B[1], B[0]);
  B[3] = B3;

  det = estimate(4, B);
  errbound = ccwerrboundB * detsum;
  if ((det >= errbound) || (-det >= errbound)) {
    return det;
  }

  Two_Diff_Tail(pa[0], pc[0], acx, acxtail);
  Two_Diff_Tail(pb[0], pc[0], bcx, bcxtail);
  Two_Diff_Tail(pa[1], pc[1], acy, acytail);
  Two_Diff_Tail(pb[1], pc[1], bcy, bcytail);

  if ((acxtail == 0.0) && (acytail == 0.0) && (bcxtail == 0.0) && (bcytail == 0.0)) {
    return det;
  }

  errbound = ccwerrboundC * detsum + resulterrbound * Absolute(det);
  det += (acx * bcytail + bcy * acxtail) - (acy * bcxtail + bcx * acytail);
  if ((det >= errbound) || (-det >= errbound)) {
    return det;
  }

  Two_Product(acxtail, bcy, s1, s0);
  Two_Product(acytail, bcx, t1, t0);
  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
  u[3] = u3;
  C1length = fast_expansion_sum_zeroelim(4, B, 4, u, C1);

  Two_Product(acx, bcytail, s1, s0);
  Two_Product(acy, bcxtail, t1, t0);
  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
  u[3] = u3;
  C2length = fast_expansion_sum_zeroelim(C1length, C1, 4, u, C2);

  Two_Product(acxtail, bcytail, s1, s0);
  Two_Product(acytail, bcxtail, t1, t0);
  Two_Two_Diff(s1, s0, t1, t0, u3, u[2], u[1], u[0]);
  u[3] = u3;
  Dlength = fast_expansion_sum_zeroelim(C2length, C2, 4, u, D);

  return (D[Dlength - 1]);
}

static double orient2d(const double *pa, const double *pb, const double *pc)
{
  double detleft, detright, det;
  double detsum, errbound;

  detleft = (pa[0] - pc[0]) * (pb[1] - pc[1]);
  detright = (pa[1] - pc[1]) * (pb[0] - pc[0]);
  det = detleft - detright;

  if (detleft > 0.0) {
    if (detright <= 0.0) {
      return det;
    }
    detsum = detleft + detright;
  }
  else if (detleft < 0.0) {
    if (detright >= 0.0) {
      return det;
    }
    detsum = -detleft - detright;
  }
  else {
    return det;
  }

  errbound = ccwerrboundA * detsum;
  if ((det >= errbound) || (-det >= errbound)) {
    return det;
  }

  return orient2dadapt(pa, pb, pc, detsum);
}

/*  incircle()   Adaptive exact 2D incircle test.  Robust.
 *
 *               Return a positive value if the point pd lies inside the
 *               circle passing through pa, pb, and pc; a negative value if
 *               it lies outside; and zero if the four points are cocircular.
 *               The points pa, pb, and pc must be in counterclockwise
 *               order, or the sign of the result will be reversed.
 *
 *  This uses exact arithmetic to ensure a correct answer.
 *  The result returned is the determinant of a matrix.
 *  This determinant is computed adaptively, in the sense that exact
 *  arithmetic is used only to the degree it is needed to ensure that the
 *  returned value has the correct sign.  Hence, incircle() is usually quite
 *  fast, but will run more slowly when the input points are cocircular or
 *  nearly so.
 */

static double incircleadapt(
    const double *pa, const double *pb, const double *pc, const double *pd, double permanent)
{
  INEXACT double adx, bdx, cdx, ady, bdy, cdy;
  double det, errbound;

  INEXACT double bdxcdy1, cdxbdy1, cdxady1, adxcdy1, adxbdy1, bdxady1;
  double bdxcdy0, cdxbdy0, cdxady0, adxcdy0, adxbdy0, bdxady0;
  double bc[4], ca[4], ab[4];
  INEXACT double bc3, ca3, ab3;
  double axbc[8], axxbc[16], aybc[8], ayybc[16], adet[32];
  int axbclen, axxbclen, aybclen, ayybclen, alen;
  double bxca[8], bxxca[16], byca[8], byyca[16], bdet[32];
  int bxcalen, bxxcalen, bycalen, byycalen, blen;
  double cxab[8], cxxab[16], cyab[8], cyyab[16], cdet[32];
  int cxablen, cxxablen, cyablen, cyyablen, clen;
  double abdet[64];
  int ablen;
  double fin1[1152], fin2[1152];
  double *finnow, *finother, *finswap;
  int finlength;

  double adxtail, bdxtail, cdxtail, adytail, bdytail, cdytail;
  INEXACT double adxadx1, adyady1, bdxbdx1, bdybdy1, cdxcdx1, cdycdy1;
  double adxadx0, adyady0, bdxbdx0, bdybdy0, cdxcdx0, cdycdy0;
  double aa[4], bb[4], cc[4];
  INEXACT double aa3, bb3, cc3;
  INEXACT double ti1, tj1;
  double ti0, tj0;
  double u[4], v[4];
  INEXACT double u3, v3;
  double temp8[8], temp16a[16], temp16b[16], temp16c[16];
  double temp32a[32], temp32b[32], temp48[48], temp64[64];
  int temp8len, temp16alen, temp16blen, temp16clen;
  int temp32alen, temp32blen, temp48len, temp64len;
  double axtbb[8], axtcc[8], aytbb[8], aytcc[8];
  int axtbblen, axtcclen, aytbblen, aytcclen;
  double bxtaa[8], bxtcc[8], bytaa[8], bytcc[8];
  int bxtaalen, bxtcclen, bytaalen, bytcclen;
  double cxtaa[8], cxtbb[8], cytaa[8], cytbb[8];
  int cxtaalen, cxtbblen, cytaalen, cytbblen;
  double axtbc[8], aytbc[8], bxtca[8], bytca[8], cxtab[8], cytab[8];
  int axtbclen, aytbclen, bxtcalen, bytcalen, cxtablen, cytablen;
  double axtbct[16], aytbct[16], bxtcat[16], bytcat[16], cxtabt[16], cytabt[16];
  int axtbctlen, aytbctlen, bxtcatlen, bytcatlen, cxtabtlen, cytabtlen;
  double axtbctt[8], aytbctt[8], bxtcatt[8];
  double bytcatt[8], cxtabtt[8], cytabtt[8];
  int axtbcttlen, aytbcttlen, bxtcattlen, bytcattlen, cxtabttlen, cytabttlen;
  double abt[8], bct[8], cat[8];
  int abtlen, bctlen, catlen;
  double abtt[4], bctt[4], catt[4];
  int abttlen, bcttlen, cattlen;
  INEXACT double abtt3, bctt3, catt3;
  double negate;

  INEXACT double bvirt;
  double avirt, bround, around;
  INEXACT double c;
  INEXACT double abig;
  double ahi, alo, bhi, blo;
  double err1, err2, err3;
  INEXACT double _i, _j;
  double _0;

  adx = (double)(pa[0] - pd[0]);
  bdx = (double)(pb[0] - pd[0]);
  cdx = (double)(pc[0] - pd[0]);
  ady = (double)(pa[1] - pd[1]);
  bdy = (double)(pb[1] - pd[1]);
  cdy = (double)(pc[1] - pd[1]);

  Two_Product(bdx, cdy, bdxcdy1, bdxcdy0);
  Two_Product(cdx, bdy, cdxbdy1, cdxbdy0);
  Two_Two_Diff(bdxcdy1, bdxcdy0, cdxbdy1, cdxbdy0, bc3, bc[2], bc[1], bc[0]);
  bc[3] = bc3;
  axbclen = scale_expansion_zeroelim(4, bc, adx, axbc);
  axxbclen = scale_expansion_zeroelim(axbclen, axbc, adx, axxbc);
  aybclen = scale_expansion_zeroelim(4, bc, ady, aybc);
  ayybclen = scale_expansion_zeroelim(aybclen, aybc, ady, ayybc);
  alen = fast_expansion_sum_zeroelim(axxbclen, axxbc, ayybclen, ayybc, adet);

  Two_Product(cdx, ady, cdxady1, cdxady0);
  Two_Product(adx, cdy, adxcdy1, adxcdy0);
  Two_Two_Diff(cdxady1, cdxady0, adxcdy1, adxcdy0, ca3, ca[2], ca[1], ca[0]);
  ca[3] = ca3;
  bxcalen = scale_expansion_zeroelim(4, ca, bdx, bxca);
  bxxcalen = scale_expansion_zeroelim(bxcalen, bxca, bdx, bxxca);
  bycalen = scale_expansion_zeroelim(4, ca, bdy, byca);
  byycalen = scale_expansion_zeroelim(bycalen, byca, bdy, byyca);
  blen = fast_expansion_sum_zeroelim(bxxcalen, bxxca, byycalen, byyca, bdet);

  Two_Product(adx, bdy, adxbdy1, adxbdy0);
  Two_Product(bdx, ady, bdxady1, bdxady0);
  Two_Two_Diff(adxbdy1, adxbdy0, bdxady1, bdxady0, ab3, ab[2], ab[1], ab[0]);
  ab[3] = ab3;
  cxablen = scale_expansion_zeroelim(4, ab, cdx, cxab);
  cxxablen = scale_expansion_zeroelim(cxablen, cxab, cdx, cxxab);
  cyablen = scale_expansion_zeroelim(4, ab, cdy, cyab);
  cyyablen = scale_expansion_zeroelim(cyablen, cyab, cdy, cyyab);
  clen = fast_expansion_sum_zeroelim(cxxablen, cxxab, cyyablen, cyyab, cdet);

  ablen = fast_expansion_sum_zeroelim(alen, adet, blen, bdet, abdet);
  finlength = fast_expansion_sum_zeroelim(ablen, abdet, clen, cdet, fin1);

  det = estimate(finlength, fin1);
  errbound = iccerrboundB * permanent;
  if ((det >= errbound) || (-det >= errbound)) {
    return det;
  }

  Two_Diff_Tail(pa[0], pd[0], adx, adxtail);
  Two_Diff_Tail(pa[1], pd[1], ady, adytail);
  Two_Diff_Tail(pb[0], pd[0], bdx, bdxtail);
  Two_Diff_Tail(pb[1], pd[1], bdy, bdytail);
  Two_Diff_Tail(pc[0], pd[0], cdx, cdxtail);
  Two_Diff_Tail(pc[1], pd[1], cdy, cdytail);
  if ((adxtail == 0.0) && (bdxtail == 0.0) && (cdxtail == 0.0) && (adytail == 0.0) &&
      (bdytail == 0.0) && (cdytail == 0.0)) {
    return det;
  }

  errbound = iccerrboundC * permanent + resulterrbound * Absolute(det);
  det += ((adx * adx + ady * ady) *
              ((bdx * cdytail + cdy * bdxtail) - (bdy * cdxtail + cdx * bdytail)) +
          2.0 * (adx * adxtail + ady * adytail) * (bdx * cdy - bdy * cdx)) +
         ((bdx * bdx + bdy * bdy) *
              ((cdx * adytail + ady * cdxtail) - (cdy * adxtail + adx * cdytail)) +
          2.0 * (bdx * bdxtail + bdy * bdytail) * (cdx * ady - cdy * adx)) +
         ((cdx * cdx + cdy * cdy) *
              ((adx * bdytail + bdy * adxtail) - (ady * bdxtail + bdx * adytail)) +
          2.0 * (cdx * cdxtail + cdy * cdytail) * (adx * bdy - ady * bdx));
  if ((det >= errbound) || (-det >= errbound)) {
    return det;
  }

  finnow = fin1;
  finother = fin2;

  if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0)) {
    Square(adx, adxadx1, adxadx0);
    Square(ady, adyady1, adyady0);
    Two_Two_Sum(adxadx1, adxadx0, adyady1, adyady0, aa3, aa[2], aa[1], aa[0]);
    aa[3] = aa3;
  }
  if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0)) {
    Square(bdx, bdxbdx1, bdxbdx0);
    Square(bdy, bdybdy1, bdybdy0);
    Two_Two_Sum(bdxbdx1, bdxbdx0, bdybdy1, bdybdy0, bb3, bb[2], bb[1], bb[0]);
    bb[3] = bb3;
  }
  if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0)) {
    Square(cdx, cdxcdx1, cdxcdx0);
    Square(cdy, cdycdy1, cdycdy0);
    Two_Two_Sum(cdxcdx1, cdxcdx0, cdycdy1, cdycdy0, cc3, cc[2], cc[1], cc[0]);
    cc[3] = cc3;
  }

  if (adxtail != 0.0) {
    axtbclen = scale_expansion_zeroelim(4, bc, adxtail, axtbc);
    temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, 2.0 * adx, temp16a);

    axtcclen = scale_expansion_zeroelim(4, cc, adxtail, axtcc);
    temp16blen = scale_expansion_zeroelim(axtcclen, axtcc, bdy, temp16b);

    axtbblen = scale_expansion_zeroelim(4, bb, adxtail, axtbb);
    temp16clen = scale_expansion_zeroelim(axtbblen, axtbb, -cdy, temp16c);

    temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
    temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
    finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
    finswap = finnow;
    finnow = finother;
    finother = finswap;
  }
  if (adytail != 0.0) {
    aytbclen = scale_expansion_zeroelim(4, bc, adytail, aytbc);
    temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, 2.0 * ady, temp16a);

    aytbblen = scale_expansion_zeroelim(4, bb, adytail, aytbb);
    temp16blen = scale_expansion_zeroelim(aytbblen, aytbb, cdx, temp16b);

    aytcclen = scale_expansion_zeroelim(4, cc, adytail, aytcc);
    temp16clen = scale_expansion_zeroelim(aytcclen, aytcc, -bdx, temp16c);

    temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
    temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
    finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
    finswap = finnow;
    finnow = finother;
    finother = finswap;
  }
  if (bdxtail != 0.0) {
    bxtcalen = scale_expansion_zeroelim(4, ca, bdxtail, bxtca);
    temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, 2.0 * bdx, temp16a);

    bxtaalen = scale_expansion_zeroelim(4, aa, bdxtail, bxtaa);
    temp16blen = scale_expansion_zeroelim(bxtaalen, bxtaa, cdy, temp16b);

    bxtcclen = scale_expansion_zeroelim(4, cc, bdxtail, bxtcc);
    temp16clen = scale_expansion_zeroelim(bxtcclen, bxtcc, -ady, temp16c);

    temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
    temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
    finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
    finswap = finnow;
    finnow = finother;
    finother = finswap;
  }
  if (bdytail != 0.0) {
    bytcalen = scale_expansion_zeroelim(4, ca, bdytail, bytca);
    temp16alen = scale_expansion_zeroelim(bytcalen, bytca, 2.0 * bdy, temp16a);

    bytcclen = scale_expansion_zeroelim(4, cc, bdytail, bytcc);
    temp16blen = scale_expansion_zeroelim(bytcclen, bytcc, adx, temp16b);

    bytaalen = scale_expansion_zeroelim(4, aa, bdytail, bytaa);
    temp16clen = scale_expansion_zeroelim(bytaalen, bytaa, -cdx, temp16c);

    temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
    temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
    finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
    finswap = finnow;
    finnow = finother;
    finother = finswap;
  }
  if (cdxtail != 0.0) {
    cxtablen = scale_expansion_zeroelim(4, ab, cdxtail, cxtab);
    temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, 2.0 * cdx, temp16a);

    cxtbblen = scale_expansion_zeroelim(4, bb, cdxtail, cxtbb);
    temp16blen = scale_expansion_zeroelim(cxtbblen, cxtbb, ady, temp16b);

    cxtaalen = scale_expansion_zeroelim(4, aa, cdxtail, cxtaa);
    temp16clen = scale_expansion_zeroelim(cxtaalen, cxtaa, -bdy, temp16c);

    temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
    temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
    finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
    finswap = finnow;
    finnow = finother;
    finother = finswap;
  }
  if (cdytail != 0.0) {
    cytablen = scale_expansion_zeroelim(4, ab, cdytail, cytab);
    temp16alen = scale_expansion_zeroelim(cytablen, cytab, 2.0 * cdy, temp16a);

    cytaalen = scale_expansion_zeroelim(4, aa, cdytail, cytaa);
    temp16blen = scale_expansion_zeroelim(cytaalen, cytaa, bdx, temp16b);

    cytbblen = scale_expansion_zeroelim(4, bb, cdytail, cytbb);
    temp16clen = scale_expansion_zeroelim(cytbblen, cytbb, -adx, temp16c);

    temp32alen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32a);
    temp48len = fast_expansion_sum_zeroelim(temp16clen, temp16c, temp32alen, temp32a, temp48);
    finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
    finswap = finnow;
    finnow = finother;
    finother = finswap;
  }

  if ((adxtail != 0.0) || (adytail != 0.0)) {
    if ((bdxtail != 0.0) || (bdytail != 0.0) || (cdxtail != 0.0) || (cdytail != 0.0)) {
      Two_Product(bdxtail, cdy, ti1, ti0);
      Two_Product(bdx, cdytail, tj1, tj0);
      Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
      u[3] = u3;
      negate = -bdy;
      Two_Product(cdxtail, negate, ti1, ti0);
      negate = -bdytail;
      Two_Product(cdx, negate, tj1, tj0);
      Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
      v[3] = v3;
      bctlen = fast_expansion_sum_zeroelim(4, u, 4, v, bct);

      Two_Product(bdxtail, cdytail, ti1, ti0);
      Two_Product(cdxtail, bdytail, tj1, tj0);
      Two_Two_Diff(ti1, ti0, tj1, tj0, bctt3, bctt[2], bctt[1], bctt[0]);
      bctt[3] = bctt3;
      bcttlen = 4;
    }
    else {
      bct[0] = 0.0;
      bctlen = 1;
      bctt[0] = 0.0;
      bcttlen = 1;
    }

    if (adxtail != 0.0) {
      temp16alen = scale_expansion_zeroelim(axtbclen, axtbc, adxtail, temp16a);
      axtbctlen = scale_expansion_zeroelim(bctlen, bct, adxtail, axtbct);
      temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, 2.0 * adx, temp32a);
      temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
      if (bdytail != 0.0) {
        temp8len = scale_expansion_zeroelim(4, cc, adxtail, temp8);
        temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a);
        finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
        finswap = finnow;
        finnow = finother;
        finother = finswap;
      }
      if (cdytail != 0.0) {
        temp8len = scale_expansion_zeroelim(4, bb, -adxtail, temp8);
        temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a);
        finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
        finswap = finnow;
        finnow = finother;
        finother = finswap;
      }

      temp32alen = scale_expansion_zeroelim(axtbctlen, axtbct, adxtail, temp32a);
      axtbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adxtail, axtbctt);
      temp16alen = scale_expansion_zeroelim(axtbcttlen, axtbctt, 2.0 * adx, temp16a);
      temp16blen = scale_expansion_zeroelim(axtbcttlen, axtbctt, adxtail, temp16b);
      temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
      temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
    }
    if (adytail != 0.0) {
      temp16alen = scale_expansion_zeroelim(aytbclen, aytbc, adytail, temp16a);
      aytbctlen = scale_expansion_zeroelim(bctlen, bct, adytail, aytbct);
      temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, 2.0 * ady, temp32a);
      temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;

      temp32alen = scale_expansion_zeroelim(aytbctlen, aytbct, adytail, temp32a);
      aytbcttlen = scale_expansion_zeroelim(bcttlen, bctt, adytail, aytbctt);
      temp16alen = scale_expansion_zeroelim(aytbcttlen, aytbctt, 2.0 * ady, temp16a);
      temp16blen = scale_expansion_zeroelim(aytbcttlen, aytbctt, adytail, temp16b);
      temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
      temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
    }
  }
  if ((bdxtail != 0.0) || (bdytail != 0.0)) {
    if ((cdxtail != 0.0) || (cdytail != 0.0) || (adxtail != 0.0) || (adytail != 0.0)) {
      Two_Product(cdxtail, ady, ti1, ti0);
      Two_Product(cdx, adytail, tj1, tj0);
      Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
      u[3] = u3;
      negate = -cdy;
      Two_Product(adxtail, negate, ti1, ti0);
      negate = -cdytail;
      Two_Product(adx, negate, tj1, tj0);
      Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
      v[3] = v3;
      catlen = fast_expansion_sum_zeroelim(4, u, 4, v, cat);

      Two_Product(cdxtail, adytail, ti1, ti0);
      Two_Product(adxtail, cdytail, tj1, tj0);
      Two_Two_Diff(ti1, ti0, tj1, tj0, catt3, catt[2], catt[1], catt[0]);
      catt[3] = catt3;
      cattlen = 4;
    }
    else {
      cat[0] = 0.0;
      catlen = 1;
      catt[0] = 0.0;
      cattlen = 1;
    }

    if (bdxtail != 0.0) {
      temp16alen = scale_expansion_zeroelim(bxtcalen, bxtca, bdxtail, temp16a);
      bxtcatlen = scale_expansion_zeroelim(catlen, cat, bdxtail, bxtcat);
      temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, 2.0 * bdx, temp32a);
      temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
      if (cdytail != 0.0) {
        temp8len = scale_expansion_zeroelim(4, aa, bdxtail, temp8);
        temp16alen = scale_expansion_zeroelim(temp8len, temp8, cdytail, temp16a);
        finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
        finswap = finnow;
        finnow = finother;
        finother = finswap;
      }
      if (adytail != 0.0) {
        temp8len = scale_expansion_zeroelim(4, cc, -bdxtail, temp8);
        temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a);
        finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
        finswap = finnow;
        finnow = finother;
        finother = finswap;
      }

      temp32alen = scale_expansion_zeroelim(bxtcatlen, bxtcat, bdxtail, temp32a);
      bxtcattlen = scale_expansion_zeroelim(cattlen, catt, bdxtail, bxtcatt);
      temp16alen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, 2.0 * bdx, temp16a);
      temp16blen = scale_expansion_zeroelim(bxtcattlen, bxtcatt, bdxtail, temp16b);
      temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
      temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
    }
    if (bdytail != 0.0) {
      temp16alen = scale_expansion_zeroelim(bytcalen, bytca, bdytail, temp16a);
      bytcatlen = scale_expansion_zeroelim(catlen, cat, bdytail, bytcat);
      temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, 2.0 * bdy, temp32a);
      temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;

      temp32alen = scale_expansion_zeroelim(bytcatlen, bytcat, bdytail, temp32a);
      bytcattlen = scale_expansion_zeroelim(cattlen, catt, bdytail, bytcatt);
      temp16alen = scale_expansion_zeroelim(bytcattlen, bytcatt, 2.0 * bdy, temp16a);
      temp16blen = scale_expansion_zeroelim(bytcattlen, bytcatt, bdytail, temp16b);
      temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
      temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
    }
  }
  if ((cdxtail != 0.0) || (cdytail != 0.0)) {
    if ((adxtail != 0.0) || (adytail != 0.0) || (bdxtail != 0.0) || (bdytail != 0.0)) {
      Two_Product(adxtail, bdy, ti1, ti0);
      Two_Product(adx, bdytail, tj1, tj0);
      Two_Two_Sum(ti1, ti0, tj1, tj0, u3, u[2], u[1], u[0]);
      u[3] = u3;
      negate = -ady;
      Two_Product(bdxtail, negate, ti1, ti0);
      negate = -adytail;
      Two_Product(bdx, negate, tj1, tj0);
      Two_Two_Sum(ti1, ti0, tj1, tj0, v3, v[2], v[1], v[0]);
      v[3] = v3;
      abtlen = fast_expansion_sum_zeroelim(4, u, 4, v, abt);

      Two_Product(adxtail, bdytail, ti1, ti0);
      Two_Product(bdxtail, adytail, tj1, tj0);
      Two_Two_Diff(ti1, ti0, tj1, tj0, abtt3, abtt[2], abtt[1], abtt[0]);
      abtt[3] = abtt3;
      abttlen = 4;
    }
    else {
      abt[0] = 0.0;
      abtlen = 1;
      abtt[0] = 0.0;
      abttlen = 1;
    }

    if (cdxtail != 0.0) {
      temp16alen = scale_expansion_zeroelim(cxtablen, cxtab, cdxtail, temp16a);
      cxtabtlen = scale_expansion_zeroelim(abtlen, abt, cdxtail, cxtabt);
      temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, 2.0 * cdx, temp32a);
      temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
      if (adytail != 0.0) {
        temp8len = scale_expansion_zeroelim(4, bb, cdxtail, temp8);
        temp16alen = scale_expansion_zeroelim(temp8len, temp8, adytail, temp16a);
        finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
        finswap = finnow;
        finnow = finother;
        finother = finswap;
      }
      if (bdytail != 0.0) {
        temp8len = scale_expansion_zeroelim(4, aa, -cdxtail, temp8);
        temp16alen = scale_expansion_zeroelim(temp8len, temp8, bdytail, temp16a);
        finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp16alen, temp16a, finother);
        finswap = finnow;
        finnow = finother;
        finother = finswap;
      }

      temp32alen = scale_expansion_zeroelim(cxtabtlen, cxtabt, cdxtail, temp32a);
      cxtabttlen = scale_expansion_zeroelim(abttlen, abtt, cdxtail, cxtabtt);
      temp16alen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, 2.0 * cdx, temp16a);
      temp16blen = scale_expansion_zeroelim(cxtabttlen, cxtabtt, cdxtail, temp16b);
      temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
      temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
    }
    if (cdytail != 0.0) {
      temp16alen = scale_expansion_zeroelim(cytablen, cytab, cdytail, temp16a);
      cytabtlen = scale_expansion_zeroelim(abtlen, abt, cdytail, cytabt);
      temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, 2.0 * cdy, temp32a);
      temp48len = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp32alen, temp32a, temp48);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp48len, temp48, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;

      temp32alen = scale_expansion_zeroelim(cytabtlen, cytabt, cdytail, temp32a);
      cytabttlen = scale_expansion_zeroelim(abttlen, abtt, cdytail, cytabtt);
      temp16alen = scale_expansion_zeroelim(cytabttlen, cytabtt, 2.0 * cdy, temp16a);
      temp16blen = scale_expansion_zeroelim(cytabttlen, cytabtt, cdytail, temp16b);
      temp32blen = fast_expansion_sum_zeroelim(temp16alen, temp16a, temp16blen, temp16b, temp32b);
      temp64len = fast_expansion_sum_zeroelim(temp32alen, temp32a, temp32blen, temp32b, temp64);
      finlength = fast_expansion_sum_zeroelim(finlength, finnow, temp64len, temp64, finother);
      finswap = finnow;
      finnow = finother;
      finother = finswap;
    }
  }

  return finnow[finlength - 1];
}

static double incircle(const double *pa, const double *pb, const double *pc, const double *pd)
{
  double adx, bdx, cdx, ady, bdy, cdy;
  double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
  double alift, blift, clift;
  double det;
  double permanent, errbound;

  adx = pa[0] - pd[0];
  bdx = pb[0] - pd[0];
  cdx = pc[0] - pd[0];
  ady = pa[1] - pd[1];
  bdy = pb[1] - pd[1];
  cdy = pc[1] - pd[1];

  bdxcdy = bdx * cdy;
  cdxbdy = cdx * bdy;
  alift = adx * adx + ady * ady;

  cdxady = cdx * ady;
  adxcdy = adx * cdy;
  blift = bdx * bdx + bdy * bdy;

  adxbdy = adx * bdy;
  bdxady = bdx * ady;
  clift = cdx * cdx + cdy * cdy;

  det = alift * (bdxcdy - cdxbdy) + blift * (cdxady - adxcdy) + clift * (adxbdy - bdxady);

  permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * alift +
              (Absolute(cdxady) + Absolute(adxcdy)) * blift +
              (Absolute(adxbdy) + Absolute(bdxady)) * clift;
  errbound = iccerrboundA * permanent;
  if ((det > errbound) || (-det > errbound)) {
    return det;
  }

  return incircleadapt(pa, pb, pc, pd, permanent);
}