perf(geometry): heap queue and optional cell budget for PolyLabel

PolyLabel kept its candidate cells in a sorted List, which costs
O(n^2) through Insert and RemoveAt(0). A 48 x 0.002 sliver searched
to 1e-4 took about 27 s.

Use a PriorityQueue keyed by (MaxDist, insertion order). Equal bounds
still pop newest first, and initial cells in sorted order, so every
existing caller gets bit-identical poles. A frozen copy of the
sorted-list implementation checks this on 1800 random shapes and on
symmetric shapes with tied poles. Swapping the tie order fails 25 of
the differential rows.

A new Find(outer, holes, precision, maxCells) overload stops after
about maxCells distance evaluations and returns the best cell so far.
When the starting grid would exceed a quarter of the budget, it is
coarsened so the whole shape is sampled first. The unbounded overload
behaves as before.

OpenNest.Core/Geometry/PolyLabel.cs

@@ -2,14 +2,37 @@ using System.Collections.Generic;
 
 namespace OpenNest.Geometry
 {
+    /// <summary>
+    /// Pole of inaccessibility: the interior point farthest from the outline and holes,
+    /// found by best-first quadtree refinement to within <c>precision</c>.
+    /// </summary>
     public static class PolyLabel
     {
         public static Vector Find(
             Polygon outer,
             IList<Polygon> holes = null,
             double precision = 0.5
+        ) => Find(outer, holes, precision, int.MaxValue, out _);
+
+        /// <summary>
+        /// Bounded search for callers that cannot afford refinement along a long thin shape:
+        /// stops after about <paramref name="maxCells"/> distance evaluations and returns the best
+        /// cell centre so far. When the starting grid would exceed a quarter of the budget it is
+        /// coarsened so the whole shape is sampled first. The result can lie outside a shape
+        /// thinner than the remaining refinement, so callers must check it.
+        /// </summary>
+        public static Vector Find(Polygon outer, IList<Polygon> holes, double precision, int maxCells) =>
+            Find(outer, holes, precision, maxCells, out _);
+
+        internal static Vector Find(
+            Polygon outer,
+            IList<Polygon> holes,
+            double precision,
+            int maxCells,
+            out int evaluatedCells
         )
         {
+            evaluatedCells = 0;
             if (outer.Vertices.Count < 3)
                 return outer.Vertices.Count > 0 ? outer.Vertices[0] : new Vector();
 
@@ -38,29 +61,43 @@ namespace OpenNest.Geometry
             if (cellSize == 0)
                 return new Vector((minX + maxX) / 2, (minY + maxY) / 2);
 
+            if (maxCells < int.MaxValue)
+            {
+                var gridLimit = System.Math.Max(1, maxCells / 4);
+                while (System.Math.Ceiling(width / cellSize) * System.Math.Ceiling(height / cellSize) > gridLimit)
+                    cellSize *= 2;
+            }
+
             var halfCell = cellSize / 2;
 
-            var queue = new List<Cell>();
+            var initial = new List<Cell>();
 
             for (var x = minX; x < maxX; x += cellSize)
-            for (var y = minY; y < maxY; y += cellSize)
-                queue.Add(new Cell(x + halfCell, y + halfCell, halfCell, outer, holes));
+                for (var y = minY; y < maxY; y += cellSize)
+                    initial.Add(new Cell(x + halfCell, y + halfCell, halfCell, outer, holes));
 
-            queue.Sort((a, b) => b.MaxDist.CompareTo(a.MaxDist));
+            initial.Sort((a, b) => b.MaxDist.CompareTo(a.MaxDist));
 
             var bestCell = GetCentroidCell(outer, holes);
+            evaluatedCells = initial.Count + 1;
 
-            for (var i = 0; i < queue.Count; i++)
-                if (queue[i].Dist > bestCell.Dist)
+            for (var i = 0; i < initial.Count; i++)
+                if (initial[i].Dist > bestCell.Dist)
                 {
-                    bestCell = queue[i];
+                    bestCell = initial[i];
                     break;
                 }
 
+            // Highest MaxDist first. Equal bounds pop newest first, and the initial cells in
+            // sorted order, exactly as the former sorted list did (insertion ahead of equals).
+            var queue = new PriorityQueue<Cell, (double, long)>();
+            for (var i = 0; i < initial.Count; i++)
+                queue.Enqueue(initial[i], (-initial[i].MaxDist, i - (long)initial.Count));
+            var sequence = (long)initial.Count;
+
             while (queue.Count > 0)
             {
-                var cell = queue[0];
-                queue.RemoveAt(0);
+                var cell = queue.Dequeue();
 
                 if (cell.Dist > bestCell.Dist)
                     bestCell = cell;
@@ -68,6 +105,9 @@ namespace OpenNest.Geometry
                 if (cell.MaxDist - bestCell.Dist <= precision)
                     continue;
 
+                if (evaluatedCells > maxCells - 4)
+                    break;
+
                 halfCell = cell.HalfSize / 2;
 
                 var newCells = new[]
@@ -77,25 +117,18 @@ namespace OpenNest.Geometry
                     new Cell(cell.X - halfCell, cell.Y + halfCell, halfCell, outer, holes),
                     new Cell(cell.X + halfCell, cell.Y + halfCell, halfCell, outer, holes),
                 };
+                evaluatedCells += newCells.Length;
 
                 for (var i = 0; i < newCells.Length; i++)
                 {
                     if (newCells[i].MaxDist > bestCell.Dist + precision)
-                        InsertSorted(queue, newCells[i]);
+                        queue.Enqueue(newCells[i], (-newCells[i].MaxDist, -++sequence));
                 }
             }
 
             return new Vector(bestCell.X, bestCell.Y);
         }
 
-        private static void InsertSorted(List<Cell> list, Cell cell)
-        {
-            var idx = 0;
-            while (idx < list.Count && list[idx].MaxDist > cell.MaxDist)
-                idx++;
-            list.Insert(idx, cell);
-        }
-
         private static Cell GetCentroidCell(Polygon outer, IList<Polygon> holes)
         {
             var area = 0.0;

OpenNest.Tests/Geometry/LegacyPolyLabel.cs

@@ -0,0 +1,227 @@
+using System.Collections.Generic;
+using OpenNest.Geometry;
+
+namespace OpenNest.Tests.Geometry
+{
+    /// <summary>Frozen copy of the sorted-list PolyLabel before the heap queue, used as the differential oracle.</summary>
+    internal static class LegacyPolyLabel
+    {
+        public static Vector Find(
+            Polygon outer,
+            IList<Polygon> holes = null,
+            double precision = 0.5
+        )
+        {
+            if (outer.Vertices.Count < 3)
+                return outer.Vertices.Count > 0 ? outer.Vertices[0] : new Vector();
+
+            var minX = double.MaxValue;
+            var minY = double.MaxValue;
+            var maxX = double.MinValue;
+            var maxY = double.MinValue;
+
+            for (var i = 0; i < outer.Vertices.Count; i++)
+            {
+                var v = outer.Vertices[i];
+                if (v.X < minX)
+                    minX = v.X;
+                if (v.Y < minY)
+                    minY = v.Y;
+                if (v.X > maxX)
+                    maxX = v.X;
+                if (v.Y > maxY)
+                    maxY = v.Y;
+            }
+
+            var width = maxX - minX;
+            var height = maxY - minY;
+            var cellSize = System.Math.Min(width, height);
+
+            if (cellSize == 0)
+                return new Vector((minX + maxX) / 2, (minY + maxY) / 2);
+
+            var halfCell = cellSize / 2;
+
+            var queue = new List<Cell>();
+
+            for (var x = minX; x < maxX; x += cellSize)
+                for (var y = minY; y < maxY; y += cellSize)
+                    queue.Add(new Cell(x + halfCell, y + halfCell, halfCell, outer, holes));
+
+            queue.Sort((a, b) => b.MaxDist.CompareTo(a.MaxDist));
+
+            var bestCell = GetCentroidCell(outer, holes);
+
+            for (var i = 0; i < queue.Count; i++)
+                if (queue[i].Dist > bestCell.Dist)
+                {
+                    bestCell = queue[i];
+                    break;
+                }
+
+            while (queue.Count > 0)
+            {
+                var cell = queue[0];
+                queue.RemoveAt(0);
+
+                if (cell.Dist > bestCell.Dist)
+                    bestCell = cell;
+
+                if (cell.MaxDist - bestCell.Dist <= precision)
+                    continue;
+
+                halfCell = cell.HalfSize / 2;
+
+                var newCells = new[]
+                {
+                    new Cell(cell.X - halfCell, cell.Y - halfCell, halfCell, outer, holes),
+                    new Cell(cell.X + halfCell, cell.Y - halfCell, halfCell, outer, holes),
+                    new Cell(cell.X - halfCell, cell.Y + halfCell, halfCell, outer, holes),
+                    new Cell(cell.X + halfCell, cell.Y + halfCell, halfCell, outer, holes),
+                };
+
+                for (var i = 0; i < newCells.Length; i++)
+                {
+                    if (newCells[i].MaxDist > bestCell.Dist + precision)
+                        InsertSorted(queue, newCells[i]);
+                }
+            }
+
+            return new Vector(bestCell.X, bestCell.Y);
+        }
+
+        private static void InsertSorted(List<Cell> list, Cell cell)
+        {
+            var idx = 0;
+            while (idx < list.Count && list[idx].MaxDist > cell.MaxDist)
+                idx++;
+            list.Insert(idx, cell);
+        }
+
+        private static Cell GetCentroidCell(Polygon outer, IList<Polygon> holes)
+        {
+            var area = 0.0;
+            var cx = 0.0;
+            var cy = 0.0;
+            var verts = outer.Vertices;
+
+            for (int i = 0, j = verts.Count - 1; i < verts.Count; j = i++)
+            {
+                var a = verts[i];
+                var b = verts[j];
+                var cross = a.X * b.Y - b.X * a.Y;
+                cx += (a.X + b.X) * cross;
+                cy += (a.Y + b.Y) * cross;
+                area += cross;
+            }
+
+            area *= 0.5;
+
+            if (System.Math.Abs(area) < 1e-10)
+                return new Cell(verts[0].X, verts[0].Y, 0, outer, holes);
+
+            cx /= (6 * area);
+            cy /= (6 * area);
+
+            return new Cell(cx, cy, 0, outer, holes);
+        }
+
+        private static double PointToPolygonDist(double x, double y, Polygon polygon)
+        {
+            var minDist = double.MaxValue;
+            var verts = polygon.Vertices;
+
+            for (int i = 0, j = verts.Count - 1; i < verts.Count; j = i++)
+            {
+                var a = verts[i];
+                var b = verts[j];
+
+                var dx = b.X - a.X;
+                var dy = b.Y - a.Y;
+
+                if (dx != 0 || dy != 0)
+                {
+                    var t = ((x - a.X) * dx + (y - a.Y) * dy) / (dx * dx + dy * dy);
+
+                    if (t > 1)
+                    {
+                        a = b;
+                    }
+                    else if (t > 0)
+                    {
+                        a = new Vector(a.X + dx * t, a.Y + dy * t);
+                    }
+                }
+
+                var segDx = x - a.X;
+                var segDy = y - a.Y;
+                var dist = System.Math.Sqrt(segDx * segDx + segDy * segDy);
+
+                if (dist < minDist)
+                    minDist = dist;
+            }
+
+            return minDist;
+        }
+
+        private sealed class Cell
+        {
+            public readonly double X;
+            public readonly double Y;
+            public readonly double HalfSize;
+            public readonly double Dist;
+            public readonly double MaxDist;
+
+            public Cell(double x, double y, double halfSize, Polygon outer, IList<Polygon> holes)
+            {
+                X = x;
+                Y = y;
+                HalfSize = halfSize;
+
+                var pt = new Vector(x, y);
+                var inside = outer.ContainsPoint(pt);
+
+                if (inside && holes != null)
+                {
+                    for (var i = 0; i < holes.Count; i++)
+                    {
+                        if (holes[i].ContainsPoint(pt))
+                        {
+                            inside = false;
+                            break;
+                        }
+                    }
+                }
+
+                Dist = PointToAllEdgesDist(x, y, outer, holes);
+
+                if (!inside)
+                    Dist = -Dist;
+
+                MaxDist = Dist + HalfSize * System.Math.Sqrt(2);
+            }
+        }
+
+        private static double PointToAllEdgesDist(
+            double x,
+            double y,
+            Polygon outer,
+            IList<Polygon> holes
+        )
+        {
+            var minDist = PointToPolygonDist(x, y, outer);
+
+            if (holes != null)
+            {
+                for (var i = 0; i < holes.Count; i++)
+                {
+                    var d = PointToPolygonDist(x, y, holes[i]);
+                    if (d < minDist)
+                        minDist = d;
+                }
+            }
+
+            return minDist;
+        }
+    }
+}

OpenNest.Tests/Geometry/PolyLabelDifferentialTests.cs

@@ -0,0 +1,135 @@
+using OpenNest.Geometry;
+
+namespace OpenNest.Tests.Geometry;
+
+/// <summary>
+/// The heap-queue PolyLabel must return bit-identical poles to the former sorted-list
+/// implementation (frozen in <see cref="LegacyPolyLabel"/>) for every unbounded call, so
+/// existing label positions in PlateView and report diagrams do not move.
+/// </summary>
+public class PolyLabelDifferentialTests
+{
+    public static IEnumerable<object[]> Seeds() => Enumerable.Range(0, 40).Select(seed => new object[] { seed });
+
+    [Theory]
+    [MemberData(nameof(Seeds))]
+    public void Find_MatchesTheSortedListImplementationBitForBit(int seed)
+    {
+        var random = new Random(seed);
+        for (var round = 0; round < 15; round++)
+        {
+            var (outer, holes) = RandomShape(random);
+            var box = Bounds(outer);
+            foreach (var precision in new[] { 0.5, System.Math.Max(box.Length, box.Width) * 0.005, 1e-3 })
+            {
+                var expected = LegacyPolyLabel.Find(outer, holes, precision);
+                var actual = PolyLabel.Find(outer, holes, precision);
+                Assert.True(Same(expected, actual),
+                    $"seed {seed} round {round} precision {precision}: legacy {expected.X:R},{expected.Y:R} "
+                    + $"heap {actual.X:R},{actual.Y:R}");
+            }
+        }
+    }
+
+    [Fact]
+    public void Find_SymmetricShapesWithTiedPolesMatchTheSortedList()
+    {
+        var shapes = new List<(Polygon, List<Polygon>)>
+        {
+            (Rectangle(0, 0, 10, 10), null),
+            (Rectangle(0, 0, 10, 10), new List<Polygon> { Rectangle(4, 4, 2, 2) }),
+            (Rectangle(0, 0, 30, 10), new List<Polygon> { Rectangle(9, 4, 2, 2), Rectangle(19, 4, 2, 2) }),
+            (Regular(12, 5), null),
+            (Regular(4, 100), null),
+            (Rectangle(-1e6, -1e6, 3, 1), null),
+        };
+        foreach (var (outer, holes) in shapes)
+            foreach (var precision in new[] { 0.5, 0.01, 1e-4 })
+                Assert.True(Same(LegacyPolyLabel.Find(outer, holes, precision), PolyLabel.Find(outer, holes, precision)));
+    }
+
+    [Fact]
+    public void Find_WithBudgetStaysWithinItOnALongSliverAndLandsInside()
+    {
+        // 48 x 0.002: the unbounded search refines along the whole flat ridge (seconds).
+        var sliver = Rectangle(0, 0, 48, 0.002);
+
+        var pole = PolyLabel.Find(sliver, null, 1e-4, 4096, out var evaluated);
+
+        Assert.InRange(evaluated, 1, 4096);
+        Assert.True(sliver.ContainsPoint(pole), $"{pole.X:R},{pole.Y:R}");
+        Assert.InRange(pole.Y, 0.0005, 0.0015);
+    }
+
+    [Fact]
+    public void Find_UnboundedOverloadMatchesTheDefault()
+    {
+        var shape = Rectangle(0, 0, 200, 10);
+        Assert.True(Same(PolyLabel.Find(shape), PolyLabel.Find(shape, null, 0.5, int.MaxValue)));
+    }
+
+    private static (Polygon outer, List<Polygon> holes) RandomShape(Random random)
+    {
+        var scale = System.Math.Pow(10, random.NextDouble() * 4 - 1);
+        var offset = new Vector((random.NextDouble() - 0.5) * 2e4, (random.NextDouble() - 0.5) * 2e4);
+        if (random.Next(3) == 0)
+        {
+            // Square with a few disjoint square holes on a 4x4 grid of slots.
+            var outer = Rectangle(offset.X, offset.Y, 4 * scale, 4 * scale);
+            var holes = new List<Polygon>();
+            foreach (var slot in Enumerable.Range(0, 16).OrderBy(_ => random.Next()).Take(random.Next(4)))
+            {
+                var size = scale * (0.2 + 0.6 * random.NextDouble());
+                holes.Add(Rectangle(offset.X + (slot % 4 + 0.1) * scale, offset.Y + (slot / 4 + 0.1) * scale,
+                    size, size));
+            }
+            return (outer, holes);
+        }
+
+        // Star-shaped polygon: sorted angles, random radii; sometimes very elongated.
+        var count = random.Next(3, 40);
+        var stretch = random.Next(4) == 0 ? 20 + random.NextDouble() * 100 : 1;
+        var angles = Enumerable.Range(0, count).Select(_ => random.NextDouble() * 2 * System.Math.PI)
+            .OrderBy(a => a).ToArray();
+        var polygon = new Polygon();
+        foreach (var angle in angles)
+        {
+            var radius = scale * (0.2 + 0.8 * random.NextDouble());
+            polygon.Vertices.Add(new Vector(offset.X + radius * System.Math.Cos(angle) * stretch,
+                offset.Y + radius * System.Math.Sin(angle)));
+        }
+        return (polygon, null);
+    }
+
+    private static Box Bounds(Polygon polygon)
+    {
+        var left = polygon.Vertices.Min(v => v.X);
+        var bottom = polygon.Vertices.Min(v => v.Y);
+        return new Box(left, bottom, polygon.Vertices.Max(v => v.X) - left, polygon.Vertices.Max(v => v.Y) - bottom);
+    }
+
+    private static Polygon Rectangle(double x, double y, double width, double height)
+    {
+        var polygon = new Polygon();
+        polygon.Vertices.Add(new Vector(x, y));
+        polygon.Vertices.Add(new Vector(x + width, y));
+        polygon.Vertices.Add(new Vector(x + width, y + height));
+        polygon.Vertices.Add(new Vector(x, y + height));
+        return polygon;
+    }
+
+    private static Polygon Regular(int sides, double radius)
+    {
+        var polygon = new Polygon();
+        for (var i = 0; i < sides; i++)
+        {
+            var angle = 2 * System.Math.PI * i / sides;
+            polygon.Vertices.Add(new Vector(radius * System.Math.Cos(angle), radius * System.Math.Sin(angle)));
+        }
+        return polygon;
+    }
+
+    private static bool Same(Vector a, Vector b) =>
+        BitConverter.DoubleToInt64Bits(a.X) == BitConverter.DoubleToInt64Bits(b.X)
+        && BitConverter.DoubleToInt64Bits(a.Y) == BitConverter.DoubleToInt64Bits(b.Y);
+}