diff --git a/benchmark/CDT.Benchmarks/Benchmarks.cs b/benchmark/CDT.Benchmarks/Benchmarks.cs index 0cf35b1..549fa0d 100644 --- a/benchmark/CDT.Benchmarks/Benchmarks.cs +++ b/benchmark/CDT.Benchmarks/Benchmarks.cs @@ -4,6 +4,7 @@ using BenchmarkDotNet.Attributes; using BenchmarkDotNet.Configs; +using BenchmarkDotNet.Diagnosers; using BenchmarkDotNet.Running; using CDT; @@ -59,6 +60,7 @@ public static (List> Vertices, List Edges) Read(string fileNam /// ~2 600 constraint edges) – the same dataset used in the C++ CDT benchmarks. /// [MemoryDiagnoser] +[EventPipeProfiler(EventPipeProfile.CpuSampling)] [GroupBenchmarksBy(BenchmarkLogicalGroupRule.ByCategory)] [CategoriesColumn] [ShortRunJob] @@ -127,6 +129,16 @@ public Triangulation Conforming_Auto() return cdt; } + [Benchmark(Description = "Conforming – AsProvided")] + [BenchmarkCategory("Conforming")] + public Triangulation Conforming_AsProvided() + { + var cdt = new Triangulation(VertexInsertionOrder.AsProvided); + cdt.InsertVertices(_vertices); + cdt.ConformToEdges(_edges); + return cdt; + } + // -- Full pipeline (insert + erase outer + holes) ------------------------ [Benchmark(Description = "Full pipeline – Auto")] @@ -140,6 +152,30 @@ public Triangulation FullPipeline_Auto() cdt.EraseOuterTrianglesAndHoles(); return cdt; } + + // -- Finalization APIs (EraseSuperTriangle / EraseOuterTriangles) --------- + + [Benchmark(Description = "EraseSuperTriangle – Auto")] + [BenchmarkCategory("Finalization")] + public Triangulation EraseSuperTriangle_Auto() + { + var cdt = new Triangulation(VertexInsertionOrder.Auto); + cdt.InsertVertices(_vertices); + cdt.InsertEdges(_edges); + cdt.EraseSuperTriangle(); + return cdt; + } + + [Benchmark(Description = "EraseOuterTriangles – Auto")] + [BenchmarkCategory("Finalization")] + public Triangulation EraseOuterTriangles_Auto() + { + var cdt = new Triangulation(VertexInsertionOrder.Auto); + cdt.InsertVertices(_vertices); + cdt.InsertEdges(_edges); + cdt.EraseOuterTriangles(); + return cdt; + } } // --------------------------------------------------------------------------- @@ -151,6 +187,7 @@ public Triangulation FullPipeline_Auto() /// to measure per-vertex overhead without the noise of large datasets. /// [MemoryDiagnoser] +[EventPipeProfiler(EventPipeProfile.CpuSampling)] [GroupBenchmarksBy(BenchmarkLogicalGroupRule.ByCategory)] [CategoriesColumn] [ShortRunJob] @@ -206,4 +243,37 @@ public Triangulation FloatVsDouble_Float() cdt.InsertEdges(ef); return cdt; } + + [Benchmark(Description = "Small – Conforming Auto")] + [BenchmarkCategory("SmallConforming")] + public Triangulation Small_Conforming_Auto() + { + var cdt = new Triangulation(VertexInsertionOrder.Auto); + cdt.InsertVertices(_vertices); + cdt.ConformToEdges(_edges); + return cdt; + } + + [Benchmark(Description = "Small – EraseSuperTriangle Auto")] + [BenchmarkCategory("SmallFinalization")] + public Triangulation Small_EraseSuperTriangle_Auto() + { + var cdt = new Triangulation(VertexInsertionOrder.Auto); + cdt.InsertVertices(_vertices); + cdt.InsertEdges(_edges); + cdt.EraseSuperTriangle(); + return cdt; + } + + [Benchmark(Description = "Small – EraseOuterTrianglesAndHoles Auto")] + [BenchmarkCategory("SmallFinalization")] + public Triangulation Small_EraseOuterTrianglesAndHoles_Auto() + { + var cdt = new Triangulation(VertexInsertionOrder.Auto, + IntersectingConstraintEdges.TryResolve, 0.0); + cdt.InsertVertices(_vertices); + cdt.InsertEdges(_edges); + cdt.EraseOuterTrianglesAndHoles(); + return cdt; + } } diff --git a/src/CDT.Core/KdTree.cs b/src/CDT.Core/KdTree.cs index e2c85cd..a124554 100644 --- a/src/CDT.Core/KdTree.cs +++ b/src/CDT.Core/KdTree.cs @@ -55,12 +55,14 @@ public NearestTask(int node, T minX, T minY, T maxX, T maxY, SplitDir dir, T dis private T _minX, _minY, _maxX, _maxY; private bool _boxInitialized; private int _size; + private readonly T _two; private NearestTask[] _stack = new NearestTask[InitialStackDepth]; /// Initializes an empty KD-tree with no bounding box pre-set. public KdTree() { + _two = T.One + T.One; _minX = T.MinValue; _minY = T.MinValue; _maxX = T.MaxValue; _maxY = T.MaxValue; _root = AddNewNode(); @@ -70,6 +72,7 @@ public KdTree() /// Initializes an empty KD-tree with a known bounding box. public KdTree(T minX, T minY, T maxX, T maxY) { + _two = T.One + T.One; _minX = minX; _minY = minY; _maxX = maxX; _maxY = maxY; _root = AddNewNode(); _boxInitialized = true; @@ -226,10 +229,9 @@ private static bool IsInsideBox(T px, T py, T minX, T minY, T maxX, T maxY) => px >= minX && px <= maxX && py >= minY && py <= maxY; [MethodImpl(MethodImplOptions.AggressiveInlining)] - private static T GetMid(T minX, T minY, T maxX, T maxY, SplitDir dir) + private T GetMid(T minX, T minY, T maxX, T maxY, SplitDir dir) { - T two = T.One + T.One; - return dir == SplitDir.X ? (minX + maxX) / two : (minY + maxY) / two; + return dir == SplitDir.X ? (minX + maxX) / _two : (minY + maxY) / _two; } [MethodImpl(MethodImplOptions.AggressiveInlining)] diff --git a/src/CDT.Core/Triangulation.cs b/src/CDT.Core/Triangulation.cs index 81b4a0b..2ae393d 100644 --- a/src/CDT.Core/Triangulation.cs +++ b/src/CDT.Core/Triangulation.cs @@ -100,6 +100,7 @@ public sealed class Triangulation private readonly VertexInsertionOrder _insertionOrder; private readonly IntersectingConstraintEdges _intersectingEdgesStrategy; private readonly T _minDistToConstraintEdge; + private readonly T _two; private SuperGeometryType _superGeomType; private int _nTargetVerts; @@ -137,6 +138,7 @@ public Triangulation( _insertionOrder = insertionOrder; _intersectingEdgesStrategy = intersectingEdgesStrategy; _minDistToConstraintEdge = minDistToConstraintEdge; + _two = T.One + T.One; _superGeomType = SuperGeometryType.SuperTriangle; _nTargetVerts = 0; _pieceToOriginalsView = new CovariantReadOnlyDictionary, IReadOnlyList>(_pieceToOriginals); @@ -153,6 +155,16 @@ public void InsertVertices(IReadOnlyList> newVertices) bool isFirstInsertion = _kdTree == null && _vertices.Count == 0; + // Pre-allocate backing arrays once we know the incoming vertex count. + // Euler's formula: a planar triangulation of N points has ~2N triangles. + if (isFirstInsertion) + { + int n = newVertices.Count; + _vertices.EnsureCapacity(n + Indices.SuperTriangleVertexCount); + _vertTris.EnsureCapacity(n + Indices.SuperTriangleVertexCount); + _triangles.EnsureCapacity(2 * n + 4); + } + // Build bounding box of new vertices var box = new Box2d(); box.Envelop(newVertices); @@ -187,9 +199,10 @@ public void InsertVertices(IReadOnlyList> newVertices) else { // AsProvided: sequential order, KD-tree walk-start + var stack = new Stack(4); for (int iV = insertStart; iV < _vertices.Count; iV++) { - InsertVertex(iV); + InsertVertex(iV, stack); } } } @@ -203,6 +216,10 @@ public void InsertEdges(IReadOnlyList edges) { var remaining = new List(4); var tppTasks = new List(8); + var polyL = new List(8); + var polyR = new List(8); + var outerTris = new Dictionary(); + var intersected = new List(8); foreach (var e in edges) { remaining.Clear(); @@ -211,7 +228,7 @@ public void InsertEdges(IReadOnlyList edges) { Edge edge = remaining[^1]; remaining.RemoveAt(remaining.Count - 1); - InsertEdgeIteration(edge, new Edge(e.V1 + _nTargetVerts, e.V2 + _nTargetVerts), remaining, tppTasks); + InsertEdgeIteration(edge, new Edge(e.V1 + _nTargetVerts, e.V2 + _nTargetVerts), remaining, tppTasks, polyL, polyR, outerTris, intersected); } } } @@ -227,6 +244,8 @@ public void InsertEdges(IReadOnlyList edges) public void ConformToEdges(IReadOnlyList edges) { var remaining = new List(8); + var flipStack = new Stack(4); + var flippedFixed = new List(8); foreach (var e in edges) { var shifted = new Edge(e.V1 + _nTargetVerts, e.V2 + _nTargetVerts); @@ -236,7 +255,7 @@ public void ConformToEdges(IReadOnlyList edges) { var task = remaining[^1]; remaining.RemoveAt(remaining.Count - 1); - ConformToEdgeIteration(task.Edge, task.Originals, task.Overlaps, remaining); + ConformToEdgeIteration(task.Edge, task.Originals, task.Overlaps, remaining, flipStack, flippedFixed); } } } @@ -299,18 +318,17 @@ private void AddSuperTriangle(Box2d box) _nTargetVerts = Indices.SuperTriangleVertexCount; _superGeomType = SuperGeometryType.SuperTriangle; - T two = T.One + T.One; - T cx = (box.Min.X + box.Max.X) / two; - T cy = (box.Min.Y + box.Max.Y) / two; + T cx = (box.Min.X + box.Max.X) / _two; + T cy = (box.Min.Y + box.Max.Y) / _two; T w = box.Max.X - box.Min.X; T h = box.Max.Y - box.Min.Y; T r = T.Max(w, h); - r = T.Max(two * r, T.One); + r = T.Max(_two * r, T.One); // Guard against very large numbers - while (cy <= cy - r) r = two * r; + while (cy <= cy - r) r = _two * r; - T R = two * r; + T R = _two * r; T cos30 = ParseT("0.8660254037844386"); T shiftX = R * cos30; @@ -352,16 +370,31 @@ private int AddTriangle(Triangle t) // Internal helpers – vertex insertion // ------------------------------------------------------------------------- - private void InsertVertex(int iVert, int walkStart) + private void InsertVertex(int iVert, int walkStart, Stack stack) { var (iT, iTopo) = WalkingSearchTrianglesAt(iVert, walkStart); - var stack = iTopo == Indices.NoNeighbor - ? InsertVertexInsideTriangle(iVert, iT) - : InsertVertexOnEdge(iVert, iT, iTopo, handleFixedSplitEdge: true); + if (iTopo == Indices.NoNeighbor) + InsertVertexInsideTriangle(iVert, iT, stack); + else + InsertVertexOnEdge(iVert, iT, iTopo, handleFixedSplitEdge: true, stack); EnsureDelaunayByEdgeFlips(iVert, stack); TryAddVertexToLocator(iVert); } + private void InsertVertex(int iVert, int walkStart) + { + var stack = new Stack(4); + InsertVertex(iVert, walkStart, stack); + } + + private void InsertVertex(int iVert, Stack stack) + { + int near = _kdTree != null + ? _kdTree.Nearest(_vertices[iVert].X, _vertices[iVert].Y, _vertices) + : 0; + InsertVertex(iVert, near, stack); + } + private void InsertVertex(int iVert) { // Walk-start from KD-tree nearest point, or vertex 0 as fallback @@ -381,7 +414,8 @@ private void InsertVertices_Randomized(int superGeomVertCount) int j = Random.Shared.Next(i + 1); (indices[i], indices[j]) = (indices[j], indices[i]); } - foreach (int iV in indices) { InsertVertex(iV); } + var stack = new Stack(4); + foreach (int iV in indices) { InsertVertex(iV, stack); } } private void InsertVertices_KDTreeBFS(int superGeomVertCount, Box2d box) @@ -395,28 +429,29 @@ private void InsertVertices_KDTreeBFS(int superGeomVertCount, Box2d box) var queue = new Queue<(int lo, int hi, T boxMinX, T boxMinY, T boxMaxX, T boxMaxY, int parent)>(); queue.Enqueue((0, vertexCount, box.Min.X, box.Min.Y, box.Max.X, box.Max.Y, 0)); + var stack = new Stack(4); while (queue.Count > 0) { var (lo, hi, boxMinX, boxMinY, boxMaxX, boxMaxY, parent) = queue.Dequeue(); int len = hi - lo; if (len == 0) { continue; } - if (len == 1) { InsertVertex(indices[lo], parent); continue; } + if (len == 1) { InsertVertex(indices[lo], parent, stack); continue; } int midPos = lo + len / 2; if (T.CreateChecked(boxMaxX - boxMinX) >= T.CreateChecked(boxMaxY - boxMinY)) { - NthElement(indices, lo, midPos, hi, (a, b) => _vertices[a].X.CompareTo(_vertices[b].X)); + NthElement(indices, lo, midPos, hi, new VertexXComparer(_vertices)); T split = _vertices[indices[midPos]].X; - InsertVertex(indices[midPos], parent); + InsertVertex(indices[midPos], parent, stack); if (lo < midPos) { queue.Enqueue((lo, midPos, boxMinX, boxMinY, split, boxMaxY, indices[midPos])); } if (midPos + 1 < hi) { queue.Enqueue((midPos + 1, hi, split, boxMinY, boxMaxX, boxMaxY, indices[midPos])); } } else { - NthElement(indices, lo, midPos, hi, (a, b) => _vertices[a].Y.CompareTo(_vertices[b].Y)); + NthElement(indices, lo, midPos, hi, new VertexYComparer(_vertices)); T split = _vertices[indices[midPos]].Y; - InsertVertex(indices[midPos], parent); + InsertVertex(indices[midPos], parent, stack); if (lo < midPos) { queue.Enqueue((lo, midPos, boxMinX, boxMinY, boxMaxX, split, indices[midPos])); } if (midPos + 1 < hi) { queue.Enqueue((midPos + 1, hi, boxMinX, split, boxMaxX, boxMaxY, indices[midPos])); } } @@ -427,12 +462,27 @@ private void InsertVertices_KDTreeBFS(int superGeomVertCount, Box2d box) // nth_element — O(n) average quickselect for spatial BFS partitioning // ------------------------------------------------------------------------- + private readonly struct VertexXComparer : IComparer + { + private readonly IReadOnlyList> _vertices; + public VertexXComparer(IReadOnlyList> vertices) => _vertices = vertices; + public int Compare(int a, int b) => _vertices[a].X.CompareTo(_vertices[b].X); + } + + private readonly struct VertexYComparer : IComparer + { + private readonly IReadOnlyList> _vertices; + public VertexYComparer(IReadOnlyList> vertices) => _vertices = vertices; + public int Compare(int a, int b) => _vertices[a].Y.CompareTo(_vertices[b].Y); + } + /// /// Rearranges in [, ) /// so the element at position is the one that would be there /// after a full sort; elements before it are ≤ it and elements after are ≥ it. /// - private static void NthElement(int[] arr, int lo, int nth, int hi, Comparison cmp) + private static void NthElement(int[] arr, int lo, int nth, int hi, TComparer cmp) + where TComparer : struct, IComparer { while (lo < hi - 1) { @@ -442,7 +492,7 @@ private static void NthElement(int[] arr, int lo, int nth, int hi, Comparison InsertVertex_FlipFixedEdges(int iV) + private void InsertVertex_FlipFixedEdges(int iV, Stack stack, List flipped) { - var flipped = new List(); + flipped.Clear(); // Use KD-tree if available, otherwise fall back to vertex 0 (first super-triangle vertex) int near = _kdTree != null ? _kdTree.Nearest(_vertices[iV].X, _vertices[iV].Y, _vertices) : 0; var (iT, iTopo) = WalkingSearchTrianglesAt(iV, near); - var stack = iTopo == Indices.NoNeighbor - ? InsertVertexInsideTriangle(iV, iT) - : InsertVertexOnEdge(iV, iT, iTopo, handleFixedSplitEdge: false); + if (iTopo == Indices.NoNeighbor) + InsertVertexInsideTriangle(iV, iT, stack); + else + InsertVertexOnEdge(iV, iT, iTopo, handleFixedSplitEdge: false, stack); int _dbgFlipIter2 = 0; while (stack.Count > 0) @@ -487,7 +538,6 @@ private List InsertVertex_FlipFixedEdges(int iV) } } TryAddVertexToLocator(iV); - return flipped; } private void EnsureDelaunayByEdgeFlips(int iV1, Stack triStack) @@ -588,7 +638,7 @@ private int WalkTriangles(int startVertex, V2d pos) // Insert vertex inside triangle or on edge // ------------------------------------------------------------------------- - private Stack InsertVertexInsideTriangle(int v, int iT) + private void InsertVertexInsideTriangle(int v, int iT, Stack stack) { int iNewT1 = AddTriangle(); int iNewT2 = AddTriangle(); @@ -606,14 +656,13 @@ private Stack InsertVertexInsideTriangle(int v, int iT) ChangeNeighbor(n2, iT, iNewT1); ChangeNeighbor(n3, iT, iNewT2); - var stack = new Stack(3); + stack.Clear(); stack.Push(iT); stack.Push(iNewT1); stack.Push(iNewT2); - return stack; } - private Stack InsertVertexOnEdge(int v, int iT1, int iT2, bool handleFixedSplitEdge) + private void InsertVertexOnEdge(int v, int iT1, int iT2, bool handleFixedSplitEdge, Stack stack) { int iTnew1 = AddTriangle(); int iTnew2 = AddTriangle(); @@ -649,12 +698,11 @@ private Stack InsertVertexOnEdge(int v, int iT1, int iT2, bool handleFixedS SplitFixedEdge(sharedEdge, v); } - var stack = new Stack(4); + stack.Clear(); stack.Push(iT1); stack.Push(iTnew2); stack.Push(iT2); stack.Push(iTnew1); - return stack; } // ------------------------------------------------------------------------- @@ -701,9 +749,11 @@ private void EdgeFlipInfo( n4 = tOpo.GetNeighbor(CdtUtils.Cw(oi)); } + [MethodImpl(MethodImplOptions.AggressiveInlining)] private bool IsFlipNeeded(int iV1, int iV2, int iV3, int iV4) { - if (_fixedEdges.Contains(new Edge(iV2, iV4))) return false; + // Skip HashSet lookup when there are no fixed edges (pure vertex-insertion path). + if (_fixedEdges.Count > 0 && _fixedEdges.Contains(new Edge(iV2, iV4))) return false; var v1 = _vertices[iV1]; var v2 = _vertices[iV2]; @@ -760,7 +810,11 @@ private void FlipEdge( private void InsertEdgeIteration( Edge edge, Edge originalEdge, List remaining, - List tppIterations) + List tppIterations, + List polyL, + List polyR, + Dictionary outerTris, + List intersected) { int iA = edge.V1, iB = edge.V2; if (iA == iB) return; @@ -786,14 +840,16 @@ private void InsertEdgeIteration( return; } - var polyL = new List(8) { iA, iVL }; - var polyR = new List(8) { iA, iVR }; - var outerTris = new Dictionary - { - [new Edge(iA, iVL)] = CdtUtils.EdgeNeighbor(_triangles[iT], iA, iVL), - [new Edge(iA, iVR)] = CdtUtils.EdgeNeighbor(_triangles[iT], iA, iVR), - }; - var intersected = new List(8) { iT }; + polyL.Clear(); + polyR.Clear(); + outerTris.Clear(); + intersected.Clear(); + + polyL.Add(iA); polyL.Add(iVL); + polyR.Add(iA); polyR.Add(iVR); + outerTris[new Edge(iA, iVL)] = CdtUtils.EdgeNeighbor(_triangles[iT], iA, iVL); + outerTris[new Edge(iA, iVR)] = CdtUtils.EdgeNeighbor(_triangles[iT], iA, iVR); + intersected.Add(iT); int iV = iA; var t = _triangles[iT]; @@ -906,7 +962,9 @@ private void HandleIntersectingEdgeStrategy( private void ConformToEdgeIteration( Edge edge, List originals, ushort overlaps, - List remaining) + List remaining, + Stack flipStack, + List flippedFixed) { int iA = edge.V1, iB = edge.V2; if (iA == iB) return; @@ -965,10 +1023,9 @@ private void ConformToEdgeIteration( int iMid = _vertices.Count; var start = _vertices[iA]; var end = _vertices[iB]; - T two = T.One + T.One; - AddNewVertex(new V2d((start.X + end.X) / two, (start.Y + end.Y) / two), Indices.NoNeighbor); + AddNewVertex(new V2d((start.X + end.X) / _two, (start.Y + end.Y) / _two), Indices.NoNeighbor); - var flippedFixed = InsertVertex_FlipFixedEdges(iMid); + InsertVertex_FlipFixedEdges(iMid, flipStack, flippedFixed); remaining.Add(new ConformToEdgeTask(new Edge(iMid, iB), originals, overlaps)); remaining.Add(new ConformToEdgeTask(new Edge(iA, iMid), originals, overlaps)); @@ -1241,7 +1298,8 @@ private int SplitFixedEdgeAt(Edge edge, V2d splitVert, int iT, int iTopo) { int iSplit = _vertices.Count; AddNewVertex(splitVert, Indices.NoNeighbor); - var stack = InsertVertexOnEdge(iSplit, iT, iTopo, handleFixedSplitEdge: false); + var stack = new Stack(4); + InsertVertexOnEdge(iSplit, iT, iTopo, handleFixedSplitEdge: false, stack); TryAddVertexToLocator(iSplit); EnsureDelaunayByEdgeFlips(iSplit, stack); SplitFixedEdge(edge, iSplit);