delaunator.cs

using System;
using System.Collections.Generic;
using System.Linq;

namespace Delaunator
{
    internal static class Util
    {
        public static readonly double EPSILON = double.Epsilon;
        public static readonly int EMPTY = int.MaxValue;

        public static int NextHalfedge(int i)
        {
            return i % 3 == 2 ? i - 2 : i + 1;
        }

        public static int PrevHalfedge(int i)
        {
            return i % 3 == 0 ? i + 2 : i - 1;
        }
    }

    public struct Point
    {
        public double x;
        public double y;


        public override string ToString()
        {
            return "{x: " + x + ", y: " + y + "}";
        }

        public double Dist2(Point p)
        {
            double dx = x - p.x;
            double dy = y - p.y;
            return (dx * dx) + (dy * dy);
        }

        public bool Orient(Point q, Point r)
        {
            return ((q.y - y) * (r.x - q.x)) - ((q.x - x) * (r.y - q.y)) < 0.0;
        }

        public Tuple<double, double> Circumdelta(Point b, Point c)
        {
            double dx = b.x - x;
            double dy = b.y - y;
            double ex = c.x - x;
            double ey = c.y - y;

            double bl = (dx * dx) + (dy * dy);
            double cl = (ex * ex) + (ey * ey);
            double t_d = (dx * ey - dy * ex);
            double d = 0.5 / (dx * ey - dy * ex);

            double rx = (ey * bl - dy * cl) * d;
            double ry = (dx * cl - ex * bl) * d;
            return new Tuple<double, double>(rx, ry);
        }

        public double Circumradius2(Point b, Point c)
        {
            (double rx, double ry) = Circumdelta(b, c);
            return rx * rx + ry * ry;
        }

        public Point Circumcenter(Point b, Point c)
        {
            (double rx, double ry) = Circumdelta(b, c);
            return new Point
            {
                x = x + rx,
                y = y + ry
            };
        }

        public bool InCircle(Point b, Point c, Point p)
        {
            double dx = x - p.x;
            double dy = y - p.y;
            double ex = b.x - p.x;
            double ey = b.y - p.y;
            double fx = c.x - p.x;
            double fy = c.y - p.y;

            double ap = dx * dx + dy * dy;
            double bp = ex * ex + ey * ey;
            double cp = fx * fx + fy * fy;

            return dx * (ey * cp - bp * fy) - dy * (ex * cp - bp * fx) + ap * (ex * fy - ey * fx) < 0.0;
        }

        public bool NearlyEquals(Point p)
        {
            return Math.Abs(x - p.x) <= Util.EPSILON && Math.Abs(y - p.y) <= Util.EPSILON;
        }
    }

    public class Trianglulation
    {
        public List<int> triangles;
        public List<int> halfedges;
        public List<int> hull;

        public Trianglulation(int n)
        {
            int maxTris = ((2 * n) - 5) * 3;

            triangles = new List<int>(maxTris);
            halfedges = new List<int>(maxTris);
            hull = new List<int>();
        }

        public int AddTriangle(int i0, int i1, int i2, int a, int b, int c)
        {
            int t = triangles.Count;

            triangles.Add(i0);
            triangles.Add(i1);
            triangles.Add(i2);

            halfedges.Add(a);
            halfedges.Add(b);
            halfedges.Add(c);

            if (a != Util.EMPTY)
            {
                halfedges[a] = t;
            }

            if (b != Util.EMPTY)
            {
                halfedges[b] = t + 1;
            }

            if (c != Util.EMPTY)
            {
                halfedges[c] = t + 2;
            }

            return t;
        }

        public int Legalize(int a, IList<Point> points, Hull hull)
        {
            int b = halfedges[a];
            int ar = Util.PrevHalfedge(a);

            if (b == Util.EMPTY)
                return ar;

            int al = Util.NextHalfedge(a);
            int bl = Util.PrevHalfedge(b);

            int p0 = triangles[ar];
            int pr = triangles[a];
            int pl = triangles[al];
            int p1 = triangles[bl];

            bool illegal = points[p0].InCircle(points[pr], points[pl], points[p1]);
            if(illegal)
            {
                triangles[a] = p1;
                triangles[b] = p0;

                int hbl = halfedges[bl];
                int har = halfedges[ar];

                if(hbl == Util.EMPTY)
                {
                    int e = hull.start;
                    while(true)
                    {
                        if (hull.tri[e] == bl)
                        {
                            hull.tri[e] = a;
                            break;
                        }

                        e = hull.next[e];
                        if (e == hull.start)
                            break;
                    }
                }

                halfedges[a] = hbl;
                halfedges[b] = har;
                halfedges[ar] = bl;

                if (hbl != Util.EMPTY)
                    halfedges[hbl] = a;
                if (har != Util.EMPTY)
                    halfedges[har] = b;
                if (bl != Util.EMPTY)
                    halfedges[bl] = ar;

                int br = Util.NextHalfedge(b);

                Legalize(a, points, hull);
                return Legalize(br, points, hull);
            }

            return ar;
        }
    }

    public class Hull
    {
        public List<int> prev;
        public List<int> next;
        public List<int> tri;
        public List<int> hash;
        public int start;
        public Point center;

        public Hull(int n, Point center, int i0, int i1, int i2, IList<Point> points)
        {
            prev = Enumerable.Repeat(0, n).ToList();
            next = Enumerable.Repeat(0, n).ToList();
            tri = Enumerable.Repeat(0, n).ToList();
            hash = Enumerable.Repeat(Util.EMPTY, (int)Math.Sqrt(n)).ToList();
            start = i0;
            this.center = center;

            next[i0] = i1;
            prev[i2] = i1;
            next[i1] = i2;
            prev[i0] = i2;
            next[i2] = i0;
            prev[i1] = i0;

            tri[i0] = 0;
            tri[i1] = 1;
            tri[i2] = 2;

            HashEdge(points[i0], i0);
            HashEdge(points[i1], i1);
            HashEdge(points[i2], i2);
        }

        public int HashKey(Point p)
        {
            double dx = p.x - center.x;
            double dy = p.y - center.y;

            double _p = dx / (Math.Abs(dx) + Math.Abs(dy));
            double a = ((dy > 0.0) ? (3.0 - _p) : (1.0 + _p))  / 4.0; // [0..1]

            int len = hash.Count;
            return ((int)Math.Floor(len * a)) % len;
        }

        public void HashEdge(Point p, int i)
        {
            int key = HashKey(p);
            hash[key] = i;
        }

        public Tuple<int, bool> FindVisibleEdge(Point p, IList<Point> points)
        {
            int _start = 0;
            int key = HashKey(p);
            int len = hash.Count;
            foreach (int idx in Enumerable.Range(0, len))
            {
                _start = hash[(key + idx) % len];
                if(_start != Util.EMPTY && next[_start] != Util.EMPTY)
                {
                    break;
                }
            }

            _start = prev[_start];
            int e = _start;

            while(!p.Orient(points[e], points[next[e]]))
            {
                e = next[e];
                if (e == _start)
                {
                    return new Tuple<int, bool>(Util.EMPTY, false);
                }
            }
            return new Tuple<int, bool>(e, e == _start);
        }
    }

    public static class Delaunator
    {
        public static Point CalcBboxCenter(IList<Point> points)
        {
            double min_x = double.PositiveInfinity;
            double min_y = double.PositiveInfinity;
            double max_x = double.NegativeInfinity;
            double max_y = double.NegativeInfinity;
            foreach(Point p in points) {
                min_x = Math.Min(min_x, p.x);
                min_y = Math.Min(min_y, p.y);
                max_x = Math.Max(max_x, p.x);
                max_y = Math.Max(max_y, p.y);
            }
            return new Point
            {
                x = (min_x + max_x) / 2.0,
                y = (min_y + max_y) / 2.0
            };
        }

        public static Tuple<bool, int> FindClosestPoint(IList<Point> points, Point p0)
        {
            double min_dist = double.PositiveInfinity;
            int k = 0;

            foreach (var it in points.Select((p, i) => new { p, i }))
            {
                double d = p0.Dist2(it.p);
                if (d > 0.0 && d < min_dist)
                {
                    k = it.i;
                    min_dist = d;
                }
            }

            if (double.IsPositiveInfinity(min_dist))
            {
                return new Tuple<bool, int>(false, 0);
            }
            return new Tuple<bool, int>(true, k);
        }

        public static Tuple<bool, int, int, int> FindSeedTriangle(IList<Point> points)
        {
            Point bbox_center = CalcBboxCenter(points);
            Tuple<bool, int> i0 = FindClosestPoint(points, bbox_center);
            // if i0 is false, throw...
            Point p0 = points[i0.Item2];

            Tuple<bool, int> i1 = FindClosestPoint(points, p0);
            // ^ same
            Point p1 = points[i1.Item2];

            double min_radius = double.PositiveInfinity;
            int i2 = 0;
            foreach (var it in points.Select((p, i) => new { p, i }))
            {
                if (it.i == i0.Item2 || it.i == i1.Item2)
                    continue;
                double r = p0.Circumradius2(p1, it.p);
                if (r < min_radius)
                {
                    i2 = it.i;
                    min_radius = r;
                }
            }

            if (double.IsPositiveInfinity(min_radius))
                return new Tuple<bool, int, int, int>(false, 0, 0, 0);
            if (p0.Orient(p1, points[i2]))
                return new Tuple<bool, int, int, int>(true, i0.Item2, i2, i1.Item2);
            return new Tuple<bool, int, int, int>(true, i0.Item2, i1.Item2, i2);
        }

        private static int ComparePointDists(Tuple<int, double> a, Tuple<int, double> b)
        {
            // 1 = a is greater
            // 0 = a and b are equal
            // -1 = b is greater
            if (a == null)
            {
                if (b == null)
                {
                    return 0;
                }

                return -1;
            }

            if (b == null)
                return 1;

            return a.Item2.CompareTo(b.Item2);
        }

        public static Trianglulation Triangulate(IList<Point> points)
        {
            int pcount = points.Count;

            Tuple<bool, int, int, int> i012 = FindSeedTriangle(points);
            int i0 = i012.Item2;
            int i1 = i012.Item3;
            int i2 = i012.Item4;
            Point center = points[i0].Circumcenter(points[i1], points[i2]);

            Trianglulation triangulation = new Trianglulation(pcount);
            triangulation.AddTriangle(i0, i1, i2, Util.EMPTY, Util.EMPTY, Util.EMPTY);

            // sort the points by distance from seed triangle
            List<Tuple<int, double>> dists = points.Select((point, i) => new Tuple<int, double>( i, center.Dist2(point) )).ToList();

            dists.Sort(ComparePointDists);

            Hull hull = new Hull(pcount, center, i0, i1, i2, points);

            foreach (var it in dists.Select((t, idx) => new { k = idx, i = t.Item1 }))
            {
                Point p = points[it.i];

                if (it.k > 0 && p.NearlyEquals(points[dists[it.k - 1].Item1]))
                {
                    continue;
                }

                if (it.i == i0 || it.i == i1 || it.i == i2)
                {
                    continue;
                }

                (int e, bool walk_back) = hull.FindVisibleEdge(p, points);
                if(e == Util.EMPTY)
                {
                    continue;
                }

                int t = triangulation.AddTriangle(e, it.i, hull.next[e], Util.EMPTY, Util.EMPTY, hull.tri[e]);

                hull.tri[it.i] = triangulation.Legalize(t + 2, points, hull);
                hull.tri[e] = t;

                int m = hull.next[e];
                while(true)
                {
                    int q = hull.next[m];
                    if(!p.Orient(points[m], points[q]))
                    {
                        break;
                    }
                    int _t = triangulation.AddTriangle(m, it.i, q, hull.tri[it.i], Util.EMPTY, hull.tri[m]);
                    hull.tri[it.i] = triangulation.Legalize(_t + 2, points, hull);
                    hull.next[m] = Util.EMPTY; // mark as removed
                    m = q;
                }

                if(walk_back)
                {
                    while(true)
                    {
                        int q = hull.prev[e];
                        if(!p.Orient(points[q], points[e])) {
                            break;
                        }
                        int _t = triangulation.AddTriangle(q, it.i, e, Util.EMPTY, hull.tri[e], hull.tri[q]);
                        triangulation.Legalize(_t + 2, points, hull);
                        hull.tri[q] = _t;
                        hull.next[e] = Util.EMPTY; // mark as removed
                        e = q;
                    }
                }

                hull.prev[it.i] = e;
                hull.next[it.i] = m;
                hull.prev[m] = it.i;
                hull.next[e] = it.i;
                hull.start = e;

                hull.HashEdge(p, it.i);
                hull.HashEdge(points[e], e);
            }

            int ex = hull.start;
            while(true)
            {
                triangulation.hull.Add(ex);
                ex = hull.next[ex];
                if(ex == hull.start)
                {
                    break;
                }
            }

            triangulation.triangles.TrimExcess();
            triangulation.halfedges.TrimExcess();

            return triangulation;
        }
    }
}