diff --git a/OpenNest.Core/Geometry/PolyLabel.cs b/OpenNest.Core/Geometry/PolyLabel.cs index e8de5680..a0312ff2 100644 --- a/OpenNest.Core/Geometry/PolyLabel.cs +++ b/OpenNest.Core/Geometry/PolyLabel.cs @@ -2,14 +2,37 @@ using System.Collections.Generic; namespace OpenNest.Geometry { + /// + /// Pole of inaccessibility: the interior point farthest from the outline and holes, + /// found by best-first quadtree refinement to within precision. + /// public static class PolyLabel { public static Vector Find( Polygon outer, IList holes = null, double precision = 0.5 + ) => Find(outer, holes, precision, int.MaxValue, out _); + + /// + /// Bounded search for callers that cannot afford refinement along a long thin shape: + /// stops after about 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. + /// + public static Vector Find(Polygon outer, IList holes, double precision, int maxCells) => + Find(outer, holes, precision, maxCells, out _); + + internal static Vector Find( + Polygon outer, + IList 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(); + var initial = new List(); 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(); + 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 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 holes) { var area = 0.0; diff --git a/OpenNest.Tests/Geometry/LegacyPolyLabel.cs b/OpenNest.Tests/Geometry/LegacyPolyLabel.cs new file mode 100644 index 00000000..93b2c706 --- /dev/null +++ b/OpenNest.Tests/Geometry/LegacyPolyLabel.cs @@ -0,0 +1,227 @@ +using System.Collections.Generic; +using OpenNest.Geometry; + +namespace OpenNest.Tests.Geometry +{ + /// Frozen copy of the sorted-list PolyLabel before the heap queue, used as the differential oracle. + internal static class LegacyPolyLabel + { + public static Vector Find( + Polygon outer, + IList 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(); + + 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 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 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 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 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; + } + } +} diff --git a/OpenNest.Tests/Geometry/PolyLabelDifferentialTests.cs b/OpenNest.Tests/Geometry/PolyLabelDifferentialTests.cs new file mode 100644 index 00000000..95144815 --- /dev/null +++ b/OpenNest.Tests/Geometry/PolyLabelDifferentialTests.cs @@ -0,0 +1,135 @@ +using OpenNest.Geometry; + +namespace OpenNest.Tests.Geometry; + +/// +/// The heap-queue PolyLabel must return bit-identical poles to the former sorted-list +/// implementation (frozen in ) for every unbounded call, so +/// existing label positions in PlateView and report diagrams do not move. +/// +public class PolyLabelDifferentialTests +{ + public static IEnumerable 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)> + { + (Rectangle(0, 0, 10, 10), null), + (Rectangle(0, 0, 10, 10), new List { Rectangle(4, 4, 2, 2) }), + (Rectangle(0, 0, 30, 10), new List { 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 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(); + 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); +}