Moved to BowyerWatson Algorithm

DelauneyTriangulation/Models/Geometries.cs

@@ -1,4 +1,5 @@
 using System;
+using System.Collections;
 using System.Collections.Generic;
 using Avalonia.Media;
 using CommunityToolkit.Mvvm.Input;
@@ -49,17 +50,59 @@ public static class Geometries
 
         public Avalonia.Point Start => new(X1, Y1);
         public Avalonia.Point End => new(X2, Y2);
+
+        private static bool Almost(double a, double b, double eps = 1e-9) => Math.Abs(a - b) <= eps;
+
+        public bool EqualsUndirected(Edge other, double eps = 1e-9)
+        {
+            bool same = Almost(X1, other.X1, eps) && Almost(Y1, other.Y1, eps) &&
+                        Almost(X2, other.X2, eps) && Almost(Y2, other.Y2, eps);
+
+            bool swapped = Almost(X1, other.X2, eps) && Almost(Y1, other.Y2, eps) &&
+                           Almost(X2, other.X1, eps) && Almost(Y2, other.Y1, eps);
+
+            return same || swapped;
+        }
     }
 
-    public sealed class Circle(double x, double y, double r)
+    public sealed class Circle
     {
-        public double X { get; } = x;
+        public double X { get; }
-        public double Y { get; } = y;
+        public double Y { get; }
-        public double R { get; } = r;
+        public double R { get; }
 
         public double Left => X - R;
         public double Top => Y - R;
         public double Diameter => 2 * R;
+
+        public Circle(double x, double y, double r)
+        {
+            X = x;
+            Y = y;
+            R = r;
+        }
+
+        public Circle(Point p1, Point p2, Point p3)
+        {
+            double x1 = p1.X, y1 = p1.Y;
+            double x2 = p2.X, y2 = p2.Y;
+            double x3 = p3.X, y3 = p3.Y;
+
+            var d = 2.0 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));
+
+            if (Math.Abs(d) < 1e-12) throw new ArgumentException("Points are collinear; circumcircle undefined.");
+
+            var x1p = x1 * x1 + y1 * y1;
+            var x2p = x2 * x2 + y2 * y2;
+            var x3p = x3 * x3 + y3 * y3;
+
+            var ux = (x1p * (y2 - y3) + x2p * (y3 - y1) + x3p * (y1 - y2)) / d;
+            var uy = (x1p * (x3 - x2) + x2p * (x1 - x3) + x3p * (x2 - x1)) / d;
+
+            X = ux;
+            Y = uy;
+            R = Math.Sqrt((ux - x1) * (ux - x1) + (uy - y1) * (uy - y1));
+        }
     }
 
     public sealed class Triangle(Point p1, Point p2, Point p3)

DelauneyTriangulation/ViewModels/MainWindowViewModel.cs

@@ -49,11 +49,53 @@ public partial class MainWindowViewModel : ViewModelBase
     public IAsyncRelayCommand GenerateCommand { get; }
 
     // gets distance of a and b edge to calc hypotenuse aka dist of two points
-    double PointDistance(Geometries.Point p1, Geometries.Point p2) => Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));
+    double PointDistance(Geometries.Point p1, Geometries.Point p2) =>
+        Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));
+
+    static bool Almost(double a, double b, double eps = 1e-9) => Math.Abs(a - b) <= eps;
+
+    bool IsPointInCircumcircle(Geometries.Triangle t, Geometries.Point p)
+    {
+        var a = new Geometries.Point(t.X.X, t.X.Y);
+        var b = new Geometries.Point(t.Y.X, t.Y.Y);
+        var c = new Geometries.Point(t.Z.X, t.Z.Y);
+
+        var circ = new Geometries.Circle(a, b, c);
+        var d = PointDistance(new(circ.X, circ.Y), p);
+
+        return d <= circ.R + 1e-9;
+    }
+
+    static bool UsesVertex(Geometries.Triangle t, Geometries.Point v)
+    {
+        bool isV(Avalonia.Point q) => Math.Abs(q.X - v.X) <= 1e-9 && Math.Abs(q.Y - v.Y) <= 1e-9;
+        return isV(t.X) || isV(t.Y) || isV(t.Z);
+    }
+
+    (Geometries.Point A, Geometries.Point B, Geometries.Point C) BuildSuperTriangle()
+    {
+        var minX = Points.Min(p => p.X);
+        var minY = Points.Min(p => p.Y);
+        var maxX = Points.Max(p => p.X);
+        var maxY = Points.Max(p => p.Y);
+
+        var dx = maxX - minX;
+        var dy = maxY - minY;
+        var delta = Math.Max(dx, dy);
+        var cx = (minX + maxX) * 0.5;
+        var cy = (minY + maxY) * 0.5;
+
+        var A = new Geometries.Point(cx - 2 * delta, cy - delta * 3);
+        var B = new Geometries.Point(cx, cy + 3 * delta);
+        var C = new Geometries.Point(cx + 2 * delta, cy - delta * 3);
+
+        return (A, B, C);
+    }
 
     async Task GenerateRandomPoints(int count, CancellationToken ct)
     {
-        bool CheckAdd(Geometries.Point p1, IList<Geometries.Point> p) => !p.Any(x => PointDistance(p1, x) <= MinPointDistance);
+        bool CheckAdd(Geometries.Point p1, IList<Geometries.Point> p) =>
+            !p.Any(x => PointDistance(p1, x) <= MinPointDistance);
 
         Points.Clear();
         var rand = new Random();
@@ -85,7 +127,9 @@ public partial class MainWindowViewModel : ViewModelBase
             {
                 i--;
                 retry++;
-            };
+            }
+
+            ;
         }
 
         l.Sort(Geometries.PointComparer.Instance);
@@ -97,7 +141,8 @@ public partial class MainWindowViewModel : ViewModelBase
         }
     }
 
-    void ConnectPoint(Geometries.Point p1, Geometries.Point p2) => Edges.Add(new Geometries.Edge(p1.X, p1.Y, p2.X, p2.Y));
+    void ConnectPoint(Geometries.Point p1, Geometries.Point p2) =>
+        Edges.Add(new Geometries.Edge(p1.X, p1.Y, p2.X, p2.Y));
 
     async Task ConnectPoints(CancellationToken ct)
     {
@@ -122,6 +167,7 @@ public partial class MainWindowViewModel : ViewModelBase
     }
 
     void AddTriangle(Geometries.Point p1, Geometries.Point p2, Geometries.Point p3) => Triangles.Add(new(p1, p2, p3));
+
     async Task Triangulate(CancellationToken ct)
     {
         Triangles.Clear();
@@ -133,7 +179,8 @@ public partial class MainWindowViewModel : ViewModelBase
         var p2 = toTriangulate.Where(x => !ReferenceEquals(x, p1)).MinBy(x => PointDistance(p1, x));
         if (p2 is null) return;
 
-        static double TwiceArea(Geometries.Point a, Geometries.Point b, Geometries.Point c) => Math.Abs((b.X - a.X) * (c.Y - a.Y) - (b.Y - a.Y) * (c.X - a.X));
+        static double TwiceArea(Geometries.Point a, Geometries.Point b, Geometries.Point c) =>
+            Math.Abs((b.X - a.X) * (c.Y - a.Y) - (b.Y - a.Y) * (c.X - a.X));
 
         while (toTriangulate.Count >= 3)
         {
@@ -144,7 +191,8 @@ public partial class MainWindowViewModel : ViewModelBase
             Geometries.Point? p3 = null;
             foreach (var c in candidates.OrderBy(x => PointDistance(p1, x) + PointDistance(p2, x)))
             {
-                if (TwiceArea(p1, p2, c) > 1e-6) p3 = c; break;
+                if (TwiceArea(p1, p2, c) > 1e-6) p3 = c;
+                break;
             }
 
             if (p3 is null) break;
@@ -159,6 +207,125 @@ public partial class MainWindowViewModel : ViewModelBase
         }
     }
 
+    List<Geometries.Edge> ExtractTriangleEdges(Geometries.Triangle t)
+    {
+        var e1 = new Geometries.Edge(t.X.X, t.X.Y, t.Y.X, t.Y.Y);
+        var e2 = new Geometries.Edge(t.Y.X, t.Y.Y, t.Z.X, t.Z.Y);
+        var e3 = new Geometries.Edge(t.Z.X, t.Z.Y, t.X.X, t.X.Y);
+        return [e1, e2, e3];
+    }
+
+    static string CanonKey(Geometries.Edge e)
+    {
+        var abFirst = (e.X1 < e.X2) || (Almost(e.X1, e.X2) && e.Y1 <= e.Y2);
+        var ax = abFirst ? e.X1 : e.X2;
+        var ay = abFirst ? e.Y1 : e.Y2;
+        var bx = abFirst ? e.X2 : e.X1;
+        var by = abFirst ? e.Y2 : e.Y1;
+        static double Q(double v) => Math.Round(v, 9, MidpointRounding.AwayFromZero);
+
+        return $"{Q(ax)},{Q(ay)}|{Q(bx)},{Q(by)}";
+    }
+
+    void DrawCircles(IEnumerable<Geometries.Triangle> triangles, Geometries.Point superA, Geometries.Point superB,
+        Geometries.Point superC)
+    {
+        Circles.Clear();
+        foreach (var t in triangles)
+        {
+            if (UsesVertex(t, superA) || UsesVertex(t, superB) || UsesVertex(t, superC)) continue;
+            try
+            {
+                var a = new Geometries.Point(t.X.X, t.X.Y);
+                var b = new Geometries.Point(t.Y.X, t.Y.Y);
+                var c = new Geometries.Point(t.Z.X, t.Z.Y);
+
+                var cc = new Geometries.Circle(a, b, c);
+                if (double.IsFinite(cc.X) && double.IsFinite(cc.Y) && double.IsFinite(cc.R) && cc.R > 0)
+                    Circles.Add(cc);
+            }
+            catch (ArgumentException)
+            {
+            }
+        }
+    }
+
+    async Task BowyerWatson(CancellationToken ct)
+    {
+        Triangles.Clear();
+
+        if (Points.Count < 3) return;
+
+        var (sa, sb, sc) = BuildSuperTriangle();
+        Triangles.Add(new Geometries.Triangle(sa, sb, sc));
+        await DelayAsync(ct);
+
+        foreach (var p in Points)
+        {
+            ct.ThrowIfCancellationRequested();
+
+            var bad = new List<Geometries.Triangle>();
+            foreach (var t in Triangles)
+            {
+                try
+                {
+                    if (IsPointInCircumcircle(t, p)) bad.Add(t);
+                }
+                catch (ArgumentException)
+                {
+                    // kollinear -> explicit ignore
+                }
+            }
+
+            DrawCircles(bad, sa, sb, sc);
+            await DelayAsync(ct);
+
+            var edgeCount = new Dictionary<string, Geometries.Edge>();
+            var counter = new Dictionary<string, int>();
+
+            foreach (var bt in bad)
+            {
+                foreach (var e in ExtractTriangleEdges(bt))
+                {
+                    string key = CanonKey(e);
+                    if (!edgeCount.ContainsKey(key)) edgeCount[key] = e;
+                    counter.TryGetValue(key, out int c);
+                    counter[key] = c + 1;
+                }
+            }
+
+            var polygon = new List<Geometries.Edge>();
+            foreach (var kv in counter)
+                if (kv.Value == 1)
+                    polygon.Add(edgeCount[kv.Key]);
+
+            // remove bad triangles form triangulation
+            foreach (var bt in bad) Triangles.Remove(bt);
+
+            // re-triangulate the polygonal hole
+            foreach (var e in polygon)
+            {
+                var a = new Geometries.Point(e.X1, e.Y1);
+                var b = new Geometries.Point(e.X2, e.Y2);
+
+                Triangles.Add(new Geometries.Triangle(a, b, p));
+                await DelayAsync(ct);
+            }
+        }
+
+        // remove triangles containing super triangle vertex
+        var toRemove = new List<Geometries.Triangle>();
+        foreach (var t in Triangles)
+        {
+            if (UsesVertex(t, sa) || UsesVertex(t, sb) || UsesVertex(t, sc)) toRemove.Add(t);
+        }
+
+        foreach (var t in toRemove) Triangles.Remove(t);
+
+        DrawCircles(Triangles, sa, sb, sc);
+        await DelayAsync(ct);
+    }
+
     async Task DelaunayTriangulateAsync(CancellationToken ct)
     {
         Triangles.Clear();
@@ -198,6 +365,7 @@ public partial class MainWindowViewModel : ViewModelBase
             p1h = p2h;
             await DelayAsync(ct);
         }
+
         Edges.Add(new Geometries.Edge(first.X, first.Y, last.X, last.Y));
 
         // Voronoi Diagram Extraction
@@ -209,7 +377,8 @@ public partial class MainWindowViewModel : ViewModelBase
             {
                 var s = cell.Points[i];
                 VoronoiDiagram.Add(new Geometries.Edge(f.X, f.Y, s.X, s.Y));
-                Circles.Add(new Geometries.Circle(f.X, f.Y, Points.Min(x => PointDistance(new Geometries.Point(f.X, f.Y), x))));
+                Circles.Add(new Geometries.Circle(f.X, f.Y,
+                    Points.Min(x => PointDistance(new Geometries.Point(f.X, f.Y), x))));
 
                 f = s;
                 await DelayAsync(ct);
@@ -217,15 +386,20 @@ public partial class MainWindowViewModel : ViewModelBase
         }
     }
 
-    private Task DelayAsync(CancellationToken ct) => GenerationDelay > 0 ? Task.Delay(GenerationDelay ?? 0, ct) : Task.CompletedTask;
+    private Task DelayAsync(CancellationToken ct) =>
+        GenerationDelay > 0 ? Task.Delay(GenerationDelay ?? 0, ct) : Task.CompletedTask;
 
     async Task RunGenerateAsync()
     {
         _genCts?.CancelAsync();
         _genCts = new CancellationTokenSource();
-        try { await Generate(_genCts.Token); }
+        try
-        catch (OperationCanceledException) { }
+        {
-
+            await Generate(_genCts.Token);
+        }
+        catch (OperationCanceledException)
+        {
+        }
     }
 
     private async Task Generate(CancellationToken ct)
@@ -239,7 +413,8 @@ public partial class MainWindowViewModel : ViewModelBase
         await GenerateRandomPoints(PointCount ?? 20, ct);
         // await ConnectPoints(ct);
         // await Triangulate(ct);
-        await DelaunayTriangulateAsync(ct);
+        // await DelaunayTriangulateAsync(ct);
+        await BowyerWatson(ct);
     }
 
     private bool CanGenerate() => true;