diff --git a/docs/vision-model.md b/docs/vision-model.md index 2ec91974..f780a87f 100644 --- a/docs/vision-model.md +++ b/docs/vision-model.md @@ -773,6 +773,9 @@ Changes, all exact (the polygon is bitwise the same): it; the polygon is assembled serially afterwards in the original order. Workers spin briefly between the runs of one query and sleep between frames. `ICARUS_HEIGHT_THREADS` overrides the worker count for diagnosis. + (Removed 2026-10-05: the query now files the edges near the eye by angle + once and casts each ray against its own bin, so casting is a tenth of the + query and runs on the calling thread.) * Event generation culls vertices outside the aperture before any trig. * `SvgHeightVisibility` caches the wall activity mask per eye height. * A native result keeps its packed doubles; the cone outline path is built @@ -789,3 +792,58 @@ Instruments: `tool/svg_height_drag_bench_test.dart` and `ICARUS_SVG_NATIVE_LIBRARY` to the built `icarus_height.dll`), `integration_test/view_cone_drag_performance_test.dart` and `view_cone_drag_timeline_test.dart` under `flutter drive --profile -d windows`. + +## How a cone is computed (2026-10-05) + +A cone's outline is a fan of rays from the eye, joined by straight lines. It +is exact when every place the visible wall changes has a ray: each corner the +eye can see (with rays 1e-8 radians either side where the wall turns away), +each crossing of two strokes, and each wall's crossing of the range circle. +Rays aimed at corners nobody can see only add points in the middle of a wall +that is already in the outline. The query's job is to cast the first kind and +skip the second, without ever skipping the first. + +Every ray starts at the same eye, so the query works in angles from it, the +way a 2D renderer does: + +* **Angular bins.** The edges the eye may see are filed once per query into + bins about 2π/2048 radians wide, each bin sorted nearest first. A ray tests + only its own bin's edges and stops at the first one that starts beyond its + hit. Ties go to the lowest edge id. +* **Depths proven by walls.** A run of consecutive edges along a wall ring that + crosses a bin from one boundary to the next, without leaving the bin, is an + unbroken wall across it. Every ray in the bin stops no farther than that + run's farthest point there, so that is the bin's depth: a one-dimensional + depth buffer whose values are proofs rather than samples. +* **Front-to-back culling.** The map's edge tree is walked nearest node first. + A node, an edge or a corner that begins beyond the depth of every bin it + spans is provably hidden and skipped: no filing, no events, no rays. Depths + are re-proven as the walk gets twice as far out (four times, after the first + wave), so nearby walls hide most of the map before it is touched. +* **Margins.** Every proof uses a margin far larger than the rounding in it + (1e-9 in angle and relative distance). Anything the depths cannot rule out + is tested exactly, as before. + +Dart (the web) and `native/height` (desktop) run the same algorithm; native +runs on the calling thread, the old thread pool is gone. + +Checked on 2026-10-05: + +* Against the previous native query on a grid over all 26 map sides, three + apertures and two ranges: 145,872 cones, none whose outline differs by more + than 1e-5 SVG units. The comparison skips the 1e-8 sliver beside each + silhouette, which any ray-built outline draws as a chord. +* `test/svg_cone_exact_test.dart` checks the outline between every pair of + points against an exact ray on two busy maps and five tight spots, and fails + if the hidden-corner test is made even 3% too eager. +* Per cone on Lotus and Breeze: rays fall from about 690 and 1,020 on + average to about 290 and 260, edge tests from about 17,600 and 32,300 to + about 2,000 and 1,400. +* Web, dragging a 103° cone at full length through Lotus and Breeze in Edge: + the query's p99 went from 11–13 ms and 18 ms to 2.5 ms and 4.4 ms, and the + worst query from 15–18 ms and 22 ms to about 3 ms and 6 ms. +* Native, over the same grid: p99 from 3.5 ms to 1.6 ms, p50 from 0.26 ms to + 0.18 ms. + +The slowest cones left are full circles in open areas such as Breeze mid, +where the visible outline itself has some 2,800 corners. diff --git a/lib/view_cone/svg_height_visibility.dart b/lib/view_cone/svg_height_visibility.dart index a2a70b20..027b7e3b 100644 --- a/lib/view_cone/svg_height_visibility.dart +++ b/lib/view_cone/svg_height_visibility.dart @@ -8,6 +8,7 @@ import 'svg_floor_visibility.dart'; /// Offline SVG ink footprints with vertical bands above connected local ground. /// Ground slopes are not geometry here. The caller clips the result to SVG fill. + class SvgHeightVisibility { SvgHeightVisibility._( this.walls, @@ -320,6 +321,11 @@ class SvgHeightVisibility { List? _seamMask; Int8List? _seams; + // Reused by every Dart cone query, so a dragged cone allocates little. + final _bins = _ConeBins(); + Int32List? _seenVertices; + var _seenStamp = 0; + bool _findSeam(int vertex, List active) { final topology = _topology; final live = [ @@ -963,48 +969,27 @@ class SvgHeightVisibility { nativeMicros: result.queryMicros), eyeElevationMeters: ground == null ? null : eye); } - // The native query's algorithm, in Dart for the web: the same events, - // rays and outline as ish_query in native/height. Every structure is a - // list or an index: hashing points and angles dominated the cost once - // Dart is compiled to JavaScript. + // Every ray starts at the eye, so the walls near it are filed by the + // angle they cover (see _ConeBins): a ray tests only the few edges in its + // bin, and a corner that a nearer wall provably hides casts no rays. final half = apertureRadians / 2; - final arcAngles = Float64List(arcSteps + 1); - for (var i = 0; i <= arcSteps; i++) { - arcAngles[i] = -half + apertureRadians * i / arcSteps; - } - final arcHits = [ - for (final angle in arcAngles) - _cast( - origin, - Offset(math.cos(directionRadians + angle), - math.sin(directionRadians + angle)), - range, - active, - stats), + final whole = half + 1e-6 >= math.pi; + final bins = _bins; + bins.file(this, origin, directionRadians, range, active, + lowest: whole ? -math.pi : -half, + span: whole ? 2 * math.pi : apertureRadians); + + final angles = [ + for (var i = 0; i <= arcSteps; i++) -half + apertureRadians * i / arcSteps ]; - bool hiddenEvent(double angle, Offset delta) { - if (apertureRadians / arcSteps >= math.pi || - angle <= -half || - angle >= half) return false; - final interval = ((angle + half) / apertureRadians * arcSteps) - .floor() - .clamp(0, arcSteps - 1); - final first = arcHits[interval]?._edgeIndex; - if (first == null || - first != arcHits[interval + 1]?._edgeIndex || - angle - _cornerOffset < arcAngles[interval] || - angle + _cornerOffset > arcAngles[interval + 1]) return false; - final distance = delta.distance; - final hit = - _edges[first].intersection(origin, delta / distance, distance); - return hit != null && hit < distance - 1e-7; - } - - final angles = [...arcAngles]; // Rays aimed at a vertex or at a wall's crossing of the range circle // stay in the outline even when their neighbours meet the same edge. final vertexAngles = []; void add(double event, {required bool vertex}) { + // Behind a full circle, a ray beside the seam wraps round to the + // other end: -pi and pi are the same direction. + if (whole && event < -half) event += 2 * math.pi; + if (whole && event > half) event -= 2 * math.pi; if (event < -half || event > half) return; angles.add(event); if (vertex) vertexAngles.add(event); @@ -1016,40 +1001,33 @@ class SvgHeightVisibility { if (beside) add(angle + _cornerOffset, vertex: false); } - // A vertex outside the aperture, with slack for the corner offsets, - // cannot start a ray inside it: skip its trigonometry. - final facing = - Offset(math.cos(directionRadians), math.sin(directionRadians)); - final cullSector = half + 1e-6 < math.pi; - final cosSlack = math.cos(math.min(math.pi, half + 1e-6)); - bool inSector(Offset delta) => - !cullSector || _dot(delta, facing) >= cosSlack * delta.distance; - double angleOf(Offset delta) { - final relative = math.atan2(delta.dy, delta.dx) - directionRadians; - return math.atan2(math.sin(relative), math.cos(relative)); + // Whether an event at [angle] could start a ray in the aperture and is + // not provably behind a nearer wall. + bool open(double angle, double distance) { + if (!whole && (angle < -half - 1e-6 || angle > half + 1e-6)) { + return false; + } + return !bins.hides(angle, distance); } final rangeSquared = range * range; // Overlapping painted strokes create visibility corners at their crossing. - // Events behind a proven nearer straight wall cannot change the boundary. for (final crossing in _crossings) { if (!active[crossing.$2] || !active[crossing.$3]) continue; final delta = crossing.$1 - origin; - if (delta.distanceSquared > rangeSquared || delta == Offset.zero) { - continue; - } - if (!inSector(delta)) continue; - final angle = angleOf(delta); - if (hiddenEvent(angle, delta)) continue; + final distanceSquared = delta.distanceSquared; + if (distanceSquared > rangeSquared || delta == Offset.zero) continue; + final angle = bins.angleOf(delta.dx, delta.dy); + if (!open(angle, math.sqrt(distanceSquared))) continue; emit(angle); } final preparationMicros = timer.elapsedMicroseconds; - final candidates = []; - _tree?.query(Rect.fromCircle(center: origin, radius: range), candidates); final topology = _topology; - for (final id in candidates) { + final seen = _seenVertices ??= Int32List(topology.points.length); + final stamp = ++_seenStamp; + for (var slot = 0; slot < bins.count; slot++) { + final id = bins.edgeIds[slot]; final edge = _edges[id]; - if (!active[edge.wall]) continue; // A long wall can cross the range circle without either endpoint being // in range. Seed that exact transition so the polygon follows the wall // all the way to the circle instead of cutting diagonally short of it. @@ -1063,12 +1041,13 @@ class SvgHeightVisibility { final remaining = rangeSquared - closest.distanceSquared; if (remaining >= 0) { final offset = math.sqrt(remaining / lengthSquared); - for (final t in [projection - offset, projection + offset]) { + for (var end = 0; end < 2; end++) { + final t = end == 0 ? projection - offset : projection + offset; if (t < 0 || t > 1) continue; final delta = relative + segment * t; - if (!inSector(delta)) continue; - final angle = angleOf(delta); - if (angle >= -half && angle <= half && !hiddenEvent(angle, delta)) { + final angle = bins.angleOf(delta.dx, delta.dy); + if (!open(angle, range)) continue; + if (angle >= -half && angle <= half) { angles.add(angle); vertexAngles.add(angle); } @@ -1076,36 +1055,18 @@ class SvgHeightVisibility { } for (var end = 0; end < 2; end++) { final vertex = end == 0 ? topology.aVertex[id] : topology.bVertex[id]; + if (seen[vertex] == stamp) continue; + seen[vertex] = stamp; + final distance = bins.vertexDistance[vertex]; + if (distance > range || distance == 0) continue; + final angle = bins.vertexAngle[vertex]; + if (!open(angle, distance)) continue; if (_seam(vertex, active)) continue; - final delta = (end == 0 ? edge.a : edge.b) - origin; - final distanceSquared = delta.distanceSquared; - if (distanceSquared > rangeSquared || delta == Offset.zero) continue; - if (!inSector(delta)) continue; - final angle = angleOf(delta); - // If both bounding arc rays hit the same straight segment, that - // segment covers the angular interval. A vertex strictly behind it - // cannot change the visible boundary. Nearer corners still add rays. - if (apertureRadians / arcSteps < math.pi && - angle > -half && - angle < half) { - final interval = ((angle + half) / apertureRadians * arcSteps) - .floor() - .clamp(0, arcSteps - 1); - final first = arcHits[interval]?._edgeIndex; - final last = arcHits[interval + 1]?._edgeIndex; - if (first != null && - first == last && - angle - _cornerOffset >= arcAngles[interval] && - angle + _cornerOffset <= arcAngles[interval + 1]) { - final distance = math.sqrt(distanceSquared); - final hit = - _edges[first].intersection(origin, delta / distance, distance); - if (hit != null && hit < distance - 1e-7) continue; - } - } // Rays just beside a vertex the wall runs straight across meet its // two edges, so only the vertex ray adds a corner. - emit(angle, beside: !_passThrough(vertex, delta, active)); + emit(angle, + beside: !_passThrough( + vertex, topology.points[vertex] - origin, active)); } } final sorted = _sortedUnique(angles); @@ -1119,10 +1080,7 @@ class SvgHeightVisibility { for (final angle in sorted) { final direction = Offset(math.cos(directionRadians + angle), math.sin(directionRadians + angle)); - final arc = _indexOf(arcAngles, angle); - final hit = arc >= 0 - ? arcHits[arc] - : _cast(origin, direction, range, active, stats); + final hit = bins.cast(this, origin, direction, angle, range, stats); final point = hit?.point ?? origin + direction * range; // Rays meeting the same edge in a row lie on one straight line; the // middle ones add nothing to the outline. Vertex rays always stay. @@ -1840,3 +1798,612 @@ bool _betweenOnSameLine(Offset first, Offset middle, Offset last) { final roundoff = 64 * 2.220446049250313e-16 * math.max(1.0, length); return _cross(offset, span).abs() <= roundoff * length; } + +/// The active edges a cone's eye may see, filed by the angles they cover. +/// +/// Every ray of a cone leaves the same eye. So rather than walk the map's +/// edge tree once per ray, a query files the edges it may meet once: each +/// gets its nearest distance to the eye and the span of angles it covers, and +/// goes into every bin of that span. A ray then tests only its bin's edges. +/// +/// Each bin also gets a depth, the farthest any ray in it can travel. A run +/// of consecutive ring edges that crosses a bin from one boundary to the next +/// is an unbroken wall across it, so every ray in the bin stops no farther +/// than the run's farthest point there. The nearest such bound is the bin's +/// depth. Whatever begins beyond the depths it spans is hidden: an edge there +/// is not filed, and a corner there needs no rays, since they would land in +/// the middle of whatever hides it. +/// +/// The edge tree is walked nearest node first, the way a renderer draws front +/// to back, and a node whose bounds begin beyond the depths of every bin they +/// span is skipped whole: most of a map is behind the walls nearest the eye. +/// Depths are only ever used with a margin, and whatever they let through is +/// tested exactly, so the outline is the one every edge would give. +class _ConeBins { + /// Filed edges, by slot: edge id, nearest distance to the eye, and the span + /// of angles from the cone's direction they cover, which may run past pi. + int count = 0; + Int32List edgeIds = Int32List(64); + Float64List near = Float64List(64); + Float64List from = Float64List(64); + Float64List to = Float64List(64); + + /// Bins split [lowest, lowest + binCount * width] evenly. + int binCount = 0; + double lowest = 0, width = 1; + + /// Bin k's filed slots are items[offsets[k]] up to items[offsets[k + 1]]. + Int32List offsets = Int32List(1); + Int32List items = Int32List(0); + + /// The farthest any ray in a bin can travel; infinite when unproven. + Float64List depth = Float64List(0); + + /// The largest depth over ranges of bins: a segment tree whose leaves, + /// from [_leaves] on, are the depths. + Float64List _peaks = Float64List(0); + int _leaves = 0; + + /// The world direction of each bin boundary. + Float64List boundaryX = Float64List(0), boundaryY = Float64List(0); + + double _cosine = 1, _sine = 0; + var _whole = false; + double _relativeLowest = double.nan, _relativeWidth = double.nan; + Float64List _relativeX = Float64List(0), _relativeY = Float64List(0); + + /// Per vertex this query, where [vertexMark] is [_vertexStamp]: its angle + /// from the cone's direction and its distance from the eye. + Float64List vertexAngle = Float64List(0), vertexDistance = Float64List(0); + Int32List vertexMark = Int32List(0); + var _vertexStamp = 0; + double _ox = 0, _oy = 0; + late _VertexTopology _topology; + + void _see(int vertex) { + if (vertexMark[vertex] == _vertexStamp) return; + vertexMark[vertex] = _vertexStamp; + final point = _topology.points[vertex]; + final dx = point.dx - _ox, dy = point.dy - _oy; + vertexAngle[vertex] = angleOf(dx, dy); + vertexDistance[vertex] = math.sqrt(dx * dx + dy * dy); + } + + /// Margins on the proofs: rounding in the angles and distances here is + /// some 1e-15, far inside these. + static const _angleMargin = 1e-9; + static const _relativeMargin = 1e-9, _absoluteMargin = 1e-9; + + /// Rays beside a corner leave this far from its angle. + static const _cornerSlack = SvgHeightVisibility._cornerOffset + _angleMargin; + + /// The narrowest bin; a cone's bins are about this wide. + static const _binWidth = 2 * math.pi / 2048; + + /// The angle of [dx], [dy] from the cone's direction, in [-pi, pi]. + double angleOf(double dx, double dy) => + math.atan2(-dx * _sine + dy * _cosine, dx * _cosine + dy * _sine); + + /// [index] held to the boundaries just outside the cone; NaN, from an + /// aperture so small its bins have no width, to the first. + double _boundaryIndex(double index) => + index >= -1 ? (index <= binCount + 1 ? index : binCount + 1.0) : -1.0; + + int binOf(double angle) { + // Clamped as a double: a narrow cone's bins are so fine that angles + // outside it can lie past any integer. + final k = ((angle - lowest) / width).floorToDouble(); + if (!(k >= 0)) return 0; + return k >= binCount ? binCount - 1 : k.toInt(); + } + + /// Whether every ray within a corner offset of [angle] provably stops + /// before [distance]. + bool hides(double angle, double distance) { + final last = binOf(angle + _cornerSlack); + for (var k = binOf(angle - _cornerSlack); k <= last; k++) { + if (!_beyond(distance, depth[k])) return false; + } + // Behind a full circle, rays beside an angle at the seam wrap round. + if (_whole) { + if (angle - _cornerSlack < lowest && + !_beyond(distance, depth[binCount - 1])) { + return false; + } + if (angle + _cornerSlack > lowest + binCount * width && + !_beyond(distance, depth[0])) { + return false; + } + } + return true; + } + + static bool _beyond(double distance, double depth) => + distance > depth * (1 + _relativeMargin) + _absoluteMargin; + + void file(SvgHeightVisibility model, Offset origin, double direction, + double range, List active, + {required double lowest, required double span}) { + _cosine = math.cos(direction); + _sine = math.sin(direction); + this.lowest = lowest; + _topology = model._topology; + _ox = origin.dx; + _oy = origin.dy; + final vertices = _topology.points.length; + if (vertexMark.length < vertices) { + vertexAngle = Float64List(vertices); + vertexDistance = Float64List(vertices); + vertexMark = Int32List(vertices); + } + _vertexStamp++; + final edgeCount = model._edges.length; + if (_slotMark.length < edgeCount) { + _slotMark = Int32List(edgeCount); + _dirty = Int32List(edgeCount); + } + _mark++; + final bins = math.max(64, (span / _binWidth).ceil()); + binCount = bins; + width = span / bins; + _whole = span >= 2 * math.pi; + if (depth.length < bins) { + depth = Float64List(bins); + boundaryX = Float64List(bins + 1); + boundaryY = Float64List(bins + 1); + offsets = Int32List(bins + 1); + } + // Boundary directions relative to the cone's direction depend only on + // the bins, so they are kept between queries and turned to face it. + if (_relativeLowest != lowest || + _relativeWidth != width || + _relativeX.length < bins + 1) { + _relativeLowest = lowest; + _relativeWidth = width; + _relativeX = Float64List(bins + 1); + _relativeY = Float64List(bins + 1); + for (var j = 0; j <= bins; j++) { + _relativeX[j] = math.cos(lowest + j * width); + _relativeY[j] = math.sin(lowest + j * width); + } + } + for (var j = 0; j <= bins; j++) { + final x = _relativeX[j], y = _relativeY[j]; + boundaryX[j] = x * _cosine - y * _sine; + boundaryY[j] = x * _sine + y * _cosine; + } + depth.fillRange(0, bins, double.infinity); + _leaves = 1; + while (_leaves < bins) { + _leaves *= 2; + } + if (_peaks.length < 2 * _leaves) _peaks = Float64List(2 * _leaves); + _peaks.fillRange(0, 2 * _leaves, double.infinity); + count = 0; + + final root = model._tree; + if (root != null) { + _walk(model, root, origin, range, active); + } + + // File each edge in the bins it may be seen in: count, then place. + offsets.fillRange(0, bins + 1, 0); + _place(false); + for (var k = 0; k < bins; k++) { + offsets[k + 1] += offsets[k]; + } + if (items.length < offsets[bins]) { + items = Int32List(math.max(offsets[bins], items.length * 2)); + } + _place(true); + // Placing advanced each bin's offset to the next bin's start. + for (var k = bins; k > 0; k--) { + offsets[k] = offsets[k - 1]; + } + offsets[0] = 0; + // Nearest first within each bin, so a ray stops at the first edge that + // begins beyond its hit. The walk filed them nearly in this order. + for (var k = 0; k < bins; k++) { + final end = offsets[k + 1]; + for (var i = offsets[k] + 1; i < end; i++) { + final slot = items[i]; + final key = near[slot]; + var j = i - 1; + while (j >= offsets[k] && near[items[j]] > key) { + items[j + 1] = items[j]; + j--; + } + items[j + 1] = slot; + } + } + } + + // Tree nodes waiting to be visited, nearest first: a binary heap. + final _heapNodes = <_EdgeNode>[]; + Float64List _heapKeys = Float64List(64); + + void _push(_EdgeNode node, double key) { + var i = _heapNodes.length; + _heapNodes.add(node); + if (_heapKeys.length <= i) { + _heapKeys = Float64List(_heapKeys.length * 2)..setAll(0, _heapKeys); + } + while (i > 0) { + final parent = (i - 1) >> 1; + if (_heapKeys[parent] <= key) break; + _heapNodes[i] = _heapNodes[parent]; + _heapKeys[i] = _heapKeys[parent]; + i = parent; + } + _heapNodes[i] = node; + _heapKeys[i] = key; + } + + /// Removes the nearest node; its key is left in [_popped]. + _EdgeNode _pop() { + final top = _heapNodes[0]; + _popped = _heapKeys[0]; + final node = _heapNodes.removeLast(); + final n = _heapNodes.length; + if (n > 0) { + final key = _heapKeys[n]; + var i = 0; + while (true) { + var child = 2 * i + 1; + if (child >= n) break; + if (child + 1 < n && _heapKeys[child + 1] < _heapKeys[child]) child++; + if (_heapKeys[child] >= key) break; + _heapNodes[i] = _heapNodes[child]; + _heapKeys[i] = _heapKeys[child]; + i = child; + } + _heapNodes[i] = node; + _heapKeys[i] = key; + } + return top; + } + + double _popped = 0; + + static double _distanceTo(Rect bounds, double x, double y) { + final dx = math.max(0.0, math.max(bounds.left - x, x - bounds.right)); + final dy = math.max(0.0, math.max(bounds.top - y, y - bounds.bottom)); + return math.sqrt(dx * dx + dy * dy); + } + + void _walk(SvgHeightVisibility model, _EdgeNode root, Offset origin, + double range, List active) { + final ox = origin.dx, oy = origin.dy; + _heapNodes.clear(); + _push(root, _distanceTo(root.bounds, ox, oy)); + // Depths are proven from the edges filed so far, again each time the + // walk passes twice as far from the eye. + var wave = math.max(range / 16, 1e-3); + var proven = 0; + while (_heapNodes.isNotEmpty) { + final node = _pop(); + final distance = _popped; + if (distance > range) break; + if (distance > wave) { + if (count > proven) { + _walkChains(model, origin, _whole, proven); + _raisePeaks(); + proven = count; + } + while (wave < distance) { + wave *= 4; + } + } + if (!_boundsMayShow(node.bounds, distance, ox, oy)) continue; + final ids = node.ids; + if (ids != null) { + for (final id in ids) { + _fileEdge(model, id, ox, oy, range, active); + } + } else { + final left = node.left!, right = node.right!; + final leftKey = _distanceTo(left.bounds, ox, oy); + if (leftKey <= range) _push(left, leftKey); + final rightKey = _distanceTo(right.bounds, ox, oy); + if (rightKey <= range) _push(right, rightKey); + } + } + if (count > proven) _walkChains(model, origin, _whole, proven); + } + + /// Whether anything in [bounds], [distance] from the eye at its nearest, + /// could show in the cone. + bool _boundsMayShow(Rect bounds, double distance, double ox, double oy) { + final column = ox < bounds.left ? 0 : (ox > bounds.right ? 2 : 1); + final row = oy < bounds.top ? 0 : (oy > bounds.bottom ? 2 : 1); + if (distance < _absoluteMargin || (column == 1 && row == 1)) return true; + // Seen from outside, a box spans less than half a turn, between the two + // corners that bound its silhouette. Which two depends only on where the + // eye is around the box: per region, each corner as (right?, bottom?). + final corners = _silhouettes[row * 3 + column]; + final ax = ((corners & 8) != 0 ? bounds.right : bounds.left) - ox; + final ay = ((corners & 4) != 0 ? bounds.bottom : bounds.top) - oy; + final bx = ((corners & 2) != 0 ? bounds.right : bounds.left) - ox; + final by = ((corners & 1) != 0 ? bounds.bottom : bounds.top) - oy; + final ta = angleOf(ax, ay), tb = angleOf(bx, by); + final turn = ax * by - ay * bx; + final first = turn > 0 ? ta : (turn < 0 ? tb : math.min(ta, tb)); + var last = turn > 0 ? tb : (turn < 0 ? ta : math.max(ta, tb)); + if (last < first) last += 2 * math.pi; + return _mayShow(first - _cornerSlack, last + _cornerSlack, distance); + } + + /// For each region around a box, rows top to bottom and columns left to + /// right, the two corners of its silhouette as bits: first corner right, + /// first bottom, second right, second bottom. + static const _silhouettes = [ + 0x9, 0x2, 0x3, // + 0x1, 0x0, 0xb, // + 0x3, 0x7, 0x9, + ]; + + /// Whether something [distance] from the eye, spanning angles [first] to + /// [last] (which may run past pi), falls in the cone and is not provably + /// behind the depths there. + bool _mayShow(double first, double last, double distance) { + const margin = 1e-6; + final highest = lowest + binCount * width; + for (var shift = 2 * math.pi; shift > -4 * math.pi; shift -= 2 * math.pi) { + final start = first + shift, end = last + shift; + if (end < lowest - margin || start > highest + margin) continue; + if (!_beyond(distance, _peak(binOf(start), binOf(end)))) return true; + } + return false; + } + + /// The largest depth over bins [first] to [last]. + double _peak(int first, int last) { + var result = 0.0; + var low = first + _leaves, high = last + _leaves + 1; + while (low < high) { + if (low.isOdd) result = math.max(result, _peaks[low++]); + if (high.isOdd) result = math.max(result, _peaks[--high]); + low >>= 1; + high >>= 1; + } + return result; + } + + void _raisePeaks() { + _peaks.setRange(_leaves, _leaves + binCount, depth); + for (var i = _leaves - 1; i > 0; i--) { + _peaks[i] = math.max(_peaks[2 * i], _peaks[2 * i + 1]); + } + } + + void _fileEdge(SvgHeightVisibility model, int id, double ox, double oy, + double range, List active) { + final edge = model._edges[id]; + if (!active[edge.wall]) return; + final ax = edge.a.dx - ox, ay = edge.a.dy - oy; + final bx = edge.b.dx - ox, by = edge.b.dy - oy; + final ex = edge._ex, ey = edge._ey; + final lengthSquared = ex * ex + ey * ey; + var t = lengthSquared == 0 ? 0.0 : -(ax * ex + ay * ey) / lengthSquared; + t = t < 0 ? 0.0 : (t > 1 ? 1.0 : t); + final cx = ax + ex * t, cy = ay + ey * t; + final distance = math.sqrt(cx * cx + cy * cy); + if (distance > range) return; + final va = _topology.aVertex[id], vb = _topology.bVertex[id]; + _see(va); + _see(vb); + var first = -math.pi, last = math.pi; + if (distance >= _absoluteMargin) { + final ta = vertexAngle[va], tb = vertexAngle[vb]; + final turn = ax * by - ay * bx; + if (turn > 0) { + first = ta; + last = tb; + } else if (turn < 0) { + first = tb; + last = ta; + } else { + first = math.min(ta, tb); + last = math.max(ta, tb); + } + if (last < first) last += 2 * math.pi; + if (!_mayShow(first - _cornerSlack, last + _cornerSlack, distance)) { + return; + } + } + if (count == edgeIds.length) { + final size = count * 2; + edgeIds = Int32List(size)..setAll(0, edgeIds); + near = Float64List(size)..setAll(0, near); + from = Float64List(size)..setAll(0, from); + to = Float64List(size)..setAll(0, to); + } + if (distance >= _absoluteMargin) _slotMark[id] = _mark; + edgeIds[count] = id; + near[count] = distance; + from[count] = first; + to[count] = last; + count++; + } + + /// Edges filed this query that a run can pass through: [_slotMark] is + /// [_mark]. An edge through the eye stops nothing reliably beside it. + Int32List _slotMark = Int32List(0); + var _mark = 0; + + // The runs to walk this wave: [_dirty] is [_dirtyStamp] on their edges. + Int32List _dirty = Int32List(0); + var _dirtyStamp = 0; + final _heads = []; + + /// Lowers each bin's depth to the farthest point of any run of consecutive + /// ring edges that crosses the bin from one boundary to the next. Such a run + /// is an unbroken wall across the bin, so it stops every ray in it. Single + /// edges rarely span a bin: Riot's outlines are many short strokes. + void _walkChains( + SvgHeightVisibility model, Offset origin, bool whole, int firstNew) { + final edges = model._edges; + final mark = _mark; + bool chained(int id) => + id >= 0 && id < edges.length && _slotMark[id] == mark; + final ox = origin.dx, oy = origin.dy; + final aVertex = _topology.aVertex, bVertex = _topology.bVertex; + // Runs already walked have proven all they can; walk again only those + // with an edge filed since, each from its first edge. + final stamp = ++_dirtyStamp; + _heads.clear(); + for (var slot = firstNew; slot < count; slot++) { + var id = edgeIds[slot]; + if (_slotMark[id] != mark) continue; + while (_dirty[id] != stamp) { + _dirty[id] = stamp; + if (!chained(id - 1) || bVertex[id - 1] != aVertex[id]) { + _heads.add(id); + break; + } + id--; + } + } + for (final head in _heads) { + var id = head; + var raw = vertexAngle[aVertex[id]]; + // The run's angle, unwrapped so it never jumps by a turn. + var angle = raw; + // The last boundary crossed, and since then: the farthest point and + // the angles the run has reached. + var hasLast = false; + var last = 0; + var farthest = 0.0, lowAngle = 0.0, highAngle = 0.0; + while (true) { + final edge = edges[id]; + final next = bVertex[id]; + final nextRaw = vertexAngle[next]; + var turn = nextRaw - raw; + if (turn > math.pi) { + turn -= 2 * math.pi; + } else if (turn <= -math.pi) { + turn += 2 * math.pi; + } + final end = angle + turn; + final ex = edge._ex, ey = edge._ey; + final along = (edge.a.dx - ox) * ey - (edge.a.dy - oy) * ex; + // Boundaries the edge plainly crosses, in the order it crosses them. + final step = end > angle ? 1 : -1; + var start = step > 0 + ? ((angle + _angleMargin - lowest) / width).ceilToDouble() + : ((angle - _angleMargin - lowest) / width).floorToDouble(); + var finish = step > 0 + ? ((end - _angleMargin - lowest) / width).floorToDouble() + : ((end + _angleMargin - lowest) / width).ceilToDouble(); + if (!whole) { + // Boundaries outside the cone are skipped below, and a narrow + // cone's bins are so fine those can lie past any integer. + start = _boundaryIndex(start); + finish = _boundaryIndex(finish); + } + var j = start.toInt(); + final stop = finish.toInt(); + for (; step > 0 ? j <= stop : j >= stop; j += step) { + if (!whole && (j < 0 || j > binCount)) continue; + final boundary = whole ? j % binCount : j; + final distance = + along / (boundaryX[boundary] * ey - boundaryY[boundary] * ex); + final at = lowest + j * width; + if (!(distance >= 0) || distance.isInfinite) { + hasLast = false; + continue; + } + if (hasLast && (j - last).abs() == 1) { + final previous = lowest + last * width; + // Between the two crossings the run stayed inside the bin. + if (lowAngle >= math.min(at, previous) - 1e-8 && + highAngle <= math.max(at, previous) + 1e-8) { + final k = + whole ? math.min(j, last) % binCount : math.min(j, last); + final bound = math.max(farthest, distance); + if (bound < depth[k]) depth[k] = bound; + } + } + hasLast = true; + last = j; + farthest = distance; + lowAngle = highAngle = at; + } + farthest = math.max(farthest, vertexDistance[next]); + lowAngle = math.min(lowAngle, end); + highAngle = math.max(highAngle, end); + if (!chained(id + 1) || aVertex[id + 1] != next) break; + id++; + raw = nextRaw; + angle = end; + } + } + } + + /// Counts each slot into offsets[k + 1] for every bin k it is filed in, or, + /// with [fill], writes it at offsets[k] and advances that. + void _place(bool fill) { + final highest = lowest + binCount * width; + for (var slot = 0; slot < count; slot++) { + final distance = near[slot]; + if (distance < _absoluteMargin) { + // An edge through the eye can stop a ray in any direction. + for (var k = 0; k < binCount; k++) { + if (fill) { + items[offsets[k]++] = slot; + } else { + offsets[k + 1]++; + } + } + continue; + } + // Rays are admitted a sliver past an edge's ends; see _Edge.intersection. + final pad = _angleMargin + 1e-12 / distance; + // A span starting just below -pi also covers the bins just below pi. + for (var shift = 2 * math.pi; + shift > -4 * math.pi; + shift -= 2 * math.pi) { + final start = from[slot] + shift - pad, end = to[slot] + shift + pad; + if (end < lowest) break; + if (start > highest) continue; + final first = binOf(start), last = binOf(end); + for (var k = first; k <= last; k++) { + if (distance > depth[k] * (1 + _relativeMargin) + _absoluteMargin) { + continue; + } + if (fill) { + items[offsets[k]++] = slot; + } else { + offsets[k + 1]++; + } + } + } + } + } + + /// The nearest hit along [direction], [angle] from the cone's direction, + /// within [range]; ties go to the lowest edge id. + SvgVisibilityHit? cast(SvgHeightVisibility model, Offset origin, + Offset direction, double angle, double range, _Counters stats) { + final k = binOf(angle); + stats.cells++; + var best = range; + var bestId = -1; + for (var i = offsets[k]; i < offsets[k + 1]; i++) { + final slot = items[i]; + if (near[slot] > best) break; + final id = edgeIds[slot]; + stats.edgeTests++; + final distance = model._edges[id].intersection(origin, direction, best); + if (distance == null) continue; + if (bestId < 0 || distance < best || (distance == best && id < bestId)) { + best = distance; + bestId = id; + } + } + if (bestId < 0) return null; + final wall = model._edges[bestId].wall; + if (model.walls[wall].unknownHeight) stats.unknownHits++; + return model._hit(origin, direction, best, wall, bestId); + } +} diff --git a/native/height/icarus_svg_height.cpp b/native/height/icarus_svg_height.cpp index b62bee52..7c811771 100644 --- a/native/height/icarus_svg_height.cpp +++ b/native/height/icarus_svg_height.cpp @@ -7,25 +7,15 @@ #include #include #include +#include #include #include -#include -#include -#include -#include -#include #include #include +#include #include #include -#ifdef _WIN32 -#ifndef NOMINMAX -#define NOMINMAX -#endif -#include -#endif - static_assert(sizeof(ISHResult) == 72); static_assert(offsetof(ISHResult, points) == 8); static_assert(offsetof(ISHResult, queryMicros) == 64); @@ -33,8 +23,8 @@ static_assert(offsetof(ISHResult, queryMicros) == 64); namespace { using Clock = std::chrono::steady_clock; constexpr double pi = 3.141592653589793238462643383279502884; +constexpr double infinity = std::numeric_limits::infinity(); constexpr double cornerOffset = 1e-8; -constexpr double slabPadding = 1e-10; constexpr uint32_t maximumEdges = 1u << 19; constexpr uint32_t maximumWalls = 1u << 20; constexpr uint32_t maximumArcSteps = 4096; @@ -153,44 +143,16 @@ bool overlaps(const Bounds &a, const Bounds &b) { a.bottom < b.top); } -struct PreparedRay { - Point origin, direction; - double inverseX, inverseY; -}; - -bool entry(const Bounds &bounds, const PreparedRay &ray, double range, - double &result) { - double lo = 0.0, hi = range; - if (ray.direction.x == 0) { - if (ray.origin.x < bounds.left - slabPadding || - ray.origin.x > bounds.right + slabPadding) - return false; - } else { - const double a = - (bounds.left - slabPadding - ray.origin.x) * ray.inverseX; - const double b = - (bounds.right + slabPadding - ray.origin.x) * ray.inverseX; - lo = std::max(lo, std::min(a, b)); - hi = std::min(hi, std::max(a, b)); - if (lo > hi) - return false; - } - if (ray.direction.y == 0) { - if (ray.origin.y < bounds.top - slabPadding || - ray.origin.y > bounds.bottom + slabPadding) - return false; +void collectCandidates(const Node *node, const Bounds &area, + std::vector &output) { + if (!overlaps(node->bounds, area)) + return; + if (!node->ids.empty()) { + output.insert(output.end(), node->ids.begin(), node->ids.end()); } else { - const double a = - (bounds.top - slabPadding - ray.origin.y) * ray.inverseY; - const double b = - (bounds.bottom + slabPadding - ray.origin.y) * ray.inverseY; - lo = std::max(lo, std::min(a, b)); - hi = std::min(hi, std::max(a, b)); - if (lo > hi) - return false; + collectCandidates(node->left.get(), area, output); + collectCandidates(node->right.get(), area, output); } - result = lo; - return true; } struct Hit { @@ -199,144 +161,615 @@ struct Hit { uint32_t edge = 0; }; -struct Counters { - uint64_t edgeTests = 0, nodes = 0; -}; +struct Crossing { Point point; uint32_t first, second; }; -// One per chunk, padded to a cache line: neighbouring chunks run on different -// threads and must not bounce the same line while counting. -struct alignas(64) ChunkCounters { - Counters value; - char padding[64 - sizeof(Counters)]; -}; +// Moves a per-query stamp on, clearing its marks the one time in 2^32 it +// wraps, so a stale mark never matches. +void advance(uint32_t &stamp, std::vector &marks) { + if (++stamp == 0) { + std::fill(marks.begin(), marks.end(), 0u); + stamp = 1; + } +} -void collectCandidates(const Node *node, const Bounds &area, - std::vector &output); -struct Crossing { Point point; uint32_t first, second; }; +// Margins on the proofs: rounding in the angles and distances here is some +// 1e-15, far inside these. +constexpr double angleMargin = 1e-9; +constexpr double relativeMargin = 1e-9, absoluteMargin = 1e-9; +// Rays beside a corner leave this far from its angle. +constexpr double cornerSlack = cornerOffset + angleMargin; +// The narrowest bin; a cone's bins are about this wide. +constexpr double binWidth = 2 * pi / 2048; -// A persistent pool for the per-query ray casts. Rays are independent and -// the polygon is assembled serially afterwards, so the result is bitwise the -// same as the single-threaded query; only the wall-clock time changes. The -// caller thread works too, so a query never waits on a sleeping worker. +bool beyond(double distance, double depth) { + return distance > depth * (1 + relativeMargin) + absoluteMargin; +} + +double distanceTo(const Bounds &bounds, double x, double y) { + const double dx = std::max(0.0, std::max(bounds.left - x, x - bounds.right)); + const double dy = std::max(0.0, std::max(bounds.top - y, y - bounds.bottom)); + return std::sqrt(dx * dx + dy * dy); +} + +// The active edges a cone's eye may see, filed by the angles they cover. This +// is the Dart _ConeBins in lib/view_cone/svg_height_visibility.dart, which +// carries the proofs. // -// A run is published as one 64-bit ticket: the chunk count in the high half -// and the next chunk index in the low half. A worker claims a chunk with a -// single fetch-add on that word, so the index it receives is always paired -// with the count of the same run. A stale claim from an earlier run carries -// that run's count, fails the bounds test, and touches nothing; a claim -// within range keeps the run alive until the chunk is done, so the callback -// and the remaining counter it then reads belong to that run. -// A query is a few milliseconds of work split across threads, and the -// caller waits for every chunk. A thread the scheduler sets aside for a -// time slice (15 ms on Windows) holding one chunk stalls the whole query, so -// the threads doing a query run above normal priority while they do it. -struct Boost { -#ifdef _WIN32 - Boost() : thread(GetCurrentThread()), previous(GetThreadPriority(thread)) { - if (previous < THREAD_PRIORITY_ABOVE_NORMAL) - SetThreadPriority(thread, THREAD_PRIORITY_ABOVE_NORMAL); - } - ~Boost() { - if (previous < THREAD_PRIORITY_ABOVE_NORMAL) - SetThreadPriority(thread, previous); - } - HANDLE thread; - int previous; -#endif -}; +// Every ray of a cone leaves the same eye. So rather than walk the edge tree +// once per ray, a query files the edges it may meet once: each gets its +// nearest distance to the eye and the span of angles it covers, and goes into +// every bin of that span. A ray then tests only its bin's edges. +// +// Each bin also gets a depth, the farthest any ray in it can travel. A run of +// consecutive ring edges that crosses a bin from one boundary to the next is +// an unbroken wall across it, so every ray in the bin stops no farther than +// the run's farthest point there. Whatever begins beyond the depths it spans +// is hidden: an edge there is not filed, and a corner there needs no rays. +// +// The edge tree is walked nearest node first, and a node whose bounds begin +// beyond the depths of every bin they span is skipped whole. Depths are only +// ever used with a margin, and whatever they let through is tested exactly, +// so the outline is the one every edge would give. +struct ConeBins { + const std::vector *edges = nullptr; + const std::vector *points = nullptr; + const uint8_t *active = nullptr; + + // Filed edges, by slot: edge id, nearest distance to the eye, and the span + // of angles from the cone's direction they cover, which may run past pi. + std::vector edgeIds; + std::vector near, from, to; + + // Bins split [lowest, lowest + binCount * width] evenly. Bin k's filed + // slots are items[offsets[k]] up to items[offsets[k + 1]], nearest first. + uint32_t binCount = 0; + double lowest = 0, width = 1; + std::vector offsets, items; + + // The farthest any ray in a bin can travel; infinite when unproven. + std::vector depth; + // The largest depth over ranges of bins: a segment tree whose leaves, from + // index leaves on, are the depths. + std::vector peaks; + size_t leaves = 0; + // The world direction of each bin boundary. + std::vector boundaryX, boundaryY; + + double cosine = 1, sine = 0, ox = 0, oy = 0; + bool whole = false; + + // Per vertex this query, where vertexMark is vertexStamp: its angle from + // the cone's direction and its distance from the eye. + std::vector vertexAngle, vertexDistance; + std::vector vertexMark; + uint32_t vertexStamp = 0; + + // Edges filed this query that a run can pass through: slotMark is mark. + // An edge through the eye stops nothing reliably beside it. + std::vector slotMark; + uint32_t mark = 0; + // The runs to walk this wave: dirty is dirtyStamp on their edges. + std::vector dirty; + uint32_t dirtyStamp = 0; + std::vector heads; + + // Tree nodes waiting to be visited, nearest first: a binary heap. + std::vector heapNodes; + std::vector heapKeys; + + uint64_t nodes = 0; + + void reserve(size_t vertexCount, size_t edgeCount) { + vertexAngle.resize(vertexCount); + vertexDistance.resize(vertexCount); + vertexMark.assign(vertexCount, 0); + slotMark.assign(edgeCount, 0); + dirty.assign(edgeCount, 0); + } + + size_t count() const { return edgeIds.size(); } + + // The angle of dx, dy from the cone's direction, in [-pi, pi]. + double angleOf(double dx, double dy) const { + return std::atan2(-dx * sine + dy * cosine, dx * cosine + dy * sine); + } + + uint32_t binOf(double angle) const { + const double k = std::floor((angle - lowest) / width); + // NaN, from an aperture so small its bins have no width, goes first. + if (!(k >= 0)) return 0; + return k >= binCount ? binCount - 1 : uint32_t(k); + } -struct Pool { - explicit Pool(unsigned workers) { - for (unsigned i = 0; i < workers; ++i) - threads.emplace_back([this] { -#ifdef _WIN32 - SetThreadPriority(GetCurrentThread(), THREAD_PRIORITY_ABOVE_NORMAL); -#endif - loop(); - }); - } - ~Pool() { - { - std::lock_guard lock(mutex); - stop.store(true, std::memory_order_release); + // index held to the boundaries just outside the cone; NaN, from an + // aperture so small its bins have no width, to the first. + static double boundaryIndex(double index, double bins) { + return index >= -1 ? (index <= bins + 1 ? index : bins + 1) : -1; + } + + // Whether every ray within a corner offset of angle provably stops before + // distance. + bool hides(double angle, double distance) const { + const uint32_t last = binOf(angle + cornerSlack); + for (uint32_t k = binOf(angle - cornerSlack); k <= last; ++k) + if (!beyond(distance, depth[k])) + return false; + // Behind a full circle, rays beside an angle at the seam wrap round. + if (whole) { + if (angle - cornerSlack < lowest && + !beyond(distance, depth[binCount - 1])) + return false; + if (angle + cornerSlack > lowest + binCount * width && + !beyond(distance, depth[0])) + return false; } - wake.notify_all(); - for (std::thread &thread : threads) thread.join(); + return true; } - void run(size_t count, const std::function &task) { - if (count == 0) return; - if (threads.empty() || count == 1 || count > kMaximumChunks) { - for (size_t i = 0; i < count; ++i) task(i); - return; + + void file(const std::vector &edgeList, + const std::vector &vertexPoints, const Node *root, + Point origin, double direction, double range, + const uint8_t *activeWalls, double lowestAngle, double span) { + edges = &edgeList; + points = &vertexPoints; + active = activeWalls; + cosine = std::cos(direction); + sine = std::sin(direction); + lowest = lowestAngle; + ox = origin.x; + oy = origin.y; + advance(vertexStamp, vertexMark); + advance(mark, slotMark); + const uint32_t bins = + std::max(64u, uint32_t(std::ceil(span / binWidth))); + binCount = bins; + width = span / bins; + whole = span >= 2 * pi; + boundaryX.resize(bins + 1); + boundaryY.resize(bins + 1); + for (uint32_t j = 0; j <= bins; ++j) { + const double world = direction + lowest + j * width; + boundaryX[j] = std::cos(world); + boundaryY[j] = std::sin(world); } - { - std::lock_guard lock(mutex); - job = &task; - remaining.store(count, std::memory_order_relaxed); - ticket.store(uint64_t(count) << 32, std::memory_order_release); + depth.assign(bins, infinity); + leaves = 1; + while (leaves < bins) + leaves *= 2; + peaks.assign(2 * leaves, infinity); + edgeIds.clear(); + near.clear(); + from.clear(); + to.clear(); + nodes = 0; + if (root) + walk(root, range); + + // File each edge in the bins it may be seen in: count, then place. + offsets.assign(bins + 1, 0); + place(false); + for (uint32_t k = 0; k < bins; ++k) + offsets[k + 1] += offsets[k]; + if (items.size() < offsets[bins]) + items.resize(std::max(offsets[bins], items.size() * 2)); + place(true); + // Placing advanced each bin's offset to the next bin's start. + for (uint32_t k = bins; k > 0; --k) + offsets[k] = offsets[k - 1]; + offsets[0] = 0; + // Nearest first within each bin, so a ray stops at the first edge that + // begins beyond its hit. The walk filed them nearly in this order. + for (uint32_t k = 0; k < bins; ++k) { + const uint32_t begin = offsets[k], end = offsets[k + 1]; + for (uint32_t i = begin + 1; i < end; ++i) { + const uint32_t slot = items[i]; + const double key = near[slot]; + uint32_t j = i; + while (j > begin && near[items[j - 1]] > key) { + items[j] = items[j - 1]; + --j; + } + items[j] = slot; + } + } + } + + // The nearest hit along direction, angle from the cone's direction, within + // range; ties go to the lowest edge id. + Hit cast(Point origin, Point direction, double angle, double range, + uint64_t &edgeTests) const { + const uint32_t k = binOf(angle); + double best = range; + bool found = false; + uint32_t bestId = 0; + for (uint32_t i = offsets[k]; i < offsets[k + 1]; ++i) { + const uint32_t slot = items[i]; + if (near[slot] > best) + break; + const uint32_t id = edgeIds[slot]; + ++edgeTests; + double distance; + if (!(*edges)[id].intersection(origin, direction, best, distance)) + continue; + if (!found || distance < best || (distance == best && id < bestId)) { + best = distance; + bestId = id; + found = true; + } } - wake.notify_all(); - work(); - // The caller spins on the last chunks: they finish within microseconds - // and a condition-variable sleep here would cost more than the work. - // Every claimed chunk is counted, so once remaining reaches zero no - // thread is inside the callback and the stack-owned task may go. - while (remaining.load(std::memory_order_acquire) != 0) - std::this_thread::yield(); + return found ? Hit{true, best, bestId} : Hit{}; } private: - static constexpr size_t kMaximumChunks = size_t(1) << 31; - static uint64_t countOf(uint64_t ticket) { return ticket >> 32; } - static uint64_t indexOf(uint64_t ticket) { return ticket & 0xffffffffu; } - void work() { - for (;;) { - const uint64_t claim = ticket.fetch_add(1, std::memory_order_acq_rel); - const uint64_t index = indexOf(claim); - if (index >= countOf(claim)) return; - (*job)(size_t(index)); - remaining.fetch_sub(1, std::memory_order_acq_rel); + void see(uint32_t vertex) { + if (vertexMark[vertex] == vertexStamp) + return; + vertexMark[vertex] = vertexStamp; + const Point point = (*points)[vertex]; + const double dx = point.x - ox, dy = point.y - oy; + vertexAngle[vertex] = angleOf(dx, dy); + vertexDistance[vertex] = std::sqrt(dx * dx + dy * dy); + } + + void push(const Node *node, double key) { + size_t i = heapNodes.size(); + heapNodes.push_back(node); + heapKeys.push_back(key); + while (i > 0) { + const size_t parent = (i - 1) >> 1; + if (heapKeys[parent] <= key) + break; + heapNodes[i] = heapNodes[parent]; + heapKeys[i] = heapKeys[parent]; + i = parent; } + heapNodes[i] = node; + heapKeys[i] = key; } - bool pending() const { - const uint64_t current = ticket.load(std::memory_order_acquire); - return indexOf(current) < countOf(current); - } - void loop() { - for (;;) { - // A query issues three runs a few hundred microseconds apart, and a - // drag issues a query every frame. Spin briefly before sleeping so the - // next run finds the workers awake; sleep for real between frames. - bool ready = false; - for (int spin = 0; spin < 4000 && !ready; ++spin) { - ready = pending() || stop.load(std::memory_order_acquire); - if (!ready) std::this_thread::yield(); + + // Removes the nearest node and its key. + const Node *pop(double &popped) { + const Node *top = heapNodes[0]; + popped = heapKeys[0]; + const Node *node = heapNodes.back(); + const double key = heapKeys.back(); + heapNodes.pop_back(); + heapKeys.pop_back(); + const size_t n = heapNodes.size(); + if (n > 0) { + size_t i = 0; + for (;;) { + size_t child = 2 * i + 1; + if (child >= n) + break; + if (child + 1 < n && heapKeys[child + 1] < heapKeys[child]) + ++child; + if (heapKeys[child] >= key) + break; + heapNodes[i] = heapNodes[child]; + heapKeys[i] = heapKeys[child]; + i = child; } - if (!ready) { - std::unique_lock lock(mutex); - wake.wait(lock, [this] { return stop.load(std::memory_order_acquire) || pending(); }); + heapNodes[i] = node; + heapKeys[i] = key; + } + return top; + } + + void walk(const Node *root, double range) { + heapNodes.clear(); + heapKeys.clear(); + push(root, distanceTo(root->bounds, ox, oy)); + // Depths are proven from the edges filed so far, again each time the + // walk passes four times as far from the eye. + double wave = std::max(range / 16, 1e-3); + size_t proven = 0; + while (!heapNodes.empty()) { + double distance; + const Node *node = pop(distance); + if (distance > range) + break; + if (distance > wave) { + if (count() > proven) { + walkChains(proven); + raisePeaks(); + proven = count(); + } + while (wave < distance) + wave *= 4; + } + ++nodes; + if (!boundsMayShow(node->bounds, distance)) + continue; + if (!node->ids.empty()) { + for (uint32_t id : node->ids) + fileEdge(id, range); + } else { + const Node *left = node->left.get(), *right = node->right.get(); + const double leftKey = distanceTo(left->bounds, ox, oy); + if (leftKey <= range) + push(left, leftKey); + const double rightKey = distanceTo(right->bounds, ox, oy); + if (rightKey <= range) + push(right, rightKey); } - if (stop.load(std::memory_order_acquire)) return; - work(); } + if (count() > proven) + walkChains(proven); } - std::vector threads; - std::mutex mutex; - std::condition_variable wake; - const std::function *job = nullptr; - std::atomic ticket{0}; - std::atomic remaining{0}; - std::atomic stop{false}; -}; -unsigned poolWorkers() { - if (const char *override = std::getenv("ICARUS_HEIGHT_THREADS")) { - const long value = std::strtol(override, nullptr, 10); - if (value >= 0 && value <= 64) return unsigned(value); + // Whether anything in bounds, distance from the eye at its nearest, could + // show in the cone. + bool boundsMayShow(const Bounds &bounds, double distance) const { + const int column = ox < bounds.left ? 0 : (ox > bounds.right ? 2 : 1); + const int row = oy < bounds.top ? 0 : (oy > bounds.bottom ? 2 : 1); + if (distance < absoluteMargin || (column == 1 && row == 1)) + return true; + // Seen from outside, a box spans less than half a turn, between the two + // corners that bound its silhouette. Which two depends only on where the + // eye is around the box: per region, each corner as (right?, bottom?). + static constexpr bool silhouette[3][3][4] = { + {{1, 0, 0, 1}, {0, 0, 1, 0}, {0, 0, 1, 1}}, + {{0, 0, 0, 1}, {0, 0, 0, 0}, {1, 0, 1, 1}}, + {{0, 0, 1, 1}, {0, 1, 1, 1}, {1, 0, 0, 1}}}; + const bool *corners = silhouette[row][column]; + const double ax = (corners[0] ? bounds.right : bounds.left) - ox; + const double ay = (corners[1] ? bounds.bottom : bounds.top) - oy; + const double bx = (corners[2] ? bounds.right : bounds.left) - ox; + const double by = (corners[3] ? bounds.bottom : bounds.top) - oy; + const auto [first, last] = + span(ax, ay, bx, by, angleOf(ax, ay), angleOf(bx, by)); + return mayShow(first - cornerSlack, last + cornerSlack, distance); + } + + // The angles the eye sees between points a and b, at angles ta and tb: + // first to last, last unwrapped past first. + static std::pair span(double ax, double ay, double bx, + double by, double ta, double tb) { + const double turn = ax * by - ay * bx; + double first, last; + if (turn > 0) { + first = ta; + last = tb; + } else if (turn < 0) { + first = tb; + last = ta; + } else { + first = std::min(ta, tb); + last = std::max(ta, tb); + } + if (last < first) + last += 2 * pi; + return {first, last}; } - const unsigned cores = std::thread::hardware_concurrency(); - return cores > 2 ? std::min(5u, cores - 1) : 0; -} + + // Whether something distance from the eye, spanning angles first to last + // (which may run past pi), falls in the cone and is not provably behind + // the depths there. + bool mayShow(double first, double last, double distance) const { + constexpr double margin = 1e-6; + const double highest = lowest + binCount * width; + for (double shift = 2 * pi; shift > -4 * pi; shift -= 2 * pi) { + const double start = first + shift, end = last + shift; + if (end < lowest - margin || start > highest + margin) + continue; + if (!beyond(distance, peak(binOf(start), binOf(end)))) + return true; + } + return false; + } + + // The largest depth over bins first to last. + double peak(uint32_t first, uint32_t last) const { + double result = 0; + size_t low = first + leaves, high = last + leaves + 1; + while (low < high) { + if (low & 1) + result = std::max(result, peaks[low++]); + if (high & 1) + result = std::max(result, peaks[--high]); + low >>= 1; + high >>= 1; + } + return result; + } + + void raisePeaks() { + std::copy(depth.begin(), depth.begin() + binCount, + peaks.begin() + leaves); + for (size_t i = leaves - 1; i > 0; --i) + peaks[i] = std::max(peaks[2 * i], peaks[2 * i + 1]); + } + + void fileEdge(uint32_t id, double range) { + const Edge &edge = (*edges)[id]; + if (!active[edge.wall]) + return; + const double ax = edge.a.x - ox, ay = edge.a.y - oy; + const double bx = edge.b.x - ox, by = edge.b.y - oy; + const double ex = edge.b.x - edge.a.x, ey = edge.b.y - edge.a.y; + const double lengthSquared = ex * ex + ey * ey; + double t = lengthSquared == 0 ? 0.0 : -(ax * ex + ay * ey) / lengthSquared; + t = t < 0 ? 0.0 : (t > 1 ? 1.0 : t); + const double cx = ax + ex * t, cy = ay + ey * t; + const double distance = std::sqrt(cx * cx + cy * cy); + if (distance > range) + return; + see(edge.aVertex); + see(edge.bVertex); + double first = -pi, last = pi; + if (distance >= absoluteMargin) { + std::tie(first, last) = span(ax, ay, bx, by, vertexAngle[edge.aVertex], + vertexAngle[edge.bVertex]); + if (!mayShow(first - cornerSlack, last + cornerSlack, distance)) + return; + slotMark[id] = mark; + } + edgeIds.push_back(id); + near.push_back(distance); + from.push_back(first); + to.push_back(last); + } + + // A boundary index of a whole turn's bins, brought into [0, binCount): + // a run's unwrapped angle stays within a few turns. + int64_t wrap(int64_t j) const { + while (j < 0) + j += binCount; + while (j >= binCount) + j -= binCount; + return j; + } + + bool chained(int64_t id) const { + return id >= 0 && id < int64_t(edges->size()) && slotMark[size_t(id)] == mark; + } + + // Lowers each bin's depth to the farthest point of any run of consecutive + // ring edges that crosses the bin from one boundary to the next. Such a run + // is an unbroken wall across the bin, so it stops every ray in it. Single + // edges rarely span a bin: Riot's outlines are many short strokes. + void walkChains(size_t firstNew) { + const std::vector &list = *edges; + // Runs already walked have proven all they can; walk again only those + // with an edge filed since, each from its first edge. + advance(dirtyStamp, dirty); + const uint32_t stamp = dirtyStamp; + heads.clear(); + for (size_t slot = firstNew; slot < count(); ++slot) { + uint32_t id = edgeIds[slot]; + if (slotMark[id] != mark) + continue; + while (dirty[id] != stamp) { + dirty[id] = stamp; + if (!chained(int64_t(id) - 1) || list[id - 1].bVertex != list[id].aVertex) { + heads.push_back(id); + break; + } + --id; + } + } + const int64_t bins = binCount; + for (uint32_t head : heads) { + uint32_t id = head; + double raw = vertexAngle[list[id].aVertex]; + // The run's angle, unwrapped so it never jumps by a turn. + double angle = raw; + // The last boundary crossed, and since then: the farthest point and the + // angles the run has reached. + bool hasLast = false; + int64_t last = 0; + double farthest = 0, lowAngle = 0, highAngle = 0; + for (;;) { + const Edge &edge = list[id]; + const uint32_t next = edge.bVertex; + const double nextRaw = vertexAngle[next]; + double turn = nextRaw - raw; + if (turn > pi) + turn -= 2 * pi; + else if (turn <= -pi) + turn += 2 * pi; + const double end = angle + turn; + const double ex = edge.b.x - edge.a.x, ey = edge.b.y - edge.a.y; + const double along = (edge.a.x - ox) * ey - (edge.a.y - oy) * ex; + // Boundaries the edge plainly crosses, in the order it crosses them. + const int64_t step = end > angle ? 1 : -1; + double start = step > 0 + ? std::ceil((angle + angleMargin - lowest) / width) + : std::floor((angle - angleMargin - lowest) / width); + double finish = step > 0 + ? std::floor((end - angleMargin - lowest) / width) + : std::ceil((end + angleMargin - lowest) / width); + if (!whole) { + // Boundaries outside the cone are skipped below and leave the run + // alone, so only those just outside need walking. A narrow cone's + // bins are so fine the others can lie past any integer. + start = boundaryIndex(start, double(bins)); + finish = boundaryIndex(finish, double(bins)); + } + const int64_t stop = int64_t(finish); + for (int64_t j = int64_t(start); step > 0 ? j <= stop : j >= stop; + j += step) { + if (!whole && (j < 0 || j > bins)) + continue; + const int64_t boundary = whole ? wrap(j) : j; + const double distance = + along / (boundaryX[size_t(boundary)] * ey - + boundaryY[size_t(boundary)] * ex); + const double at = lowest + double(j) * width; + if (!(distance >= 0 && distance < infinity)) { + hasLast = false; + continue; + } + if (hasLast && std::abs(j - last) == 1) { + const double previous = lowest + double(last) * width; + // Between the two crossings the run stayed inside the bin. + if (lowAngle >= std::min(at, previous) - 1e-8 && + highAngle <= std::max(at, previous) + 1e-8) { + const int64_t low = std::min(j, last); + const size_t k = size_t(whole ? wrap(low) : low); + const double bound = std::max(farthest, distance); + if (bound < depth[k]) + depth[k] = bound; + } + } + hasLast = true; + last = j; + farthest = distance; + lowAngle = highAngle = at; + } + farthest = std::max(farthest, vertexDistance[next]); + lowAngle = std::min(lowAngle, end); + highAngle = std::max(highAngle, end); + if (!chained(int64_t(id) + 1) || list[id + 1].aVertex != next) + break; + ++id; + raw = nextRaw; + angle = end; + } + } + } + + // Counts each slot into offsets[k + 1] for every bin k it is filed in, or, + // with fill, writes it at offsets[k] and advances that. + void place(bool fill) { + const double highest = lowest + binCount * width; + const uint32_t slots = uint32_t(count()); + for (uint32_t slot = 0; slot < slots; ++slot) { + const double distance = near[slot]; + if (distance < absoluteMargin) { + // An edge through the eye can stop a ray in any direction. + for (uint32_t k = 0; k < binCount; ++k) { + if (fill) + items[offsets[k]++] = slot; + else + ++offsets[k + 1]; + } + continue; + } + // Rays are admitted a sliver past an edge's ends; see Edge::intersection. + const double pad = angleMargin + 1e-12 / distance; + // A span starting just below -pi also covers the bins just below pi. + for (double shift = 2 * pi; shift > -4 * pi; shift -= 2 * pi) { + const double start = from[slot] + shift - pad; + const double end = to[slot] + shift + pad; + if (end < lowest) + break; + if (start > highest) + continue; + const uint32_t first = binOf(start), last = binOf(end); + for (uint32_t k = first; k <= last; ++k) { + if (beyond(distance, depth[k])) + continue; + if (fill) + items[offsets[k]++] = slot; + else + ++offsets[k + 1]; + } + } + } + } +}; struct Handle { std::vector edges; @@ -346,16 +779,15 @@ struct Handle { std::vector output; std::vector activeScratch; ISHResult resultScratch{}; - std::vector angles, arcAngles, vertexAngles; - std::vector arcHits, hits; - std::vector chunkCounters; - std::vector> chunkAngles, chunkVertexAngles; - std::vector candidates; + std::vector angles, vertexAngles; // The edges meeting at each vertex, to tell a corner from a seam. std::vector> vertexEdges; std::vector vertexPoints; + ConeBins bins; + // Vertices already turned into events this query: eventMark is eventStamp. + std::vector eventMark; + uint32_t eventStamp = 0; bool interiorSides = false; - Pool pool{poolWorkers()}; std::mutex mutex; std::string error; @@ -386,11 +818,10 @@ struct Handle { vertexEdges[edges[i].aVertex].push_back(i); vertexEdges[edges[i].bVertex].push_back(i); } - candidates.reserve(edges.size()); + bins.reserve(vertices.size(), edges.size()); + eventMark.assign(vertices.size(), 0); angles.reserve(std::min(maximumPoints, size_t(4097) + vertices.size() * 3)); - arcAngles.reserve(maximumArcSteps + 1); vertexAngles.reserve(vertices.size()); - arcHits.reserve(maximumArcSteps + 1); if (!edges.empty()) { std::vector ids(edges.size()); for (uint32_t i = 0; i < ids.size(); ++i) @@ -486,81 +917,6 @@ bool passThrough(const Handle &handle, uint32_t vertex, Point delta, return count == 2 && sides[0] * sides[1] < 0; } -Hit castRay(const Handle &handle, Point origin, Point direction, double range, - const uint8_t *active, Counters &counters) { - Hit result; - if (!handle.tree || range == 0) - return result; - double best = range; - const PreparedRay ray{origin, direction, - direction.x == 0 ? 0 : 1 / direction.x, - direction.y == 0 ? 0 : 1 / direction.y}; - // A balanced binary tree with at most 2^19 leaves needs fewer than 64 - // pending siblings. Keep ray traversal off the allocator hot path. - struct Pending { - const Node *node; - double entryDistance; - }; - std::array stack{}; - size_t stackSize = 0; - double rootEntry; - if (entry(handle.tree->bounds, ray, best, rootEntry)) - stack[stackSize++] = {handle.tree.get(), rootEntry}; - else - ++counters.nodes; - while (stackSize) { - const Pending pending = stack[--stackSize]; - const Node *node = pending.node; - ++counters.nodes; - if (pending.entryDistance > best) - continue; - if (!node->ids.empty()) { - for (uint32_t id : node->ids) { - const Edge &edge = handle.edges[id]; - if (!active[edge.wall]) - continue; - ++counters.edgeTests; - double distance; - if (edge.intersection(origin, direction, best, distance) && - (distance < best || !result.found)) { - best = distance; - result = {true, distance, id}; - } - } - continue; - } - double leftEntry, rightEntry; - const bool hasLeft = entry(node->left->bounds, ray, best, leftEntry); - const bool hasRight = entry(node->right->bounds, ray, best, rightEntry); - if (hasLeft && hasRight) { - if (leftEntry <= rightEntry) { - stack[stackSize++] = {node->right.get(), rightEntry}; - stack[stackSize++] = {node->left.get(), leftEntry}; - } else { - stack[stackSize++] = {node->left.get(), leftEntry}; - stack[stackSize++] = {node->right.get(), rightEntry}; - } - } else if (hasLeft) { - stack[stackSize++] = {node->left.get(), leftEntry}; - } else if (hasRight) { - stack[stackSize++] = {node->right.get(), rightEntry}; - } - } - return result; -} - -void collectCandidates(const Node *node, const Bounds &area, - std::vector &output) { - if (!overlaps(node->bounds, area)) - return; - if (!node->ids.empty()) { - output.insert(output.end(), node->ids.begin(), node->ids.end()); - } else { - collectCandidates(node->left.get(), area, output); - collectCandidates(node->right.get(), area, output); - } -} - void copyText(const std::string &text, char *target, uint32_t capacity) { if (!target || capacity == 0) return; @@ -651,7 +1007,6 @@ int32_t ish_query(void *opaque, double originX, double originY, out->status = ISH_BUSY; return ISH_BUSY; } - Boost boost; const auto started = Clock::now(); try { const double values[] = {originX, originY, directionRadians, range, @@ -665,188 +1020,102 @@ int32_t ish_query(void *opaque, double originX, double originY, return failure(handle, out, ISH_INVALID, "invalid SVG query range, aperture, arc steps or mask"); + // Every ray starts at the eye, so the walls near it are filed by the + // angle they cover (see ConeBins): a ray tests only the few edges in its + // bin, and a corner that a nearer wall provably hides casts no rays. const Point origin{originX, originY}; const double half = apertureRadians / 2; + const bool whole = half + 1e-6 >= pi; + ConeBins &bins = handle.bins; + bins.file(handle.edges, handle.vertexPoints, handle.tree.get(), origin, + directionRadians, range, active, whole ? -pi : -half, + whole ? 2 * pi : apertureRadians); + auto &angles = handle.angles; - auto &arcAngles = handle.arcAngles; auto &vertexAngles = handle.vertexAngles; - auto &arcHits = handle.arcHits; angles.clear(); - arcAngles.clear(); vertexAngles.clear(); - arcHits.clear(); - Counters counters; - constexpr size_t chunkSize = 32; - constexpr size_t eventChunk = 256; - for (uint32_t i = 0; i <= arcSteps; ++i) { - const double angle = -half + apertureRadians * i / arcSteps; - angles.push_back(angle); - arcAngles.push_back(angle); - } - arcHits.resize(arcAngles.size()); - { - const size_t chunks = (arcAngles.size() + chunkSize - 1) / chunkSize; - handle.chunkCounters.assign(chunks, ChunkCounters{}); - handle.pool.run(chunks, [&](size_t chunk) { - Counters local; - const size_t end = std::min(arcAngles.size(), (chunk + 1) * chunkSize); - for (size_t i = chunk * chunkSize; i < end; ++i) { - const double world = directionRadians + arcAngles[i]; - arcHits[i] = castRay(handle, origin, {std::cos(world), std::sin(world)}, - range, active, local); - } - handle.chunkCounters[chunk].value = local; - }); - for (const ChunkCounters &local : handle.chunkCounters) { - counters.edgeTests += local.value.edgeTests; - counters.nodes += local.value.nodes; - } + for (uint32_t i = 0; i <= arcSteps; ++i) + angles.push_back(-half + apertureRadians * i / arcSteps); + // Rays aimed at a vertex or at a wall's crossing of the range circle stay + // in the outline even when their neighbours meet the same edge. + auto add = [&](double event, bool vertex) { + // Behind a full circle, a ray beside the seam wraps round to the other + // end: -pi and pi are the same direction. + if (whole && event < -half) event += 2 * pi; + if (whole && event > half) event -= 2 * pi; + if (event < -half || event > half) return; + angles.push_back(event); + if (vertex) vertexAngles.push_back(event); + }; + auto emit = [&](double angle, bool beside) { + if (beside) add(angle - cornerOffset, false); + add(angle, true); + if (beside) add(angle + cornerOffset, false); + }; + // Whether an event at angle could start a ray in the aperture and is not + // provably behind a nearer wall. + auto open = [&](double angle, double distance) { + if (!whole && (angle < -half - 1e-6 || angle > half + 1e-6)) + return false; + return !bins.hides(angle, distance); + }; + + const double rangeSquared = range * range; + // Overlapping painted strokes create visibility corners at their crossing. + for (const Crossing &crossing : handle.crossings) { + if (!active[crossing.first] || !active[crossing.second]) continue; + const Point delta = crossing.point - origin; + const double distanceSquared = dot(delta, delta); + if (distanceSquared > rangeSquared || (delta.x == 0 && delta.y == 0)) + continue; + const double angle = bins.angleOf(delta.x, delta.y); + if (!open(angle, std::sqrt(distanceSquared))) continue; + emit(angle, true); } const auto prepared = Clock::now(); - auto &candidates = handle.candidates; - candidates.clear(); - if (handle.tree) { - const Bounds area{origin.x - range, origin.y - range, origin.x + range, - origin.y + range}; - collectCandidates(handle.tree.get(), area, candidates); - } - const double rangeSquared = range * range; - // Sector culling: a vertex outside the aperture (with slack for the - // corner offsets) cannot start a ray inside it, so skip its trig. - const Point facing{std::cos(directionRadians), std::sin(directionRadians)}; - const bool cullSector = half + 1e-6 < pi; - const double cosSlack = std::cos(std::min(pi, half + 1e-6)); - auto inSector = [&](Point delta) { - if (!cullSector) return true; - const double length = std::sqrt(dot(delta, delta)); - return dot(delta, facing) >= cosSlack * length; - }; - auto hiddenEvent = [&](double angle, Point delta) { - if (apertureRadians / arcSteps >= pi || angle <= -half || angle >= half) return false; - int interval=int(std::floor((angle+half)/apertureRadians*double(arcSteps))); - interval=std::max(0,std::min(interval,int(arcSteps)-1)); - const Hit &first=arcHits[size_t(interval)], &last=arcHits[size_t(interval)+1]; - if (!first.found || !last.found || first.edge!=last.edge || - angle-cornerOffsetarcAngles[size_t(interval)+1]) return false; - const double distance=std::sqrt(dot(delta,delta)); - double hitDistance; - return handle.edges[first.edge].intersection(origin,{delta.x/distance,delta.y/distance},distance,hitDistance) - && hitDistance localAngles = std::move(handle.chunkAngles[chunk]); - std::vector localVertex = std::move(handle.chunkVertexAngles[chunk]); - localAngles.clear(); - localVertex.clear(); - struct Store { - std::vector &angles, &vertex, &outAngles, &outVertex; - ~Store() { outAngles = std::move(angles); outVertex = std::move(vertex); } - } store{localAngles, localVertex, handle.chunkAngles[chunk], handle.chunkVertexAngles[chunk]}; - auto emit = [&](double angle, bool vertex, bool beside = true) { - for (double event : {angle - cornerOffset, angle, angle + cornerOffset}) { - if (event != angle && !beside) continue; - if (event >= -half && event <= half) { - localAngles.push_back(event); - if (vertex && event == angle) localVertex.push_back(event); - } - } - }; - if (chunk < crossingChunks) { - const size_t end = std::min(handle.crossings.size(), (chunk + 1) * eventChunk); - for (size_t c = chunk * eventChunk; c < end; ++c) { - const Crossing &crossing = handle.crossings[c]; - if (!active[crossing.first] || !active[crossing.second]) continue; - const Point delta = crossing.point - origin; - if (dot(delta, delta) > rangeSquared || (delta.x == 0 && delta.y == 0)) continue; - if (!inSector(delta)) continue; - const double relative = std::atan2(delta.y, delta.x) - directionRadians; - const double angle = std::atan2(std::sin(relative), std::cos(relative)); - if (hiddenEvent(angle, delta)) continue; - emit(angle, true); - } - return; - } - const size_t first = (chunk - crossingChunks) * eventChunk; - const size_t end = std::min(candidates.size(), first + eventChunk); - for (size_t c = first; c < end; ++c) { - const uint32_t id = candidates[c]; + advance(handle.eventStamp, handle.eventMark); + const uint32_t stamp = handle.eventStamp; + for (size_t slot = 0; slot < bins.count(); ++slot) { + const uint32_t id = bins.edgeIds[slot]; const Edge &edge = handle.edges[id]; - if (!active[edge.wall]) - continue; - // Preserve the exact wall/range-circle transition, even when neither - // endpoint lies in range. Otherwise a polygon chord clips the wall early. + // A long wall can cross the range circle without either endpoint being + // in range. Seed that exact transition so the polygon follows the wall + // all the way to the circle instead of cutting diagonally short of it. const Point segment = edge.b - edge.a; - const Point relativeStart = edge.a - origin; + const Point relative = edge.a - origin; const double lengthSquared = dot(segment, segment); - const double projection = -dot(relativeStart, segment) / lengthSquared; - const Point closest{relativeStart.x + segment.x * projection, - relativeStart.y + segment.y * projection}; + const double projection = -dot(relative, segment) / lengthSquared; + const Point closest{relative.x + segment.x * projection, + relative.y + segment.y * projection}; const double remaining = rangeSquared - dot(closest, closest); if (remaining >= 0) { const double offset = std::sqrt(remaining / lengthSquared); for (double t : {projection - offset, projection + offset}) { if (t < 0 || t > 1) continue; - const Point delta{relativeStart.x + segment.x * t, - relativeStart.y + segment.y * t}; - if (!inSector(delta)) continue; - const double relative = std::atan2(delta.y, delta.x) - directionRadians; - const double angle = std::atan2(std::sin(relative), std::cos(relative)); - if (angle >= -half && angle <= half && !hiddenEvent(angle,delta)) { - localAngles.push_back(angle); - localVertex.push_back(angle); - } - } - } - const Point endpoints[] = {edge.a, edge.b}; - const uint32_t endpointVertices[] = {edge.aVertex, edge.bVertex}; - for (int endpoint = 0; endpoint < 2; ++endpoint) { - if (seam(handle, endpointVertices[endpoint], active)) continue; - const Point point = endpoints[endpoint]; - const Point delta = point - origin; - const double distanceSquared = dot(delta, delta); - if (distanceSquared > rangeSquared || - (delta.x == 0 && delta.y == 0)) - continue; - if (!inSector(delta)) continue; - const double relative = std::atan2(delta.y, delta.x) - directionRadians; - const double angle = std::atan2(std::sin(relative), std::cos(relative)); - if (apertureRadians / arcSteps < pi && angle > -half && angle < half) { - int interval = int(std::floor((angle + half) / apertureRadians * - double(arcSteps))); - interval = std::max(0, std::min(interval, int(arcSteps) - 1)); - const Hit &first = arcHits[size_t(interval)]; - const Hit &last = arcHits[size_t(interval) + 1]; - if (first.found && last.found && first.edge == last.edge && - angle - cornerOffset >= arcAngles[size_t(interval)] && - angle + cornerOffset <= arcAngles[size_t(interval) + 1]) { - const double distance = std::sqrt(distanceSquared); - double hitDistance; - const Point ray{delta.x / distance, delta.y / distance}; - if (handle.edges[first.edge].intersection(origin, ray, distance, - hitDistance) && - hitDistance < distance - 1e-7) - continue; + const double angle = bins.angleOf(relative.x + segment.x * t, + relative.y + segment.y * t); + if (!open(angle, range)) continue; + if (angle >= -half && angle <= half) { + angles.push_back(angle); + vertexAngles.push_back(angle); } } - emit(angle, true, - !passThrough(handle, endpointVertices[endpoint], delta, active)); } + for (uint32_t vertex : {edge.aVertex, edge.bVertex}) { + if (handle.eventMark[vertex] == stamp) continue; + handle.eventMark[vertex] = stamp; + const double distance = bins.vertexDistance[vertex]; + if (distance > range || distance == 0) continue; + const double angle = bins.vertexAngle[vertex]; + if (!open(angle, distance)) continue; + if (seam(handle, vertex, active)) continue; + // Rays just beside a vertex the wall runs straight across meet its + // two edges, so only the vertex ray adds a corner. + emit(angle, !passThrough(handle, vertex, + handle.vertexPoints[vertex] - origin, active)); } - }); - for (size_t chunk = 0; chunk < eventChunks; ++chunk) { - angles.insert(angles.end(), handle.chunkAngles[chunk].begin(), handle.chunkAngles[chunk].end()); - vertexAngles.insert(vertexAngles.end(), handle.chunkVertexAngles[chunk].begin(), handle.chunkVertexAngles[chunk].end()); } std::sort(angles.begin(), angles.end()); angles.erase(std::unique(angles.begin(), angles.end()), angles.end()); @@ -862,40 +1131,14 @@ int32_t ish_query(void *opaque, double originX, double originY, handle.output.reserve((angles.size() + 1) * 2); handle.output.push_back(origin.x); handle.output.push_back(origin.y); - auto &hits = handle.hits; - hits.resize(angles.size()); - { - const size_t chunks = (angles.size() + chunkSize - 1) / chunkSize; - handle.chunkCounters.assign(chunks, ChunkCounters{}); - handle.pool.run(chunks, [&](size_t chunk) { - Counters local; - const size_t end = std::min(angles.size(), (chunk + 1) * chunkSize); - for (size_t i = chunk * chunkSize; i < end; ++i) { - const double angle = angles[i]; - const auto found = std::lower_bound(arcAngles.begin(), arcAngles.end(), angle); - if (found != arcAngles.end() && *found == angle) { - hits[i] = arcHits[size_t(found - arcAngles.begin())]; - } else { - const double world = directionRadians + angle; - hits[i] = castRay(handle, origin, {std::cos(world), std::sin(world)}, - range, active, local); - } - } - handle.chunkCounters[chunk].value = local; - }); - for (const ChunkCounters &local : handle.chunkCounters) { - counters.edgeTests += local.value.edgeTests; - counters.nodes += local.value.nodes; - } - } + uint64_t edgeTests = 0; Hit previousHit; size_t sameEdgeRun = 0; Point runAnchor{}; - for (size_t rayIndex = 0; rayIndex < angles.size(); ++rayIndex) { - const double angle = angles[rayIndex]; + for (const double angle : angles) { const double world = directionRadians + angle; const Point direction{std::cos(world), std::sin(world)}; - const Hit hit = hits[rayIndex]; + const Hit hit = bins.cast(origin, direction, angle, range, edgeTests); const double distance = hit.found ? hit.distance : range; const double x = origin.x + direction.x * distance; const double y = origin.y + direction.y * distance; @@ -942,9 +1185,10 @@ int32_t ish_query(void *opaque, double originX, double originY, out->points = handle.output.data(); out->pointCount = uint32_t(handle.output.size() / 2); out->rayCount = uint32_t(angles.size()); - out->edgeTests = counters.edgeTests; - out->spatialNodes = counters.nodes; - out->candidateEdges = uint32_t(candidates.size()); + out->edgeTests = edgeTests; + // Tree nodes the filing walk visited, and one bin per ray. + out->spatialNodes = bins.nodes + angles.size(); + out->candidateEdges = uint32_t(bins.count()); out->preparationMicros = std::chrono::duration(prepared - started).count(); out->candidateMicros = std::chrono::duration( diff --git a/native/height/svg_height_native_test.cpp b/native/height/svg_height_native_test.cpp index 3c9ac0de..fe8107f2 100644 --- a/native/height/svg_height_native_test.cpp +++ b/native/height/svg_height_native_test.cpp @@ -45,6 +45,21 @@ int main() { CHECK(std::abs(result->points[4] - 5) < 1e-12); CHECK(std::abs(result->points[5] - 5) < 1e-12); + // A cone this narrow has bins far finer than any integer count of them + // across the wall's angles; the wall still stops its rays. + result->structSize = sizeof(*result); + CHECK(ish_query(handle, 0, 0, 0, 10, 1e-300, 2, active, 1, result) == + ISH_OK); + CHECK(result->pointCount >= 2); + CHECK(std::abs(result->points[2] - 5) < 1e-12); + + // At the smallest positive aperture the bins have no width at all. + result->structSize = sizeof(*result); + CHECK(ish_query(handle, 0, 0, 0, 10, 4.9406564584124654e-324, 2, active, 1, + result) == ISH_OK); + CHECK(result->pointCount >= 2); + CHECK(std::abs(result->points[2] - 5) < 1e-12); + active[0] = 0; result->structSize = sizeof(*result); CHECK(ish_query(handle, 0, 0, 0, 10, 1.5707963267948966, 2, active, diff --git a/test/svg_cone_exact_test.dart b/test/svg_cone_exact_test.dart new file mode 100644 index 00000000..99fed0c8 --- /dev/null +++ b/test/svg_cone_exact_test.dart @@ -0,0 +1,291 @@ +import 'dart:convert'; +import 'dart:io'; +import 'dart:math' as math; + +import 'package:flutter_test/flutter_test.dart'; +import 'package:icarus/view_cone/svg_height_visibility.dart'; + +/// A cone's outline is a fan of rays from the eye joined by straight lines. +/// The cone query skips every ray it can prove would land in the middle of a +/// wall another ray already covers. If it ever skips a real corner, the line +/// joining its neighbours cuts across a wall or a gap. So in every direction +/// a single exact ray must end where the outline does: on the wall it meets, +/// or, in open sky, on the chord the outline draws along the range circle. +/// +/// Directions come from the outline (between each pair of its points) and, +/// independently of it, from just either side of every wall corner in range +/// and an even sweep, so a wall the query missed entirely is still probed. +SvgHeightVisibility _model(Map json, {String? nativeLibrary}) { + final model = SvgHeightVisibility.fromJson(json); + if (nativeLibrary != null) { + expect(model.enableNativeAcceleration(libraryPath: nativeLibrary), isTrue); + } + return model; +} + +Map _asset(String name) => jsonDecode(utf8.decode( + gzip.decode(File('assets/maps/$name.json.gz').readAsBytesSync()))) + as Map; + +List _rectangle(double x0, double y0, double x1, double y1) => + [x0, y0, x1, y0, x1, y1, x0, y1]; + +Map _walls(List> rings) => { + 'version': 1, + 'coordinateSpace': 'svg', + 'verticalSpace': 'meters-above-local-floor', + 'walls': [ + for (var i = 0; i < rings.length; i++) + { + 'id': 'wall-$i', + 'rings': [rings[i]], + 'bands': [ + [0, null] + ], + 'unknownHeight': false, + 'fillRule': 'nonzero', + } + ], + 'supports': const [], + 'receiver': const [], + 'cellSizeSvg': 4, + }; + +double _relative(Offset delta, double direction) { + final turn = math.atan2(delta.dy, delta.dx) - direction; + return math.atan2(math.sin(turn), math.cos(turn)); +} + +/// The places where [model]'s cone outline disagrees with exact rays. +List _wrong(SvgHeightVisibility model, Offset eye, double direction, + double aperture, double range) { + final outline = model + .horizontalCone( + origin: eye, + directionRadians: direction, + range: range, + apertureRadians: aperture, + supportId: model.automaticSupportAt(eye)?.id) + .polygon; + return _wrongOutline(model, outline, eye, direction, aperture, range); +} + +/// The places where [outline] disagrees with exact rays through [model]. +List _wrongOutline(SvgHeightVisibility model, List outline, + Offset eye, double direction, double aperture, double range) { + final support = model.automaticSupportAt(eye)?.id; + final points = outline.skip(1).toList(); + if (points.length < 2) return ['no outline']; + final half = aperture / 2; + // Angles in ray order, from the cone's first edge: on a full circle the + // first and last points lie on the same line, at -half and half. + double turn(double angle) => math.atan2(math.sin(angle), math.cos(angle)); + final angles = [ + -half + turn(_relative(points.first - eye, direction) + half) + ]; + for (final p in points.skip(1)) { + // Rays come in increasing angle, so each step is forward; a hair + // backward is rounding. + var step = (_relative(p - eye, direction) - angles.last) % (2 * math.pi); + if (step > 2 * math.pi - 1e-9) step -= 2 * math.pi; + angles.add(angles.last + step); + } + // The last point is the ray along the cone's far edge, a full turn on + // from the first on a full circle even where the two coincide. + if (angles.last < half - math.pi) angles.last += 2 * math.pi; + final probes = [ + for (var i = 0; i + 1 < angles.length; i++) (angles[i] + angles[i + 1]) / 2, + for (final wall in model.walls) + for (final ring in wall.rings) + for (final corner in ring) + if ((corner - eye).distance <= range) ...[ + _relative(corner - eye, direction) - 1e-6, + _relative(corner - eye, direction) + 1e-6, + ], + for (var i = 0; i < 200; i++) -half + aperture * (i + 0.5) / 200, + ]; + final wrong = []; + // Every point lies where the ray toward it ends. This needs no probe + // between points, so it also holds a cone too narrow to probe. + for (final p in points) { + final delta = p - eye; + if (delta.distance < 1e-9) continue; + final hit = model.castRay( + origin: eye, + directionRadians: math.atan2(delta.dy, delta.dx), + range: range, + supportId: support); + final truth = hit?.distance ?? range; + if ((truth - delta.distance).abs() > 1e-6) { + wrong.add('point $p at ${delta.distance.toStringAsFixed(5)}, ' + 'its ray ends at ${truth.toStringAsFixed(5)}'); + } + } + for (final angle in probes) { + if (angle <= -half || angle >= half) continue; + // The outline point pair the probe falls between, in ray order. + var after = 0; + while (after < angles.length && angles[after] < angle) { + after++; + } + if (after == 0 || after == angles.length) continue; + final gap = angles[after] - angles[after - 1]; + // Beside a silhouette the outline draws a chord between the corner ray + // and the ray 1e-8 past it, by design. + if (gap < 1e-7) continue; + final world = direction + angle; + final ray = Offset(math.cos(world), math.sin(world)); + final a = points[after - 1] - eye, b = points[after] - eye; + final edge = b - a; + final line = (a.dx * edge.dy - a.dy * edge.dx) / + (ray.dx * edge.dy - ray.dy * edge.dx); + final hit = model.castRay( + origin: eye, directionRadians: world, range: range, supportId: support); + final truth = hit?.distance ?? range; + // Open sky: the outline follows the range circle with chords, which sag + // inside it by at most this much. + final slack = hit == null ? range * (1 - math.cos(gap / 2)) + 1e-6 : 1e-6; + if (!line.isFinite || line > truth + 1e-6 || line < truth - slack) { + wrong.add('ray ${angle.toStringAsFixed(7)}: outline at ' + '${line.toStringAsFixed(5)}, ray ends at ${truth.toStringAsFixed(5)}'); + } + } + return wrong; +} + +void main() { + final nativeLibrary = Platform.environment['ICARUS_SVG_NATIVE_LIBRARY']; + + void check(String label, {String? library}) { + final skip = label == 'native' && library == null + ? 'Set ICARUS_SVG_NATIVE_LIBRARY.' + : null; + + test('$label: a wall straddling the seam behind a full circle', () { + // Found in review: an edge whose angles start just below -pi was not + // filed in the bins just below pi. + final model = _model( + _walls([ + _rectangle(10, 0, 11, 10), + _rectangle(-21, -21, -20, -20), + ]), + nativeLibrary: library); + expect(_wrong(model, Offset.zero, math.pi, 2 * math.pi, 100), isEmpty); + // The last ray, at the seam itself, meets the corner on it. + final outline = model + .horizontalCone( + origin: Offset.zero, + directionRadians: math.pi, + range: 100, + apertureRadians: 2 * math.pi) + .polygon; + expect(outline.last.distance, closeTo(10, 1e-9)); + }, skip: skip); + + test('$label: a vanishingly narrow cone returns at once', () { + // Found in review: at 5e-324 radians the bins have no width, and a + // corner a hair off the cone's direction gave NaN bin indices. + final model = _model( + _walls([ + [10, -1e-8, 11, -1e-8, 11, 10, 10, 10] + ]), + nativeLibrary: library); + for (final aperture in [1e-6, 1e-300, 5e-324]) { + expect(_wrong(model, Offset.zero, 0, aperture, 100), isEmpty, + reason: 'aperture $aperture'); + } + }, timeout: const Timeout(Duration(seconds: 20)), skip: skip); + + test('$label: a sliver of wall between two rays', () { + final model = _model(_walls([_rectangle(10, .03, 11, .05)]), + nativeLibrary: library); + expect(_wrong(model, Offset.zero, 0, 1.8, 100), isEmpty); + }, skip: skip); + + test('$label: real cones end where exact rays do', () { + // Tight spots against map edges and gaps between walls, plus an even + // spread over two busy maps. + final cases = <(String, Offset)>[ + ('haven_svg_height_defense', const Offset(76, 308)), + ('pearl_svg_height_defense', const Offset(436, 292)), + ('pearl_svg_height_attack', const Offset(76, 340)), + ('bind_svg_height_defense', const Offset(164, 276)), + ('sunset_svg_height_defense', const Offset(141.37, 428)), + for (final name in [ + 'breeze_svg_height_attack', + 'lotus_svg_height_attack' + ]) + for (var x = 20.0; x < 460; x += 48) + for (var y = 20.0; y < 460; y += 48) (name, Offset(x, y)), + ]; + final models = {}; + final failures = []; + var cones = 0; + for (final (name, at) in cases) { + final model = + models[name] ??= _model(_asset(name), nativeLibrary: library); + final eye = model.standablePointNear(at); + if (eye == null) continue; + for (final aperture in [1.8, math.pi, 2 * math.pi]) { + final direction = (at.dx * 7 + at.dy * 3 + aperture) % (2 * math.pi); + cones++; + for (final wrong in _wrong(model, eye, direction, aperture, 140)) { + failures.add('$name $eye dir $direction aperture $aperture $wrong'); + } + } + } + expect(cones, greaterThan(250)); + expect(failures.take(20), isEmpty); + }, timeout: const Timeout(Duration(minutes: 5)), skip: skip); + } + + test('the check finds a wall an outline left out', () { + // Found in review: on a full circle the first point's angle could come + // back as pi rather than -pi, and the check then probed nothing. + final wall = _walls([_rectangle(10, .03, 11, .05)]); + final real = _model(wall), empty = _model(_walls(const [])); + for (final (direction, aperture, arcSteps) in [ + (0.0, 1.8, 96), + (0.0, 2 * math.pi, 96), + (-math.pi / 2, 2 * math.pi, 96), + (math.pi, 2 * math.pi, 96), + // Found in review: steps between points wider than half a turn. + (0.0, 4.0, 1), + (0.0, 3.2, 1), + (0.0, 2 * math.pi, 1), + ]) { + final missing = empty + .horizontalCone( + origin: Offset.zero, + directionRadians: direction, + range: 100, + apertureRadians: aperture, + arcSteps: arcSteps) + .polygon; + expect( + _wrongOutline(real, missing, Offset.zero, direction, aperture, 100), + isNotEmpty, + reason: 'direction $direction aperture $aperture'); + } + // Found in review: cones too narrow to probe between their points. + final ledge = _model(_walls([ + [10, -1e-8, 11, -1e-8, 11, 10, 10, 10] + ])); + for (final aperture in [1e-6, 1e-300, 5e-324]) { + final missing = empty + .horizontalCone( + origin: Offset.zero, + directionRadians: 0, + range: 100, + apertureRadians: aperture) + .polygon; + expect(_wrongOutline(ledge, missing, Offset.zero, 0, aperture, 100), + isNotEmpty, + reason: 'aperture $aperture'); + } + }); + + check('Dart'); + // Desktop runs the same query in native/height. + check('native', library: nativeLibrary); +} diff --git a/test/svg_height_visibility_test.dart b/test/svg_height_visibility_test.dart index 276aad78..8aeac930 100644 --- a/test/svg_height_visibility_test.dart +++ b/test/svg_height_visibility_test.dart @@ -486,6 +486,29 @@ void main() { range: 100, apertureRadians: 1); + test('a vanishingly narrow cone still returns at once', () { + // Its bins are so thin that the walls beside it lie past any integer + // number of them. + final model = SvgHeightVisibility.fromJson(data([ + wall('near', [rectangle(10, -20, 11, 20)]), + wall('ring', [ + [-30, -30, 30, -30, 30, 30, -30, 30], + [-29, -29, -29, 29, 29, 29, 29, -29], + ]), + ])); + for (final aperture in [1e-6, 1e-300, 5e-324]) { + final cone = model.horizontalCone( + origin: Offset.zero, + directionRadians: 0, + range: 100, + apertureRadians: aperture); + expect(cone.polygon.length, greaterThan(1), reason: 'aperture $aperture'); + expect( + cone.polygon.skip(1).every((p) => (p.dx - 10).abs() < 1e-9), isTrue, + reason: 'aperture $aperture'); + } + }); + test('touching pieces of one wall cast the cone the whole wall casts', () { final whole = coneOf([ wall('whole', [rectangle(10, -20, 11, 20)])