diff --git a/OpenNest.Core/Geometry/EdgeGridPolygon.cs b/OpenNest.Core/Geometry/EdgeGridPolygon.cs new file mode 100644 index 0000000..f3ff00b --- /dev/null +++ b/OpenNest.Core/Geometry/EdgeGridPolygon.cs @@ -0,0 +1,401 @@ +#nullable enable +using System; + +namespace OpenNest.Geometry +{ + /// + /// Immutable flat-array polygon with a uniform edge grid, used as an outer-shell + /// clearance prefilter. Two closed polygons share positive area only when an edge pair + /// crosses/touches or one polygon's vertex lies strictly inside the other; neither + /// happening certifies the two closed regions (hence any materials inside them) are + /// clear. returns Clear only in that certified case and Unknown + /// for uncertain contacts, so it can only ever skip the exact + /// gate when the exact gate would also find no overlap - the exact gate triangulates + /// both polygons per call and dominates runtime on finely flattened arc geometry. + /// + /// A holds the shared geometry; + /// produces a placement in world coordinates in O(1) - translation leaves the grid and + /// all cell indices unchanged, only the predicate coordinates shift. + /// + /// + public sealed class EdgeGridPolygon + { + /// Vertex-on-segment / collinearity tolerance for conservative touches. + private const double TouchEps = 1e-9; + + private readonly EdgeGridPolygonTemplate _template; + + /// Translation applied to the shared template geometry. + private readonly double Dx; + + private readonly double Dy; + + private EdgeGridPolygon(EdgeGridPolygonTemplate template, double dx, double dy) + { + _template = template; + Dx = dx; + Dy = dy; + } + + private double MinX => _template.MinX + Dx; + private double MinY => _template.MinY + Dy; + private double MaxX => _template.MaxX + Dx; + private double MaxY => _template.MaxY + Dy; + + /// + /// Builds from a closed (last vertex may repeat the first). + /// Returns null when the polygon has no usable ring - callers treat that as + /// "no information" and fall through to the exact gate. + /// + public static EdgeGridPolygon? From(Polygon polygon) + { + var template = EdgeGridPolygonTemplate.Build(polygon); + return template == null ? null : new EdgeGridPolygon(template, 0, 0); + } + + /// Returns a placement sharing immutable geometry, with an added translation. + public EdgeGridPolygon Translated(double dx, double dy) => new(_template, Dx + dx, Dy + dy); + + private double X(int i) => _template.X[i] + Dx; + private double Y(int i) => _template.Y[i] + Dy; + + /// + /// Certifies disjoint filled perimeters. Any crossing, containment or uncertain + /// boundary contact returns Unknown and must defer to the exact collision test. + /// Holes need not be supplied: removing material cannot invalidate Clear. + /// + public static ShellRelation Relate(EdgeGridPolygon a, EdgeGridPolygon b) + { + if ( + a.MaxX <= b.MinX + || b.MaxX <= a.MinX + || a.MaxY <= b.MinY + || b.MaxY <= a.MinY + ) + return ShellRelation.Clear; // disjoint bounding boxes + + // One walk per direction reports the strongest edge relation: a transversal + // crossing shares a positive-area wedge (overlap); a mere touch shares zero + // area but may hide a crossing in near-degenerate coordinates (unknown). + var edge = EdgeRelation(a, b); + if (edge < 2) + { + var back = EdgeRelation(b, a); + if (back > edge) + edge = back; + } + if (edge == 2) + return ShellRelation.Unknown; + + // Only fully disjoint boundaries can certify clearance. Point touches and + // collinear/near-degenerate contacts always defer to the reference test. + switch (edge) + { + case 0: + if (ContainsPointStrictly(a, b.X(0), b.Y(0))) + return ShellRelation.Unknown; + if (ContainsPointStrictly(b, a.X(0), a.Y(0))) + return ShellRelation.Unknown; + return ShellRelation.Clear; + default: + return ShellRelation.Unknown; + } + } + + /// + /// Classifies whether any edge of crosses or touches the boundary of + /// . Walks p's grid using each query edge's own bbox cells. + /// p's grid lives in p's LOCAL frame (the template's own coordinates), so the + /// query edge is converted by subtracting p's translation first. + /// + private static int EdgeRelation(EdgeGridPolygon p, EdgeGridPolygon q) + { + var t = p._template; + var n = t.Count; + Span seen = n <= 1024 ? stackalloc int[n] : new int[n]; + seen.Clear(); + var head = t.Head; + var nodeEdge = t.NodeEdge; + var nodeNext = t.NodeNext; + var no = q._template.Count; + var strongest = 0; + + for (var e = 0; e < no; e++) + { + // Stamp per QUERY edge: a grid edge may need testing against every query + // edge; the dedupe only collapses cells an individual query edge crosses + // more than once. + var stamp = e + 1; + var i2 = (e + 1) % no; + var p0x = q.X(e) - p.Dx; + var p0y = q.Y(e) - p.Dy; + var p1x = q.X(i2) - p.Dx; + var p1y = q.Y(i2) - p.Dy; + + var c0 = ColLow(t, p0x, p1x); + if (c0 > ColHigh(t, p0x, p1x)) + continue; + var c1 = ColHigh(t, p0x, p1x); + var r0 = RowLow(t, p0y, p1y); + if (r0 > RowHigh(t, p0y, p1y)) + continue; + var r1 = RowHigh(t, p0y, p1y); + + for (var r = r0; r <= r1; r++) + for (var c = c0; c <= c1; c++) + for (var nIdx = head[r * t.Cols + c]; nIdx >= 0; nIdx = nodeNext[nIdx]) + { + var ea = nodeEdge[nIdx]; + if (seen[ea] == stamp) + continue; + seen[ea] = stamp; + var a2 = (ea + 1) % n; + var relation = SegmentRelation( + t.X[ea], t.Y[ea], t.X[a2], t.Y[a2], p0x, p0y, p1x, p1y + ); + if (relation == 2) + return 2; // transversal crossing + if (relation > strongest) + strongest = relation; + } + } + return strongest; + } + + /// + /// Segment-pair relation: 2 = transversal crossing (strict sign flips on both + /// orientations - the regions share a positive-area wedge); 1 = a clean endpoint + /// touch (zero shared area by itself; callers decide via interior-vertex tests); + /// 3 = collinear or near-degenerate contact (a shared boundary segment can hide + /// either a same-side positive overlap or an opposite-side tangency, so it must + /// defer to the exact gate); 0 = disjoint. + /// + private static int SegmentRelation( + double ax, double ay, double bx, double by, double cx, double cy, double dx, double dy + ) + { + var rx = bx - ax; + var ry = by - ay; + var sx = dx - cx; + var sy = dy - cy; + var d1 = rx * (cy - ay) - ry * (cx - ax); + var d2 = rx * (dy - ay) - ry * (dx - ax); + var d3 = sx * (ay - cy) - sy * (ax - cx); + var d4 = sx * (by - cy) - sy * (bx - cx); + + if (((d1 > 0 && d2 < 0) || (d1 < 0 && d2 > 0)) && ((d3 > 0 && d4 < 0) || (d3 < 0 && d4 > 0))) + return 2; // proper crossing + + // A near-zero orientation means the configuration is collinear or too close to + // classify; only exact-zero orientations get the clean point-touch verdict. + var scale = System.Math.Max( + 1e-30, + System.Math.Max(System.Math.Abs(rx) + System.Math.Abs(ry), System.Math.Abs(sx) + System.Math.Abs(sy)) + ); + var eps = TouchEps * scale; + var nearDegenerate = + (System.Math.Abs(d1) <= eps && d1 != 0) + || (System.Math.Abs(d2) <= eps && d2 != 0) + || (System.Math.Abs(d3) <= eps && d3 != 0) + || (System.Math.Abs(d4) <= eps && d4 != 0); + var exactDegenerate = d1 == 0 || d2 == 0 || d3 == 0 || d4 == 0; + + var touch = + (d1 == 0 && PointOnSegment(cx, cy, ax, ay, bx, by)) + || (d2 == 0 && PointOnSegment(dx, dy, ax, ay, bx, by)) + || (d3 == 0 && PointOnSegment(ax, ay, cx, cy, dx, dy)) + || (d4 == 0 && PointOnSegment(bx, by, cx, cy, dx, dy)); + + if (nearDegenerate) + return 3; + if (exactDegenerate) + // Collinear: contact along a segment (or too close to tell) must defer to + // the exact gate; collinear but disjoint edges simply do not touch. + return touch ? 3 : 0; + if (touch) + return 1; + return 0; + } + + private static bool PointOnSegment( + double px, double py, double ax, double ay, double bx, double by + ) => + System.Math.Min(ax, bx) - TouchEps <= px + && px <= System.Math.Max(ax, bx) + TouchEps + && System.Math.Min(ay, by) - TouchEps <= py + && py <= System.Math.Max(ay, by) + TouchEps; + + /// Strict ray-cast containment (boundary touches are excluded upstream). + private static bool ContainsPointStrictly(EdgeGridPolygon poly, double px, double py) + { + var t = poly._template; + var inside = false; + var n = t.Count; + for (var i = 0; i < n; i++) + { + var j = (i + 1) % n; + var yi = poly.Y(i); + var yj = poly.Y(j); + if ((yi > py) != (yj > py)) + { + var xAt = poly.X(i) + (py - yi) / (yj - yi) * (poly.X(j) - poly.X(i)); + if (px < xAt) + inside = !inside; + } + } + return inside; + } + + private static int ColLow(EdgeGridPolygonTemplate t, double a, double b) => + System.Math.Clamp((int)System.Math.Floor((System.Math.Min(a, b) - t.MinX) / t.CellSize), 0, t.Cols); + + private static int ColHigh(EdgeGridPolygonTemplate t, double a, double b) => + System.Math.Clamp((int)System.Math.Floor((System.Math.Max(a, b) - t.MinX) / t.CellSize), -1, t.Cols - 1); + + private static int RowLow(EdgeGridPolygonTemplate t, double a, double b) => + System.Math.Clamp((int)System.Math.Floor((System.Math.Min(a, b) - t.MinY) / t.CellSize), 0, t.Rows); + + private static int RowHigh(EdgeGridPolygonTemplate t, double a, double b) => + System.Math.Clamp((int)System.Math.Floor((System.Math.Max(a, b) - t.MinY) / t.CellSize), -1, t.Rows - 1); + + /// + /// Shared, immutable grid geometry for ; the grid is defined + /// relative to the shape's own local coordinates, so translated instances reuse it. + /// Per-query deduplication scratch is local, so placements may be queried concurrently. + /// + private sealed class EdgeGridPolygonTemplate + { + internal readonly double[] X; + internal readonly double[] Y; + internal readonly int Count; + internal readonly double MinX; + internal readonly double MinY; + internal readonly double MaxX; + internal readonly double MaxY; + + internal readonly double CellSize; + + internal readonly int Cols; + internal readonly int Rows; + internal readonly int[] Head; + + /// + /// Grid nodes as parallel (edge, next) arrays: an edge spanning several cells gets + /// one node PER cell - a single next-per-edge chain would corrupt the other cells' + /// chains and silently drop edges from the walk. + /// + internal readonly int[] NodeEdge; + + internal readonly int[] NodeNext; + + + private EdgeGridPolygonTemplate( + double[] x, + double[] y, + int count, + double minX, + double minY, + double maxX, + double maxY + ) + { + X = x; + Y = y; + Count = count; + MinX = minX; + MinY = minY; + MaxX = maxX; + MaxY = maxY; + + var extentX = System.Math.Max(maxX - minX, 1e-9); + var extentY = System.Math.Max(maxY - minY, 1e-9); + CellSize = System.Math.Max(System.Math.Max(extentX, extentY) / 16.0, 1e-9); + Cols = System.Math.Clamp((int)System.Math.Ceiling(extentX / CellSize) + 1, 1, 48); + Rows = System.Math.Clamp((int)System.Math.Ceiling(extentY / CellSize) + 1, 1, 48); + Head = new int[Cols * Rows]; + Array.Fill(Head, -1); + + // Pass 1: count nodes; pass 2: fill (edge, next) node arrays. + var cellsPerEdge = new int[count]; + var total = 0; + for (var e = 0; e < count; e++) + { + var i2 = (e + 1) % count; + var c0 = ClampCol(System.Math.Min(x[e], x[i2]) - minX); + var c1 = ClampCol(System.Math.Max(x[e], x[i2]) - minX); + var r0 = ClampRow(System.Math.Min(y[e], y[i2]) - minY); + var r1 = ClampRow(System.Math.Max(y[e], y[i2]) - minY); + cellsPerEdge[e] = (c1 - c0 + 1) * (r1 - r0 + 1); + total += cellsPerEdge[e]; + } + NodeEdge = new int[total]; + NodeNext = new int[total]; + var node = 0; + for (var e = 0; e < count; e++) + { + var i2 = (e + 1) % count; + var c0 = ClampCol(System.Math.Min(x[e], x[i2]) - minX); + var c1 = ClampCol(System.Math.Max(x[e], x[i2]) - minX); + var r0 = ClampRow(System.Math.Min(y[e], y[i2]) - minY); + var r1 = ClampRow(System.Math.Max(y[e], y[i2]) - minY); + for (var r = r0; r <= r1; r++) + for (var c = c0; c <= c1; c++) + { + var cell = r * Cols + c; + NodeEdge[node] = e; + NodeNext[node] = Head[cell]; + Head[cell] = node; + node++; + } + } + } + + private int ClampCol(double dx) => + System.Math.Clamp((int)System.Math.Floor(dx / CellSize), 0, Cols - 1); + + private int ClampRow(double dy) => + System.Math.Clamp((int)System.Math.Floor(dy / CellSize), 0, Rows - 1); + + internal static EdgeGridPolygonTemplate? Build(Polygon polygon) + { + var vertices = polygon.Vertices; + var n = vertices.Count; + if (n >= 2 && vertices[0].X == vertices[n - 1].X && vertices[0].Y == vertices[n - 1].Y) + n--; + if (n < 3) + return null; + var xs = new double[n]; + var ys = new double[n]; + var minX = double.MaxValue; + var minY = double.MaxValue; + var maxX = double.MinValue; + var maxY = double.MinValue; + for (var i = 0; i < n; i++) + { + var vx = vertices[i].X; + var vy = vertices[i].Y; + xs[i] = vx; + ys[i] = vy; + if (vx < minX) + minX = vx; + if (vx > maxX) + maxX = vx; + if (vy < minY) + minY = vy; + if (vy > maxY) + maxY = vy; + } + return new EdgeGridPolygonTemplate(xs, ys, n, minX, minY, maxX, maxY); + } + } + } + + /// Conservative result of an outer-perimeter prefilter. + public enum ShellRelation + { + /// Filled perimeters, and therefore their material, are disjoint. + Clear, + /// Run an exact collision test; the prefilter cannot certify clearance. + Unknown, + } +} diff --git a/OpenNest.Core/Geometry/TriangulatedRegion.cs b/OpenNest.Core/Geometry/TriangulatedRegion.cs new file mode 100644 index 0000000..42343ec --- /dev/null +++ b/OpenNest.Core/Geometry/TriangulatedRegion.cs @@ -0,0 +1,541 @@ +#nullable enable +using System; +using System.Collections.Generic; + +namespace OpenNest.Geometry +{ + /// + /// Immutable triangulation of a simple, closed, lines-only perimeter and its holes. + /// Cached triangles use the reference Collision clipping and hole-subtraction rules. + /// Translation is a parameter; preparation never retains mutable input polygons. + /// Scratch arrays and hole-piece lists are allocated per query, with a bounded + /// thread-local buffer pool. Null means the caller must use Collision.HasOverlap. + /// + public sealed class TriangulatedRegion + { + // Flat vertex pool (local frame) and triangle index triples (CCW). + private readonly double[] X; + private readonly double[] Y; + + private readonly int[] _ia; + private readonly int[] _ib; + private readonly int[] _ic; + private readonly double[] _tMinX; + private readonly double[] _tMinY; + private readonly double[] _tMaxX; + private readonly double[] _tMaxY; + + private double MinX { get; } + private double MinY { get; } + private double MaxX { get; } + private double MaxY { get; } + + /// Triangulated holes in the same local frame (null when none). + private readonly TriangulatedRegion?[]? Holes; + + // Scratch bound: clipped convex pieces stay small; anything larger bails. + private const int MaxClipVertices = 48; + private const int MaxPieces = 2048; + + private TriangulatedRegion( + double[] x, + double[] y, + int[] ia, + int[] ib, + int[] ic, + double[] tMinX, + double[] tMinY, + double[] tMaxX, + double[] tMaxY, + TriangulatedRegion?[]? holes + ) + { + X = x; + Y = y; + _ia = ia; + _ib = ib; + _ic = ic; + _tMinX = tMinX; + _tMinY = tMinY; + _tMaxX = tMaxX; + _tMaxY = tMaxY; + Holes = holes; + + var minX = double.MaxValue; + var minY = double.MaxValue; + var maxX = double.MinValue; + var maxY = double.MinValue; + for (var i = 0; i < x.Length; i++) + { + if (x[i] < minX) + minX = x[i]; + if (x[i] > maxX) + maxX = x[i]; + if (y[i] < minY) + minY = y[i]; + if (y[i] > maxY) + maxY = y[i]; + } + MinX = minX; + MinY = minY; + MaxX = maxX; + MaxY = maxY; + } + + /// + /// Ear-clips a polygon ring into cached triangles. Returns null when + /// triangulation yields nothing usable - the caller falls back to Polygon gates. + /// + public static TriangulatedRegion? Build(Polygon perimeter, IReadOnlyList? holes = null) + { + try + { + var tris = ConvexDecomposition.Triangulate(perimeter); + var count = tris.Count; + if (count == 0) + return null; + + var xs = new double[count * 3]; + var ys = new double[count * 3]; + var ia = new int[count]; + var ib = new int[count]; + var ic = new int[count]; + var minXA = new double[count]; + var minYA = new double[count]; + var maxXA = new double[count]; + var maxYA = new double[count]; + + var k = 0; + for (var t = 0; t < count; t++) + { + var v = tris[t].Vertices; // closed: prev, curr, next, prev + ia[t] = k; + xs[k] = v[0].X; + ys[k] = v[0].Y; + k++; + ib[t] = k; + xs[k] = v[1].X; + ys[k] = v[1].Y; + k++; + ic[t] = k; + xs[k] = v[2].X; + ys[k] = v[2].Y; + k++; + minXA[t] = System.Math.Min(v[0].X, System.Math.Min(v[1].X, v[2].X)); + minYA[t] = System.Math.Min(v[0].Y, System.Math.Min(v[1].Y, v[2].Y)); + maxXA[t] = System.Math.Max(v[0].X, System.Math.Max(v[1].X, v[2].X)); + maxYA[t] = System.Math.Max(v[0].Y, System.Math.Max(v[1].Y, v[2].Y)); + } + + TriangulatedRegion[]? holeSets = null; + if (holes != null && holes.Count > 0) + { + holeSets = new TriangulatedRegion[holes.Count]; + for (var h = 0; h < holes.Count; h++) + { + var holeTris = ConvexDecomposition.Triangulate(holes[h]); + if (holeTris.Count == 0) + continue; + var hx = new double[holeTris.Count * 3]; + var hy = new double[holeTris.Count * 3]; + var hia = new int[holeTris.Count]; + var hib = new int[holeTris.Count]; + var hic = new int[holeTris.Count]; + var hminX = new double[holeTris.Count]; + var hminY = new double[holeTris.Count]; + var hmaxX = new double[holeTris.Count]; + var hmaxY = new double[holeTris.Count]; + var hk = 0; + for (var t = 0; t < holeTris.Count; t++) + { + var v = holeTris[t].Vertices; + hia[t] = hk; + hx[hk] = v[0].X; + hy[hk] = v[0].Y; + hk++; + hib[t] = hk; + hx[hk] = v[1].X; + hy[hk] = v[1].Y; + hk++; + hic[t] = hk; + hx[hk] = v[2].X; + hy[hk] = v[2].Y; + hk++; + hminX[t] = System.Math.Min(v[0].X, System.Math.Min(v[1].X, v[2].X)); + hminY[t] = System.Math.Min(v[0].Y, System.Math.Min(v[1].Y, v[2].Y)); + hmaxX[t] = System.Math.Max(v[0].X, System.Math.Max(v[1].X, v[2].X)); + hmaxY[t] = System.Math.Max(v[0].Y, System.Math.Max(v[1].Y, v[2].Y)); + } + holeSets[h] = new TriangulatedRegion(hx, hy, hia, hib, hic, hminX, hminY, hmaxX, hmaxY, null); + } + } + + return new TriangulatedRegion(xs, ys, ia, ib, ic, minXA, minYA, maxXA, maxYA, holeSets); + } + catch (Exception) + { + return null; + } + } + + /// + /// Positive shared area (surviving both polygons' hole sets) between this + /// translated by (adx, ady) and other translated by (bdx, bdy). Returns null + /// when the scratch bounds are exceeded and the question cannot be decided. + /// Inputs and translations must have finite coordinates. + /// + public bool? Overlaps(TriangulatedRegion other, double adx, double ady, double bdx, double bdy) + { + // Same bbox rule as Collision.BoundingBoxesOverlap: overlap must exceed + // Tolerance.Epsilon on both axes, so a hairline box overlap never reaches the + // clip stage. + var eps = OpenNest.Math.Tolerance.Epsilon; + var overlapX = + System.Math.Min(MaxX + adx, other.MaxX + bdx) - System.Math.Max(MinX + adx, other.MinX + bdx); + var overlapY = + System.Math.Min(MaxY + ady, other.MaxY + bdy) - System.Math.Max(MinY + ady, other.MinY + bdy); + if (overlapX <= eps || overlapY <= eps) + return false; + + var areaFloor = 2 * OpenNest.Math.Tolerance.Epsilon; + var clipA = new double[MaxClipVertices * 2]; + var clipB = new double[MaxClipVertices * 2]; + var piece = new double[MaxClipVertices * 2]; + + for (var ta = 0; ta < _ia.Length; ta++) + { + var aMinX = _tMinX[ta] + adx; + var aMaxX = _tMaxX[ta] + adx; + var aMinY = _tMinY[ta] + ady; + var aMaxY = _tMaxY[ta] + ady; + for (var tb = 0; tb < other._ia.Length; tb++) + { + var bMinX = other._tMinX[tb] + bdx; + var bMaxX = other._tMaxX[tb] + bdx; + var bMinY = other._tMinY[tb] + bdy; + var bMaxY = other._tMaxY[tb] + bdy; + if ( + System.Math.Min(aMaxX, bMaxX) - System.Math.Max(aMinX, bMinX) <= eps + || System.Math.Min(aMaxY, bMaxY) - System.Math.Max(aMinY, bMinY) <= eps + ) + continue; + + var count = ClipTriangle( + ta, adx, ady, other, tb, bdx, bdy, clipA, clipB, piece + ); + if (count >= MaxClipVertices) + return null; + if (count < 3) + continue; + if (TwiceArea(piece, count) <= areaFloor) + continue; + + var (hasHoles, undecided, survived) = SubtractAllHoles( + other, adx, ady, bdx, bdy, piece, count, areaFloor + ); + if (undecided) + return null; + if (hasHoles) + { + if (survived) + return true; + } + else + { + return true; // no holes on either side: the clipped region is overlap + } + } + } + return false; + } + + /// + /// Subtracts both polygons' hole triangles from one clipped region, mirroring + /// Collision.SubtractHoles: for every hole triangle, every surviving piece is + /// split per edge into outside pieces (survivors) and the inside remainder + /// (consumed). True means a positive-area piece survived ALL holes. + /// + [ThreadStatic] + private static List? s_pool; + + [ThreadStatic] + private static double[]? s_tmpA; + + [ThreadStatic] + private static double[]? s_tmpB; + + private static double[] AcquireBuffer() + { + var pool = s_pool ??= new List(); + var n = pool.Count; + if (n == 0) + return new double[MaxClipVertices * 2]; + var buf = pool[n - 1]; + pool.RemoveAt(n - 1); + return buf; + } + + private static void ReleaseBuffer(double[] buf) + { + var pool = s_pool ??= new List(); + if (pool.Count < 64) + pool.Add(buf); + } + + private static (double[] Tmp, double[] Inside) ScratchPair() + { + s_tmpA ??= new double[MaxClipVertices * 2]; + s_tmpB ??= new double[MaxClipVertices * 2]; + return (s_tmpA, s_tmpB); + } + + private (bool hasHoles, bool undecided, bool survived) SubtractAllHoles( + TriangulatedRegion other, + double adx, + double ady, + double bdx, + double bdy, + double[] piece, + int count, + double areaFloor + ) + { + var allHoles = 0; + if (Holes != null) + allHoles += Holes.Length; + if (other.Holes != null) + allHoles += other.Holes.Length; + if (allHoles == 0) + return (false, false, false); + + // pieces[0] is the caller's own buffer - never release it back to the pool. + var pieces = new List<(double[] Buf, int Count)> { (piece, count) }; + var owned = new HashSet(); + + bool SubtractOwner(TriangulatedRegion owner, double odx, double ody) + { + if (owner.Holes == null) + return true; + for (var h = 0; h < owner.Holes.Length && pieces.Count > 0; h++) + { + var hole = owner.Holes[h]; + if (hole == null) + continue; // untriangulatable hole: nothing to subtract + for (var t = 0; t < hole._ia.Length && pieces.Count > 0; t++) + { + var hMinX = hole._tMinX[t] + odx; + var hMaxX = hole._tMaxX[t] + odx; + var hMinY = hole._tMinY[t] + ody; + var hMaxY = hole._tMaxY[t] + ody; + + var next = new List<(double[], int)>(); + for (var p = 0; p < pieces.Count; p++) + { + var (buf, pc) = pieces[p]; + + // Piece bbox (built-in uses <=: touching skips subtraction). + var pMinX = double.MaxValue; + var pMinY = double.MaxValue; + var pMaxX = double.MinValue; + var pMaxY = double.MinValue; + for (var v = 0; v < pc; v++) + { + var px = buf[v * 2]; + var py = buf[v * 2 + 1]; + if (px < pMinX) + pMinX = px; + if (px > pMaxX) + pMaxX = px; + if (py < pMinY) + pMinY = py; + if (py > pMaxY) + pMaxY = py; + } + if (pMaxX <= hMinX || hMaxX <= pMinX || pMaxY <= hMinY || hMaxY <= pMinY) + { + if (next.Count >= MaxPieces) + return false; + next.Add((buf, pc)); + continue; + } + + // Clip the piece against the hole triangle's three edges: the + // outside of each edge survives as its own piece; the inside + // remainder continues into the next edge. The remainder inside + // all three edges is consumed (the hole ate it). + var rem = AcquireBuffer(); + owned.Add(rem); + Array.Copy(buf, rem, pc * 2); + var remCount = pc; + var (tmp, insideBuf) = ScratchPair(); + for (var e = 0; e < 3 && remCount >= 3; e++) + { + var ei = e == 0 ? hole._ia[t] : e == 1 ? hole._ib[t] : hole._ic[t]; + var ej = e == 0 ? hole._ib[t] : e == 1 ? hole._ic[t] : hole._ia[t]; + var sx = hole.X[ei] + odx; + var sy = hole.Y[ei] + ody; + var ex = hole.X[ej] + odx; + var ey = hole.Y[ej] + ody; + + var outCount = + ClipHalfSpace(rem, remCount, sx, sy, ex, ey, false, tmp); + if (outCount >= MaxClipVertices) + return false; + if (outCount >= 3 && TwiceArea(tmp, outCount) > areaFloor) + { + if (next.Count >= MaxPieces) + return false; // undecided + var keep = AcquireBuffer(); + owned.Add(keep); + Array.Copy(tmp, keep, outCount * 2); + next.Add((keep, outCount)); + } + remCount = + ClipHalfSpace(rem, remCount, sx, sy, ex, ey, true, insideBuf); + if (remCount >= MaxClipVertices) + return false; // undecided + Array.Copy(insideBuf, rem, remCount * 2); + } + // The inside-all-edges remainder is consumed by the hole: drop it. + owned.Remove(rem); + ReleaseBuffer(rem); + if (owned.Remove(buf)) + ReleaseBuffer(buf); + } + pieces = next; + } + } + return true; + } + + try + { + if (!SubtractOwner(this, adx, ady) || !SubtractOwner(other, bdx, bdy)) + return (true, true, false); + + foreach (var (buf, pc) in pieces) + if (pc >= 3 && TwiceArea(buf, pc) > areaFloor) + return (true, false, true); + return (true, false, false); + } + finally + { + foreach (var buffer in owned) + ReleaseBuffer(buffer); + } + } + + /// Clip this' triangle against other's triangle; returns count into piece. + private int ClipTriangle( + int ta, + double adx, + double ady, + TriangulatedRegion other, + int tb, + double bdx, + double bdy, + double[] bufA, + double[] bufB, + double[] piece + ) + { + var ia = _ia[ta]; + var ib = _ib[ta]; + var ic = _ic[ta]; + bufA[0] = X[ia] + adx; + bufA[1] = Y[ia] + ady; + bufA[2] = X[ib] + adx; + bufA[3] = Y[ib] + ady; + bufA[4] = X[ic] + adx; + bufA[5] = Y[ic] + ady; + var count = 3; + + for (var e = 0; e < 3 && count >= 3; e++) + { + var ei = e == 0 ? other._ia[tb] : e == 1 ? other._ib[tb] : other._ic[tb]; + var ej = e == 0 ? other._ib[tb] : e == 1 ? other._ic[tb] : other._ia[tb]; + var sx = other.X[ei] + bdx; + var sy = other.Y[ei] + bdy; + var ex = other.X[ej] + bdx; + var ey = other.Y[ej] + bdy; + count = ClipHalfSpace(bufA, count, sx, sy, ex, ey, true, bufB); + if (count >= MaxClipVertices) + return count; + for (var v = 0; v < count * 2; v++) + bufA[v] = bufB[v]; + } + for (var v = 0; v < System.Math.Min(count, MaxClipVertices) * 2; v++) + piece[v] = bufA[v]; + return count; + } + + /// + /// Sutherland-Hodgman clip against one directed edge's half-plane; identical + /// classification, interpolation and dedupe to Collision.ClipHalfSpace. + /// + private static int ClipHalfSpace( + double[] verts, + int count, + double sx, + double sy, + double ex, + double ey, + bool inside, + double[] outBuf + ) + { + var kept = 0; + var cap = outBuf.Length / 2; + var edgeX = ex - sx; + var edgeY = ey - sy; + for (var i = 0; i < count; i++) + { + var j = (i + 1) % count; + var cx = verts[i * 2]; + var cy = verts[i * 2 + 1]; + var nx = verts[j * 2]; + var ny = verts[j * 2 + 1]; + var cd = edgeX * (cy - sy) - edgeY * (cx - sx); + var nd = edgeX * (ny - sy) - edgeY * (nx - sx); + if (inside ? cd >= 0 : cd <= 0) + { + if (kept >= cap) + return cap; // overflow: caller treats as undecided + kept = AddDistinct(outBuf, kept, cx, cy); + } + if ((cd < 0 && nd > 0) || (cd > 0 && nd < 0)) + { + if (kept >= cap) + return cap; // overflow + var t = cd / (cd - nd); + kept = AddDistinct( + outBuf, kept, cx + t * (nx - cx), cy + t * (ny - cy) + ); + } + } + if (kept > 1 && outBuf[0] == outBuf[(kept - 1) * 2] && outBuf[1] == outBuf[(kept - 1) * 2 + 1]) + kept--; + return kept; + } + + private static int AddDistinct(double[] buf, int count, double x, double y) + { + if (count > 0 && buf[(count - 1) * 2] == x && buf[(count - 1) * 2 + 1] == y) + return count; + buf[count * 2] = x; + buf[count * 2 + 1] = y; + return count + 1; + } + + /// Twice the area, relative to vertex 0 (cancellation-safe). + private static double TwiceArea(double[] verts, int count) + { + var twiceArea = 0.0; + for (var i = 1; i + 1 < count; i++) + twiceArea += + (verts[i * 2] - verts[0]) * (verts[(i + 1) * 2 + 1] - verts[1]) + - (verts[i * 2 + 1] - verts[1]) * (verts[(i + 1) * 2] - verts[0]); + return System.Math.Abs(twiceArea); + } + } +} diff --git a/OpenNest.Tests/Geometry/CachedCollisionEquivalenceTests.cs b/OpenNest.Tests/Geometry/CachedCollisionEquivalenceTests.cs new file mode 100644 index 0000000..fea9ccd --- /dev/null +++ b/OpenNest.Tests/Geometry/CachedCollisionEquivalenceTests.cs @@ -0,0 +1,132 @@ +using System.Diagnostics; +using OpenNest.Geometry; +using Xunit.Abstractions; +using static OpenNest.Tests.Geometry.NoFitPolygonTests; + +namespace OpenNest.Tests.Geometry; + +public class CachedCollisionEquivalenceTests +{ + private readonly ITestOutputHelper output; + + public CachedCollisionEquivalenceTests(ITestOutputHelper output) => this.output = output; + + [Fact] + public void SeededPairsMatchReference() + { + var random = new Random(73862026); + var decisions = 0; + var mismatches = 0; + var fallback = 0; + var falseClear = 0; + var clear = 0; + var overlaps = 0; + var watch = Stopwatch.StartNew(); + for (var pair = 0; pair < 200; pair++) + { + var a = Make(random, pair % 5); + var b = Make(random, (pair / 5) % 5); + var holesA = pair % 4 == 0 ? new List { Move(Square(0.7), 0.2, 0.2) } : null; + var holesB = pair % 7 == 0 ? new List { Move(Square(0.6), 0.3, 0.3) } : null; + var ta = TriangulatedRegion.Build(a, holesA); + var tb = TriangulatedRegion.Build(b, holesB); + var ga = EdgeGridPolygon.From(a); + var gb = EdgeGridPolygon.From(b); + Assert.NotNull(ta); + Assert.NotNull(tb); + Assert.NotNull(ga); + Assert.NotNull(gb); + for (var sample = 0; sample < 500; sample++) + { + var ax = random.NextDouble() * 200 - 100; + var ay = random.NextDouble() * 200 - 100; + var bx = ax + random.NextDouble() * 12 - 6; + var by = ay + random.NextDouble() * 12 - 6; + if (sample % 4 == 0) + { + // Exact, near-touching, and thin positive-area contacts at both box edges. + var gap = new[] { 0, -1e-7, 1e-7, -1e-5, 1e-5, -1e-4, 1e-4 }[(sample / 4) % 7]; + bx = ax + a.BoundingBox.Right - b.BoundingBox.Left + gap; + by = ay + a.BoundingBox.Bottom - b.BoundingBox.Bottom; + } + var reference = Collision.HasOverlap(Move(a, ax, ay), Move(b, bx, by), + holesA?.Select(h => Move(h, ax, ay)).ToList(), + holesB?.Select(h => Move(h, bx, by)).ToList()); + var cached = ta.Overlaps(tb, ax, ay, bx, by); + var relation = EdgeGridPolygon.Relate(ga.Translated(ax, ay), gb.Translated(bx, by)); + if (!cached.HasValue) + fallback++; + else if (cached.Value != reference) + { + if (mismatches < 5) + output.WriteLine($"Mismatch pair={pair} sample={sample} a=({ax:R},{ay:R}) b=({bx:R},{by:R}) expected={reference}"); + mismatches++; + } + if (relation == ShellRelation.Clear) + { + clear++; + if (reference) + falseClear++; + } + if (reference) + overlaps++; + decisions++; + } + } + output.WriteLine($"Decisions={decisions}; mismatches={mismatches}; false Clear={falseClear}; " + + $"fallback={fallback} ({100.0 * fallback / decisions:F4}%); Clear={clear}; overlaps={overlaps}; elapsed={watch.Elapsed.TotalSeconds:F3}s"); + Assert.Equal(100000, decisions); + Assert.Equal(0, mismatches); + Assert.Equal(0, falseClear); + Assert.True(overlaps > 10000); + } + + [Fact] + public void HoleContainmentAndSharedGridQueries() + { + var a = Square(10); + var b = Square(1); + var ta = TriangulatedRegion.Build(a, new[] { Move(Square(6), 2, 2) }); + var tb = TriangulatedRegion.Build(b); + Assert.NotNull(ta); + Assert.NotNull(tb); + Assert.False(ta.Overlaps(tb, 0, 0, 4, 4)); + Assert.True(ta.Overlaps(tb, 0, 0, 1.5, 4)); + var ga = EdgeGridPolygon.From(a)!; + var gb = EdgeGridPolygon.From(b)!; + Parallel.For(0, 1000, i => + { + var dx = i % 2 == 0 ? 4 : 12; + Assert.Equal(dx == 4 ? ShellRelation.Unknown : ShellRelation.Clear, + EdgeGridPolygon.Relate(ga, gb.Translated(dx, 4))); + }); + } + + [Fact] + public void EmptyRingsRequestFallback() + { + Assert.Null(TriangulatedRegion.Build(new Polygon())); + Assert.Null(EdgeGridPolygon.From(new Polygon())); + } + + private static Polygon Make(Random random, int kind) + { + var size = 2 + random.NextDouble() * 2; + switch (kind) + { + case 0: + return Square(size); + case 1: + return Star(random); + case 2: + return Ring((0, 0), (size, 0), (size, 1), (1, 1), (1, size), (0, size)); + case 3: + return Ring((0, 0), (size, 0), (0, size)); + default: + var shape = new Shape(); + shape.Entities.Add(new Arc(0, 0, size / 2, 0, System.Math.PI)); + shape.Entities.Add(new Arc(0, 0, size / 2, System.Math.PI, 2 * System.Math.PI)); + return ClipperBridge.Flatten(shape, 0.08, circumscribe: false); + } + } +}