Delaunay triangulation based on Paul Bourke's algorithm (http://paulbourke.net/papers/triangulate/)

Delaunay.cs

using OpenTK;
using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

static class Delaunay
{
    [DebuggerDisplay("({V0}, {V1}, {V2})")]
    struct Triangle
    {
        public int V0;
        public int V1;
        public int V2;

        public Triangle(int v0, int v1, int v2)
        {
            V0 = v0;
            V1 = v1;
            V2 = v2;
        }
    }

    [DebuggerDisplay("({V0}, {V1})")]
    struct Edge : IEquatable<Edge>
    {
        public int V0;
        public int V1;

        public Edge(int v0, int v1)
        {
            V0 = v0;
            V1 = v1;
        }

        public bool Equals(Edge other)
        {
            return V0 == other.V0 && V1 == other.V1
                || V0 == other.V1 && V1 == other.V0;
        }

        public override bool Equals(object obj)
        {
            if (obj is Edge)
            {
                return Equals((Edge)obj);
            }
            return false;
        }

        public override int GetHashCode()
        {
            var hash = 37;
            hash += 61 * V0.GetHashCode();
            hash += 61 * V1.GetHashCode();
            return hash;
        }
    }

    static void Circumcircle(ref Vector2 v0, ref Vector2 v1, ref Vector2 v2, out Vector2 circumcenter, out float circumradius)
    {
        var y0y1abs = Math.Abs(v0.Y - v1.Y);
        var y1y2abs = Math.Abs(v1.Y - v2.Y);
        if (y0y1abs < float.Epsilon || y1y2abs < float.Epsilon)
        {
            circumcenter = float.NaN * Vector2.One;
            circumradius = float.NaN;
        }

        if (y0y1abs < float.Epsilon)
        {
            var m2 = -(v2.X - v1.X) / (v2.Y - v1.Y);
            var mv1v2 = (v1 + v2) / 2.0f;
            circumcenter.X = (v1.X + v0.X) / 2.0f;
            circumcenter.Y = m2 * (circumcenter.X - mv1v2.X) + mv1v2.Y;
        }
        else if (y1y2abs < float.Epsilon)
        {
            var m1 = -(v1.X - v0.X) / (v1.Y - v0.Y);
            var mv0v1 = (v0 + v1) / 2.0f;
            circumcenter.X = (v2.X + v1.X) / 2.0f;
            circumcenter.Y = m1 * (circumcenter.X - mv0v1.X) + mv0v1.Y;
        }
        else
        {
            var m1 = -(v1.X - v0.X) / (v1.Y - v0.Y);
            var m2 = -(v2.X - v1.X) / (v2.Y - v1.Y);
            var mv0v1 = (v0 + v1) / 2.0f;
            var mv1v2 = (v1 + v2) / 2.0f;
            circumcenter.X = (m1 * mv0v1.X - m2 * mv1v2.X + mv1v2.Y - mv0v1.Y) / (m1 - m2);
            if (y0y1abs > y1y2abs)
            {
                circumcenter.Y = m1 * (circumcenter.X - mv0v1.X) + mv0v1.Y;
            }
            else
            {
                circumcenter.Y = m2 * (circumcenter.X - mv1v2.X) + mv1v2.Y;
            }
        }

        var d = v1 - circumcenter;
        circumradius = d.X * d.X + d.Y * d.Y;
    }

    static Triangle CreateSupertriangle(Vector2[] vertices, out Vector2[] supervertices)
    {
        var min = vertices[0];
        var max = vertices[0];
        for (int i = 1; i < vertices.Length; i++)
        {
            Vector2.ComponentMin(ref min, ref vertices[i], out min);
            Vector2.ComponentMax(ref max, ref vertices[i], out max);
        }

        var d = max - min;
        var mid = (max + min) / 2;
        var dmax = d.X > d.Y ? d.X : d.Y;
        supervertices = new Vector2[] {
            new Vector2(mid.X - 2 * dmax, mid.Y - dmax),
            new Vector2(mid.X, mid.Y + 2 * dmax),
            new Vector2(mid.X + 2 * dmax, mid.Y - dmax) };
        return new Triangle(-1, -2, -3);
    }

    public static Vector2[] Triangulate(Vector2[] vertices)
    {
        if (vertices.Length < 3) return new Vector2[0];
        Vector2[] supervertices;
        var triangles = new List<Triangle>();
        var supertriangle = CreateSupertriangle(vertices, out supervertices);
        triangles.Add(supertriangle);

        for (int i = 0; i < vertices.Length; i++)
        {
            var edges = new HashSet<Edge>();
            triangles.RemoveAll(triangle =>
            {
                float circumradius;
                Vector2 circumcenter;
                var v0 = triangle.V0 < 0 ? supervertices[-triangle.V0 - 1] : vertices[triangle.V0];
                var v1 = triangle.V1 < 0 ? supervertices[-triangle.V1 - 1] : vertices[triangle.V1];
                var v2 = triangle.V2 < 0 ? supervertices[-triangle.V2 - 1] : vertices[triangle.V2];
                Circumcircle(ref v0, ref v1, ref v2, out circumcenter, out circumradius);
                var dc = vertices[i] - circumcenter;
                var distance = dc.X * dc.X + dc.Y * dc.Y;
                if ((distance - circumradius) <= float.Epsilon)
                {
                    var edge0 = new Edge(triangle.V0, triangle.V1);
                    var edge1 = new Edge(triangle.V1, triangle.V2);
                    var edge2 = new Edge(triangle.V2, triangle.V0);
                    if (!edges.Add(edge0)) edges.Remove(edge0);
                    if (!edges.Add(edge1)) edges.Remove(edge1);
                    if (!edges.Add(edge2)) edges.Remove(edge2);
                    return true;
                }
                return false;
            });

            foreach (var edge in edges)
            {
                triangles.Add(new Triangle(edge.V0, edge.V1, i));
            }
        }

        var voffset = 0;
        triangles.RemoveAll(triangle => triangle.V0 < 0 || triangle.V1 < 0 || triangle.V2 < 0);
        var result = new Vector2[triangles.Count * 3];
        foreach (var triangle in triangles)
        {
            result[voffset++] = vertices[triangle.V0];
            result[voffset++] = vertices[triangle.V1];
            result[voffset++] = vertices[triangle.V2];
        }
        return result;
    }
}