Constrained Delaunay Triangulation from a Paper

I needed to implement Constrained Delaunay Triangulation for some side project and I decided that I will try implementing it directly from a paper. I did that in the past with other algorithms but now I decided that I will try to describe how that process looked like.

Introduction

What is the Delaunay triangulation? And first of all – what is triangulation?

Triangulation is a process of dividing a polygon into triangles. In the past I wrote an article about an algorithm called ear cut – Triangulation of Polygons.

And now, what is the Delaunay triangulation and why do we need more than one? It's a kind of triangulation which takes a list of points on a plane and creates a set of triangles which connect those points. No point is inside another triangle, and, more importantly, no point is inside the circumcircle of any triangle (circumcircle is a circle that passes through all three vertices of a triangle). So this is different from the ear cut algorithm where we generally want to slice a polygon into triangles. Here we are dealing with points spread on a plane.

Delaunay generally produces nicer looking triangulation (what does nicer mean? well, if you look at Delaunay and non-Delaunay triangulation – the former is likely going to intuitively look better). On the other hand, many other algorithms, including ear cut, can produce long and thin triangles (since they are optimizing for other properties such as how fast they can get rid of a point and move on).

Delaunay triangulation is mathematically proved to maximize the minimum angle of all the angles of the triangles in the triangulation – which intuitively feel like a property that can make it appear better. At the same time, it does not minimize the maximum angle. Also it does not minimize the length of the edges.

Constrained…?

What does constrained mean here? That apart from the list of points, we provide a list of edges that must be in the triangulation. Meaning that no triangle edge will intersect with any of them. This is useful when triangulation is used to describe some real world situation where natural constraints occur and triangles should respect that.

Papers on the topic

There are several groups of algorithms for Delaunay triangulation that emerged over the years (this list is by no means exclusive, just what I happened to observe when researching):

• Incremental algorithms - they start with a single triangle and add points one by one.
• Divide and conquer - they divide the set of points into smaller sets and then merge them back together. Common strategy in algorithm design.
• Sweep-line algorithms - they use a line that 'scans' through the points and builds the triangulation.

Each of those methods have different properties. Divide and conquer algorithms potentially paralellize very well. Incremental algorithms usually first deal with regular triangulation and then work on enforcing constraints. This means that if constrains are dynamic, it's possible to cache the initial results. Sweep-line algorithms are the newest branch and offer interesting optimizations, such as radial sweep from a random point.

When researching NPM libraries implementing CDT, I found that most of them tend to follow one of those two papers:

• A Fast Algorithm for Generating Constrained Delaunay Triangulations – S. W. Sloan, 1992.
• Sweep-line algorithm for constrained Delaunay triangulation - V. Domiter and and B. Zalik.

Existing solutions

There are some libraries which already do this:

• Delaunator - the fastest JS library for Delaunay triangulation. If you were to study one library to learn how to achieve the best performance in computational geometry algorithms – this is the one. Although the library itself does not support constrained triangulation, there's a plugin which helps with that – Constrainautor.
• cdt-js - uses W. Sloan. Not exactly a library, more like a demo, but contains functional implementation.
• poly2tri.js - based on Zalik's paper.

If you are looking for ready working code – it might be a good idea to start with one of those. I just like to implement things myself but by no means it was dictated by having no existing solutions.

Decision

I tried implementing a divide and conquer algorithm using quad edge data structure but found it quite complicated and I was moving too slow. Then I gave a try to the sweep line algorithm but found some issues in the paper and realized that in my particular usecase I will probably prefer to cache the regular triangulation before resolving constraints. So in the end I went with Sloan's algorithm which is incremental and has a clear separation between regular and constrained triangulation phases.

Algorithm from the paper

I definitely recommend at least scrolling through the whole thing: A Fast Algorithm for Generating Constrained Delaunay Triangulations.

Here's a quote of the algorithm description from the paper (warning: this is quite long and dense, has a lot of mathematical jargon):

1. (Normalize coordinates of points.) Scale the coordinates of the points so that they all lie between 0 and 1. This scaling should be uniform so that the relative positions of the points are unchanged.
2. (Sort points into bins.) Cover the region to be triangulated by a rectangular grid so that each rectangle (or bin) contains roughly $N^{1/2}$ points. Label the bins so that consecutive bins are adjacent to one another, for example by using column-by-column or row-by-row ordering, and then allocate each point to its appropriate bin. Sort the list of points in ascending sequence of their bin numbers so that consecutive points are grouped together in the x-y plane.
3. (Establish the supertriangle.) Select three dummy points to form a supertriangle that completely encompasses all of the points to be triangulated. This supertriangle initially defines a Delaunay triangulation which is compromised of a single triangle. Its vertices are defined in terms of normalized coordinates and are usually located at a considerable distance from the window which encloses the set of points.
4. (Loop over each point.) For each point P in the list of sorted points, do steps 5-7.
5. (Insert new point in triangulation.) Find an existing triangle which encloses P. Delete this triangle and form three new triangles by connecting P to each of its vertices. The net gain in the total number of triangles after this stage is two. The searching algorithm of Lawson may be used to locate the triangle containing P efficiently. Because of the bin sorting phase, only a few triangles need to be examined if the search is initiated in the triangle which has been formed most recently.
6. (Initialize stack.) Place all triangles which are adjacent to the edges opposite P on a last-in-first-out stack. There is a maximum of three such triangles.
7. (Restore Delaunay triangulation.) While the stack of triangles of triangles is not empty, execute Lawson's swapping scheme, as defined by steps 7.1-7.3.
   1. Remove a triangle which is opposite P from the top of the stack.
   2. If P is outside (or on) the circumcircle for this triangle, return to step 7.1. Else, the triangle containing P as a vertex and the unstacked triangle form a convex quadrilateral whose diagonal is drawn in the wrong direction. Swap this diagonal so that two old triangles are replaced by two new triangles and the structure of the Delaunay triangulation is locally restored.
   3. Place any triangles which are now opposite P on the stack.

So this is how it looks like. The paper contains more figures showing different cases and a bit more context and explanations which I recommend even briefly checking. But this thing above is the core of what we need.

Constructing a constrained triangulation

Next is the constrained part.

1. (Loop over each constrained edge.) Let each constrained edge to be defined by the vertices $V_i$ and $V_j$. For each of these edges, do steps 2-4.
2. (Find intersecting edges.) If the constrained edge $V_i-V_j$ is already present in the triangulation, then go to step 1. Else, search the triangulation and store all of the edges that cross $V_i-V_j$.
3. (Remove intersecting edges.) While some edges still cross the constrained edge $V_i-V_j$, do steps 3.1 and 3.2.
   1. Remove an edge from the list of edges that intersect $V_i-V_j$. Let this edge be defined by the vertices $V_k$ and $V_l$.
   2. If the two triangles that share the edge $V_k-V_l$ do not form a quadrilateral which is strictly convex, then place $V_k-V_l$ back on the list of intersecting edges and go to step 3.1. Else, swap the diagonal of this strictly convex quadrilateral os that two new triangles are substituted for two old triangles. Let the new diagonal be defined by the vertices $V_m$ and $V_n$. If $V_m-V_n$ still intersects the constrained edge $V_i-V_j$, then place it on the list of intersecting edges. If $V_m-V_n$ does not intersect $V_i-V_lj$, then place $V_m-V_n$ on a list of newly created edges.
4. (Restore Delaunay triangulation.) Repeat steps 4.1-4.3 until no further swaps take place.
   1. Loop over each edge in the list of newly created edges.
   2. Let the newly created edge be defined by the vertices $V_k$ and $V_l$. If the edge $V_k-V_l$ is equal to the constrained edge $V_i-V_j$, then skip to step 4.1.
   3. If the two triangles that share the edge $V_k-V_l$ do not satisfy the Delaunay criterion, so that a vertex of one of the triangles is inside the circumcircle of the other triangle, then these triangles form a quadrilateral with the diagonal drawn in the wrong direction. In this case, the edge $V_k-V_l$ is swapped with the other diagonal (say) $V_m-V_n$, thus substituting two new triangles for two new triangles, and $V_k-V_l$ is replaced by $V_m-V_n$ in the list of newly created edges.
5. (Remove superflouos triangles.) Remove all triangles that contain a supertriangle vertex or lie outside the domain boundary.

Observations

In my implementation I decided to skip bin sorting and start with $O(n)$ implementations of helper functions. I just wanted to make it work first. What I found at some point is that sorting points does matter since depending on the order of input points, this algorithm might not be possible to execute – when next point turns out to be in the middle of existing triangle, there's no way to split it into three. With proper sorting it won't happen.

The point of having a super triangle is to ensure that all points are inside the triangulation since the beginning, allowing to skip more complicated setup steps. If we know that generally we are interested in rectangle-shaped areas, can't we just create two triangles that cover the bounding box of points? Not really, because then some points might be collinear with one of the edges (if a point lies exactly on the rectangle's diagonal) and then it breaks triangulation since there's no way to split triangle into three.

It's also worth noting that since vertices of the supertriangle are placed so unreasonably far from any other points, the triangulation will naturally close itself. This is because if any super triangle edge were to connect with not the most outer points closest to it, then those points would end up in the circumcircle. This also means that even though we impose the edges of the bounding box as constraints, those edges are very very likely to already exist.

Which means that we can actually move the step of removing supertriangle leftovers before the constrained triangulation step. This comes in handy if we need to first cache result of the Delaunay triangulation.

Swapping diagonal

Somehow it turns out that all we need to do to ensure that Delaunay property is satisfied for all triangles is to just swap diagonals (thanks scientists!).

The step of fixing Delaunay triangulation (number 7 in the first part) is based on maintaing a stack of triangles to process. When we split a triangle, for each neighbor of the previous triangle we check the circumcircle condition and if it's not met, we swap the diagonal of the quadrilateral (four-sided polygon) formed by the new mini-triangle and the mentioned neighbor. We then consider two triangles that were neighbors of the previous neighbor and add them to the stack.

Resolving constrained edges

The paper explains the procedure of resolving constrained edges. For that we need to find all intersecting edges and then swap them. One unobvious part of the algorithm is that if it's not possible to swap a diagonal, we just continue and swap the others and somehow later on it always becomes possible.

Code steps

It's hard to give it justice in a blog post, but after multiple attempts I wrote my own summary of the algorithm:

1. Sort points by y and if equal by x.
2. Find a bounding box of the points.
3. Create a supertriangle that contains the bounding box.
4. Loop over each point p:
   1. Find a triangle that contains p.
   2. Split the triangle into three by connecting p to each of its vertices.
   3. For each of the new triangles, add a triangle neighboring the edge opposite to the point p to a stack.
   4. While the stack is not empty:
      1. Pop the triangle from the stack.
      2. If the point p is inside the circumcircle of the triangle, swap the diagonal. Add the triangles (up to two) which are now neighboring the edge opposite to p to the stack.
5. For each edge e defined as constrained:
   1. If e is already in the triangulation, continue.
   2. Find all edges intersecting e.
   3. While there are edges intersecting e, let's call them i:
      1. Remove i from the list.
      2. Check if two triangles that are on the sides of i form a convex quadrilateral. If they do, swap the diagonal.
      3. If e and the new diagonal are crossing, put the new diagonal to the list of intersecting. Else, add it to the list of created edges.
   4. Do while there are swaps:
      1. Loop over created edges:
         1. If created edge is equal to constrained edge, continue.
         2. Find quadrilateral formed by the edge and its neighboring triangles. If the two triangles do not satisfy the Delaunay condition, swap the diagonal and replace the edge from the list of shared edges.
6. Remove all triangles that have a vertex from the supertriangle.

Designing data structures

The most inventive part of the whole process is coming up with precise data structures that will make our job easier not harder. As you might have noticed, scientific papers are very concise and try to skip all details that are not absolutely necessary for reproducing the results. This means that we need to fill in the gaps and it often means coming up ourselves with needed data structures.

Let's think what operations we need to support:

• Finding a triangle that contains a point.
• Splitting a triangle into three.
• Swapping a diagonal of two neighboring triangles.
• Finding all edges that intersect a given edge.

To achieve the first one, a spatial structure like quad tree could be utilized for fast search, but I will just loop over all triangles and check if the point is inside. Lawson's algorithm mentioned by the paper involves starting in the most recently created triangle and then marching in general direction of the point.

Splitting a triangle and swapping a diagonal are somewhat related. If we want to avoid searching all triangles every time, we can store neighboring triangle for each edge.

There are popular data structures like half edges or quad edges which store a lot of information about edges and faces (triangles) positioned next to them. I initially tried implementing them, but I found it to be a lot of overhead, which brings me to what I ended up with:

class Triangle {
  n1: Triangle | null = null;
  n2: Triangle | null = null;
  n3: Triangle | null = null;

  constructor(
    public p1: Vec2,
    public p2: Vec2,
    public p3: Vec2,
  ) {}
}

Which means there's 6 fields: n1, n2, n3 for neighboring triangles and p1, p2, p3 for vertices.

I left out some helper methods. You can see them further in the interactive example's editor.

Implementation

The implementation process involved sketching the main loop of the algorithm, isolating helper functions and then implementing them one by one and writing unit tests for them. Here's the main function:

export function cdt(
  points: Array<Vec2>,
  edges: Array<[Vec2, Vec2]>,
): Set<Triangle> {
  // Left out:
  // 1. Ensure that all edge vertices are also in the points list.
  // 2. Sort points and edges.
  // 3. Set up the supertriangle.

  // 4. Loop over each point.
  for (const p of points) {
    // Find triangle `t` containing `p`.
    const t = findTriangleWithPoint(triangles, p);

    // Remove `t` and add three new triangles connecting `p` to the 
    // vertices of `t`. 
    const [t1, t2, t3] = splitTriangle(t, p);
    triangles.add(t1);
    triangles.add(t2);
    triangles.add(t3);
    triangles.delete(t);

    const stack = [t.n1, t.n2, t.n3].filter(Boolean);

    while (stack.length > 0) {
      const tr = stack.pop()!;

      if (isInCircumcircle(tr, p)) {
        const trr = tr.findNeighborWithVertex(p);
        const [nn1, nn2, swapped1, swapped2] = swapDiagonal(tr, trr, p);

        triangles.delete(tr);
        triangles.delete(trr);
        triangles.add(swapped1);
        triangles.add(swapped2);

        // Add any triangles which are now opposite of P to the 
        // stack. 
        if (nn1) {
          stack.push(nn1);
        }
        if (nn2) {
          stack.push(nn2);
        }
      }
    }
  }

  // 5. Loop over each constrained edge.
  for (let i = 0; i < edges.length; i++) {
    // If edge is already in the triangulation, skip it.
    if (doesEdgeExist(triangles, edges[i]![0], edges[i]![1])) {
      continue;
    }

    // Find intersecting edges.
    const intersecting = findIntersectingEdges(
      triangles,
      edges[i]![0],
      edges[i]![1],
    );
    const createdEdges: Array<[Vec2, Vec2]> = [];

    while (intersecting.length > 0) {
      const _e = intersecting.shift()!;
      const e = getTrianglesForEdge(triangles, ..._e);

      // If the two triangles that share the intersecting edge do 
      // not form a strictly convex quadrilateral, place the edge 
      // back on the list. 
      const v = getQuadrilateral(e);
      if (!checkIfConvexQuadrilateral(v[0], v[1], v[2], v[3])) {
        intersecting.push(_e);
        continue;
      }

      const [, , n1, n2] = swapDiagonal(e[0]!, e[1]!, v[2]);
      triangles.delete(e[0]!);
      triangles.delete(e[1]!);
      triangles.add(n1);
      triangles.add(n2);

      const newShared = findSharedEdge(n1, n2);

      // If the new diagonal still intersects an edge, place it on 
      // the list of intersecting edges. Otherwise, add it to the 
      // list of created edges. 
      if (
        areSegmentsCrossing(
          edges[i]![0],
          edges[i]![1],
          newShared[0],
          newShared[1],
        )
      ) {
        intersecting.push(newShared);
      } else {
        createdEdges.push(newShared);
      }
    }

    // Restore delaunay triangulation.
    let swapped = false;
    do {
      swapped = false;
      for (let j = 0; j < createdEdges.length; j++) {
        const _e = createdEdges[j]!;
        const e = getTrianglesForEdge(triangles, ..._e);
        if (!e) {
          continue;
        }

        const shared = [e[0].p1, e[0].p2, e[0].p3].filter((p) =>
          e[1].hasVertex(p),
        );

        // If the edge is equal to a constrained edge, skip it.
        let skip = false;
        for (let k = 0; k < edges.length; k++) {
          if (
            (equalsEpsilonVec2(edges[k]![0], shared[0]!) &&
              equalsEpsilonVec2(edges[k]![1], shared[1]!)) ||
            (equalsEpsilonVec2(edges[k]![0], shared[1]!) &&
              equalsEpsilonVec2(edges[k]![1], shared[0]!))
          ) {
            skip = true;
          }
        }
        if (skip) {
          continue;
        }

        const [p1, , p2] = getQuadrilateral(e);

        // If two triangles that share this edge do not satisfy the 
        // delaunay condition then these triangles form a 
        // quadrilateral with the diagonal in the wrong direction. 
        if (isInCircumcircle(e[1], p1) || isInCircumcircle(e[0]!, p2)) {
          // Due to how swapDiagonal works, the point `p` must come 
          // from the second triangle. 
          const p = e[1].hasVertex(p1) ? p1 : p2;
          const [, , n1, n2] = swapDiagonal(e[0], e[1], p);
          triangles.delete(e[0]);
          triangles.delete(e[1]);
          triangles.add(n1);
          triangles.add(n2);
          createdEdges[j] = findSharedEdge(n1, n2);
          swapped = true;
        }
      }
    } while (swapped);
  }

  // 6. Remove all triangles that have a vertex from the 
  // supertriangle. 
  for (const t of triangles) {
    if (
      t.hasVertex(superTriangle.p1) ||
      t.hasVertex(superTriangle.p2) ||
      t.hasVertex(superTriangle.p3)
    ) {
      triangles.delete(t);
    }
  }

  return Array.from(triangles);
}

Helper functions

Now time to implement everything that I skipped.

Sorting points

First by y and then by x.

function comparePoints(a: Vec2, b: Vec2) {
  if (equalsEpsilon(a.x, b.x)) {
    return a.y - b.y;
  }
  return a.x - b.x;
}

Does edge exist

Inefficient function for checking if given edge exists in the triangulation.

One obvious optimization would be to find triangle containing any of the two and then just check the neighbors (as the point can't exist elsewhere).

function doesEdgeExist(triangles: Set<Triangle>, a: Vec2, b: Vec2): boolean {
  for (const t of triangles) {
    if (t.hasVertex(a) && t.hasVertex(b)) {
      return true;
    }
  }
  return false;
}

Find intersecting edges

Again, very inefficient function for finding all edges that intersect a given edge. To ensure uniqueness I used a map but that's not very necessary since usually there won't be many of them and overhead might be not worth it. Plus collisions of custom hash function are definitely possible.

There are some tricks here but be careful since edges can start and stop anywhere and it's hard to predict where other points can be.

function findIntersectingEdges(
  triangles: Set<Triangle>,
  a: Vec2,
  b: Vec2,
): Array<[Vec2, Vec2]> {
  const intersecting = new Map<number, [Vec2, Vec2]>();
  for (const t of triangles) {
    if (areSegmentsCrossing(t.p1, t.p2, a, b)) {
      intersecting.set(hashEdge(t.p1, t.p2), [t.p1, t.p2]);
    }
    if (areSegmentsCrossing(t.p2, t.p3, a, b)) {
      intersecting.set(hashEdge(t.p2, t.p3), [t.p2, t.p3]);
    }
    if (areSegmentsCrossing(t.p3, t.p1, a, b)) {
      intersecting.set(hashEdge(t.p3, t.p1), [t.p3, t.p1]);
    }
  }
  return [...intersecting.values()];
}

Find triangle with a point

Another badly-written function for finding a triangle that contains a point. It can be optimized by using a spatial structure like quad tree.

function findTriangleWithPoint(triangles: Set<Triangle>, p: Vec2) {
  for (const t of triangles) {
    if (containsPoint(t, p)) {
      return t;
    }
  }

  return null;
}

Or there's a "walking" algorithm. Starting randomly in any triangle (or in case of our specific loop, in the most recently created triangle). The idea is simple: we check if the point is inside current triangle. If not, we check which of the three edges crosses with a segment that connects centroid point of the triangle with the target point of the search and we move there. If the plane is fully triangulated, it is guaranteed to find the triangle (which has an interesting consequence: once we are done and we move on to remove triangles inside "holes" defined by constrained edges, this method will stop working).

function findTriangleWithPoint(p: Vec2, start: Triangle): Triangle | null {
  let current = start;
  let _i = 0;
  while (true) {
    _i++;
    if (_i > MAX_ITERATIONS) {
      throw new Error("Infinite loop.");
    }

    if (triangleHasPoint(current, p)) {
      return current;
    }

    const p0 = centroid(current.p1, current.p2, current.p3);
    if (areSegmentsCrossing(p0, p, current.p1, current.p2) && current.n3) {
      current = current.n3!;
      continue;
    }
    if (areSegmentsCrossing(p0, p, current.p2, current.p3) && current.n1) {
      current = current.n1!;
      continue;
    }
    if (areSegmentsCrossing(p0, p, current.p3, current.p1) && current.n2) {
      current = current.n2!;
      continue;
    }

    throw new Error("Triangle not found.");
  }
}

Get triangle for edge

One more $O(n)$ function for finding a triangle that contains a given edge.

Could be similarly optimized.

function getTrianglesForEdge(
  triangles: Set<Triangle>,
  a: Vec2,
  b: Vec2,
): [Triangle, Triangle] | null {
  for (const t of triangles) {
    if (t.hasVertex(a) && t.hasVertex(b)) {
      const thirdPoint = [t.p1, t.p2, t.p3].filter(
        (p) => !equalsEpsilonVec2(p, a) && !equalsEpsilonVec2(p, b),
      );
      invariant(thirdPoint.length === 1, "Third point must be found.");
      return [t, t.getNeighbor(thirdPoint[0]!)!];
    }
  }
  return null;
}

Find shared edge

A utility function for finding a shared edge between two triangles.

function findSharedEdge(t1: Triangle, t2: Triangle): [Vec2, Vec2] {
  const shared = [t1.p1, t1.p2, t1.p3].filter((p) => t2.hasVertex(p));
  if (shared.length !== 2) {
    throw new Error("Triangles do not share an edge.");
  }
  return shared as [Vec2, Vec2];
}

Get quadrilateral of two triangles

Given two triangles, finds a quadrilateral formed by them. Returns points in counter-clockwise order.

/**
 * `[1]` and `[3]` are shared vertices.
 * `[0]` is from `t1` and `[2]` is from `t2`.
 */
function getQuadrilateral(e: [Triangle, Triangle]): [Vec2, Vec2, Vec2, Vec2] {
  const [t1, t2] = e;
  const shared = [t1.p1, t1.p2, t1.p3].filter((p) => t2.hasVertex(p));

  // Non-shared vertices.
  const q1 = equalsEpsilonVec2(t1.nextPointCCW(shared[0]!), shared[1]!)
    ? t1.nextPointCCW(shared[1]!)
    : t1.nextPointCCW(shared[0]!);
  const q2 = equalsEpsilonVec2(t2.nextPointCCW(shared[0]!), shared[1]!)
    ? t2.nextPointCCW(shared[1]!)
    : t2.nextPointCCW(shared[0]!);

  return [q1, t1.nextPointCCW(q1), q2, t2.nextPointCCW(q2)];
}

Check if quadrilateral is convex

Or in other words, if any internal angle is greater or equal to 180 degrees, the quadrilateral is not convex. It can be done by checking if cross product of all pairs of adjacent edges is of the same sign.

function checkIfConvexQuadrilateral(
  a: Vec2,
  b: Vec2,
  c: Vec2,
  d: Vec2,
): boolean {
  const ab = new Vec2(b.x - a.x, b.y - a.y);
  const bc = new Vec2(c.x - b.x, c.y - b.y);
  const cd = new Vec2(d.x - c.x, d.y - c.y);
  const da = new Vec2(a.x - d.x, a.y - d.y);

  const crossABBC = ab.cross(bc);
  const crossBCCD = bc.cross(cd);
  const crossCDDA = cd.cross(da);
  const crossDAAB = da.cross(ab);

  return (
    (crossABBC > 0 && crossBCCD > 0 && crossCDDA > 0 && crossDAAB > 0) ||
    (crossABBC < 0 && crossBCCD < 0 && crossCDDA < 0 && crossDAAB < 0)
  );
}

Swap diagonal

Thanks to storing neighbor pointers and having CCW access to triangle points, swapping a diagonal is relatively straightforward.

/**
 * Triangle `t` and its neighbor containing point `p` form a convex quadrilateral. Swap its diagonal
 * and return up to two triangles which are now opposite of `p` and two new triangles.
 *
 * @returns two triangles which are now opposite of `p` and two new triangles created by swapping
 * the diagonal.
 */
function swapDiagonal(
  t1: Triangle,
  t2: Triangle,
  p: Vec2,
): [Triangle | null, Triangle | null, Triangle, Triangle] {
  const v0 = p;
  const v1 = t2.nextPointCCW(p);
  const v2 = t1.nextPointCCW(t2.nextPointCCW(p));
  const v3 = t2.nextPointCCW(t2.nextPointCCW(p));

  const n0 = t2.getNeighbor(v3);
  const n1 = t1.getNeighbor(v3);
  const n2 = t1.getNeighbor(v1);
  const n3 = t2.getNeighbor(v1);

  // Create new triangles.
  const nn1 = new Triangle(v0, v1, v2);
  const nn2 = new Triangle(v0, v2, v3);

  // Set neighbors.
  nn1.n1 = n1;
  n1?.replaceNeighbor(t1, nn1);
  nn1.n2 = nn2;
  nn1.n3 = n0;
  n0?.replaceNeighbor(t2, nn1);
  nn2.n1 = n2;
  n2?.replaceNeighbor(t1, nn2);
  nn2.n2 = n3;
  n3?.replaceNeighbor(t2, nn2);
  nn2.n3 = nn1;

  return [n1, n2, nn1, nn2];
}

Split triangle

TO split a triangle into three, we need to create three new triangles and update their neighbors and neighbors of the original triangle.

function splitTriangle(t: Triangle, p: Vec2): [Triangle, Triangle, Triangle] {
  // Create new triangles.
  const pab = new Triangle(p, t.p1, t.p2);
  const pbc = new Triangle(p, t.p2, t.p3);
  const pca = new Triangle(p, t.p3, t.p1);

  // Update neighbors of T.
  t.n1?.replaceNeighbor(t, pbc);
  t.n2?.replaceNeighbor(t, pca);
  t.n3?.replaceNeighbor(t, pab);

  pab.n1 = t.n3;
  pab.n2 = pbc;
  pab.n3 = pca;

  pbc.n1 = t.n1;
  pbc.n2 = pca;
  pbc.n3 = pab;

  pca.n1 = t.n2;
  pca.n2 = pab;
  pca.n3 = pbc;

  return [pab, pbc, pca];
}

Is in circumcircle

Checks if point D is inside the circumcircle of triangle ABC.

function isInCircumcircle(t: Triangle, d: Vec2): boolean {
  const ax = t.p1.x - d.x;
  const ay = t.p1.y - d.y;
  const bx = t.p2.x - d.x;
  const by = t.p2.y - d.y;
  const cx = t.p3.x - d.x;
  const cy = t.p3.y - d.y;

  return (
    (ax * ax + ay * ay) * (bx * cy - cx * by) -
      (bx * bx + by * by) * (ax * cy - cx * ay) +
      (cx * cx + cy * cy) * (ax * by - bx * ay) >=
    0
  );
}

Does triangle contain point

The helper sign determines on which side of a half-plane defined by AB the point C is (it's a cross product of vectors AC and BC). If all signs are the same, the point is inside the triangle. Then the main function checks if D is inside triangle ABC by checking if it's on the same side of all edges.

function sign(a: Vec2, b: Vec2, c: Vec2): number {
  return (a.x - c.x) * (b.y - c.y) - (b.x - c.x) * (a.y - c.y);
}

function containsPoint(t: Triangle, d: Vec2): boolean {
  const d1 = sign(d, t.p1, t.p2);
  const d2 = sign(d, t.p2, t.p3);
  const d3 = sign(d, t.p3, t.p1);
  return !((d1 < 0 || d2 < 0 || d3 < 0) && (d1 > 0 || d2 > 0 || d3 > 0));
}

Is triangle counter-clockwise

Which is basically sign but in a different order (mathematically speaking one of them is cross and the other dot product of two vectors).

function isCCW(a: Vec2, b: Vec2, c: Vec2): boolean {
  return (b.x - a.x) * (c.y - a.y) - (c.x - a.x) * (b.y - a.y) > 0;
}

Result

Here is an interactive example. Click anywhere within the boundaries to create a point. Hover any triangle to see its circumcircle. "Don't constrain edges" turns off enforcing the constraints. "Don't clear super triangle" shows the edges of the gigantic triangle that starts the process.

The algorithm uses several while loops so there are ways to make it go into an infinite loop. Mostly unusual edge cases – but if you see the box freeze that's probably it. Algorithm itself runs in a matter of milliseconds even for hundreds or thousands of points.

Improvements

Next step would be to implement all the improvements that we skipped before:

• Bin sorting.
• Better algorithm for finding intersecting edges.
• Spatial structures and optimizing all $O(n)$ search cases.

Conclusion

I hope that this gives a good overview of the process of implementing an algorithm from a scientific paper. Most computer science papers follow similar conventions so knowing how to read and implement one transfers well to others.

If anything is unclear or you have any questions, feel free to DM me on Twitter.