diff --git a/DelauneyTriangulation/Models/Geometries.cs b/DelauneyTriangulation/Models/Geometries.cs index d66ea8a..a7b4b81 100644 --- a/DelauneyTriangulation/Models/Geometries.cs +++ b/DelauneyTriangulation/Models/Geometries.cs @@ -1,4 +1,5 @@ using System; +using System.Collections; using System.Collections.Generic; using Avalonia.Media; using CommunityToolkit.Mvvm.Input; @@ -12,14 +13,14 @@ public static class Geometries public double X { get; } = x; public double Y { get; } = y; public double R { get; } = r; - + public Color Color { get; set; } = color ?? new Color(255, 255, 255, 255); - + public double Left => X - R; public double Top => Y - R; public double Size => 2 * R; } - + public sealed class PointComparer : IComparer { public static readonly PointComparer Instance = new(); @@ -46,20 +47,62 @@ public static class Geometries public double X2 { get; } = x2; public double Y1 { get; } = y1; public double Y2 { get; } = y2; - + 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 Y { get; } = y; - public double R { get; } = r; - + public double X { get; } + public double Y { get; } + 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) @@ -67,7 +110,7 @@ public static class Geometries public Avalonia.Point X { get; } = new(p1.X, p1.Y); public Avalonia.Point Y { get; } = new(p2.X, p2.Y); public Avalonia.Point Z { get; } = new(p3.X, p3.Y); - - public IReadOnlyList Vertices => [ X, Y, Z ]; + + public IReadOnlyList Vertices => [X, Y, Z]; } } \ No newline at end of file diff --git a/DelauneyTriangulation/ViewModels/MainWindowViewModel.cs b/DelauneyTriangulation/ViewModels/MainWindowViewModel.cs index b8501cf..ff73718 100644 --- a/DelauneyTriangulation/ViewModels/MainWindowViewModel.cs +++ b/DelauneyTriangulation/ViewModels/MainWindowViewModel.cs @@ -24,16 +24,16 @@ public partial class MainWindowViewModel : ViewModelBase [ObservableProperty] private double _panX; [ObservableProperty] private double _panY; [ObservableProperty] private double _zoom = 1.0; - + [ObservableProperty] private bool _pointsVisible = true; [ObservableProperty] private bool _edgesVisible = true; [ObservableProperty] private bool _circlesVisible = true; [ObservableProperty] private bool _voronoiVisible = true; [ObservableProperty] private bool _trianglesVisible = true; - + [ObservableProperty] private int? _generationDelay = 20; private CancellationTokenSource? _genCts; - + public double Width { get; set; } = 1280; public double Height { get; set; } = 720; @@ -47,14 +47,56 @@ public partial class MainWindowViewModel : ViewModelBase public ObservableCollection Triangles { get; } = new(); 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 p) => !p.Any(x => PointDistance(p1, x) <= MinPointDistance); - + bool CheckAdd(Geometries.Point p1, IList p) => + !p.Any(x => PointDistance(p1, x) <= MinPointDistance); + Points.Clear(); var rand = new Random(); var rgb = (byte)(255 - PointCount); @@ -68,7 +110,7 @@ public partial class MainWindowViewModel : ViewModelBase { ct.ThrowIfCancellationRequested(); if (retry == 10) break; - + var x = rand.NextDouble() * (Width - LeftOffset) + LeftOffset; var y = rand.NextDouble() * (Height - TopOffset * 2) + TopOffset; var p = new Geometries.Point(x, y); @@ -76,8 +118,8 @@ public partial class MainWindowViewModel : ViewModelBase { p.Color = color; rgb++; - color = new Color(rgb, rgb, rgb, a); - + color = new Color(rgb, rgb, rgb, a); + l.Add(p); retry = 0; } @@ -85,9 +127,11 @@ public partial class MainWindowViewModel : ViewModelBase { i--; retry++; - }; + } + + ; } - + l.Sort(Geometries.PointComparer.Instance); foreach (var p in l) { @@ -96,25 +140,26 @@ public partial class MainWindowViewModel : ViewModelBase await DelayAsync(ct); } } - - 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) { - if(Points.Count == 0) return; + if (Points.Count == 0) return; Edges.Clear(); var toConnect = Points.ToList(); var p1 = toConnect[0]; - + while (toConnect.Count > 0) { ct.ThrowIfCancellationRequested(); var p2 = toConnect.MinBy(x => PointDistance(p1, x)); if (p2 is null) continue; - + ConnectPoint(p1, p2); - + toConnect.Remove(p1); p1 = p2; await DelayAsync(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,20 +179,22 @@ 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) { ct.ThrowIfCancellationRequested(); - + var candidates = toTriangulate.Where(x => !ReferenceEquals(x, p1) && !ReferenceEquals(x, p2)); 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; AddTriangle(p1, p2, p3); @@ -154,17 +202,136 @@ public partial class MainWindowViewModel : ViewModelBase toTriangulate.Remove(p1); p1 = p2; p2 = p3; - + await DelayAsync(ct); } } + List 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 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(); + 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(); + var counter = new Dictionary(); + + 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(); + 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(); + 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(); Edges.Clear(); - - if(Points.Count < 3) return; + + if (Points.Count < 3) return; var pts = Points.Select(p => (DelaunatorSharp.IPoint)new DelaunatorSharp.Point(p.X, p.Y)).ToArray(); @@ -185,7 +352,7 @@ public partial class MainWindowViewModel : ViewModelBase Triangles.Add(new Geometries.Triangle(p0, p1, p2)); await DelayAsync(ct); } - + // Konvex Hull Extraction var hullPoints = d.GetHullPoints(); var first = hullPoints.First(); @@ -198,8 +365,9 @@ 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 foreach (var cell in d.GetVoronoiCells()) { @@ -209,23 +377,29 @@ 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); } } } - 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); } - catch (OperationCanceledException) { } - + try + { + await Generate(_genCts.Token); + } + catch (OperationCanceledException) + { + } } private async Task Generate(CancellationToken ct) @@ -235,13 +409,14 @@ public partial class MainWindowViewModel : ViewModelBase Circles.Clear(); Triangles.Clear(); VoronoiDiagram.Clear(); - + await GenerateRandomPoints(PointCount ?? 20, ct); // await ConnectPoints(ct); // await Triangulate(ct); - await DelaunayTriangulateAsync(ct); + // await DelaunayTriangulateAsync(ct); + await BowyerWatson(ct); } - + private bool CanGenerate() => true; public MainWindowViewModel()