Complete experimental GQA bounds and visualization

Derive finite sampling bounds from prepared region trees so bounded half-space intersections and subtractions remain analyzable. Reject unbounded selections explicitly, propagate unbounded submodels, and keep infinite solids out of the Crofton sampling sphere without removing them from Boolean evaluation. Add an unsampled bounds measurement with text and JSON reporting. Use the same conservative bounds for sampling, adaptive dimensional targets, and automatic visualization scale. Restore Plot3 output for volume, gaps, overlaps, adjacent air, and exposed air. Provide exact all-rays and bounded representative policies. Representatives retain deterministic spatial samples, draw finite grid or Crofton support footprints, and preserve diagnostic severity and region identity. Publish immutable plot batches only between parallel ray batches. Add an MGED main-thread bridge that incrementally replaces or appends the named gqa::analysis overlay while keeping all view mutation on the display-owning thread. Repair named vlblock replacement so refreshed objects preserve their identity and view ownership. Document the behavior, limitations, and supporting sampling references, including Owen's replication guidance. Add the analytic and complex model comparison report. Expand regression coverage for bounds, half-spaces, submodels, plot policies, diagnostic files, and live MGED output.
fixes
2026-09-21 16:06:58 +00:00 · 2026-09-21 06:33:46 -04:00 · 2026-09-20 12:02:59 +00:00 · 2026-09-20 12:02:59 +00:00 · 2026-09-20 12:02:59 +00:00
50 changed files with 7535 additions and 274 deletions
@@ -0,0 +1,348 @@
+# GQA sampling-method comparison
+
+## Executive result
+
+The new methods are a substantial practical improvement, but Crofton is not a
+uniform numerical improvement over grids.
+
+- **Crofton QMC is the best general-purpose stochastic method.**  It is far
+  more accurate than ordinary random Crofton at the same ray count, avoids the
+  severe orientation bias seen in the legacy grid on thin or oblique models,
+  supports uncertainty intervals, and reached a 1% volume-interval target on
+  every bounded `share/db` case in the 10-second study.
+- **Rotated grid is the strongest deterministic control.**  It converges very
+  quickly on the analytic fixtures and usually agrees closely with Crofton on
+  complex-model volume, centroid, and inertia.  It is not uniformly reliable:
+  the full Generic Twin surface-area result remains about 6.5% above the two
+  Crofton results.
+- **Legacy axis-aligned GQA remains useful for compatibility, not truth.**  It
+  can be excellent on aligned or symmetric geometry, but agreement among its
+  three axial views does not detect shared grid bias.  Havoc volume is 18.2%
+  above the new-method consensus at the finest practical legacy grid tested.
+- **Ordinary random Crofton is primarily a control.**  Its observed error falls
+  close to the expected inverse-square-root rate, but QMC needs substantially
+  fewer rays for the same precision on these tests.
+
+The recommended experimental default is therefore time-first QMC with an
+explicit accuracy status, retaining rotated grid and random Crofton as
+diagnostic controls.  Legacy output should remain available as a compatibility
+and comparison control.  The experimental path now supplies conservative
+bounds, traditional Plot3 files, bounded representative plots, and incremental
+MGED visualization.
+
+## Test protocol
+
+- Analytic fixtures: a 1000 mm cube and radius-1000 mm sphere.
+- Solid density: synthetic 1 g/cm^3.  Consequently mass relative error equals
+  volume relative error.  The air fixture uses 0.001 g/cm^3 for modeled air.
+- Random and QMC convergence values are medians across ten fixed outer seeds.
+- Comparable-effort analytic table: 76,800 legacy rays, about 60,000 rotated
+  rays, and 100,000 random/QMC rays.
+- Complex-model comparison: the finest practical legacy grid, the deepest
+  completed rotated grid, and fixed high-ray Crofton runs.  Those are
+  stabilization comparisons, not strictly equal-work comparisons.
+- The four complex models are Havoc, the M35 cab, wrapped `Goliath.c`, and the
+  complete Generic Twin `all` hierarchy.
+- Complex models have no analytic truth.  “Consensus” below is the arithmetic
+  mean of the three new-method values and is only a comparison reference.
+
+For uncertainty, the implementation uses independent randomized replications,
+following the design guidance in A. B. Owen, [“Error estimation for quasi-Monte
+Carlo” (2025), section 5](https://arxiv.org/abs/2501.00150v3).  The pilot
+forecast and its safety factor are engineering heuristics, not results
+attributed to Owen.  Only fresh production replications determine the reported
+interval.
+
+## Analytic measurement accuracy
+
+All entries are errors; lower is better.  Experimental inertia tensors were
+translated from the model origin to the object centroid before comparison with
+legacy output.
+
+| Output | Fixture | Legacy axis grid | Rotated grid | Crofton random | Crofton QMC | Result |
+|---|---|---:|---:|---:|---:|---|
+| Volume | Cube | **0%** | 0.00392% | 0.18574% | 0.00674% | Legacy wins this specially aligned case |
+| Volume | Sphere | 0.01170% | 0.00679% | 0.24152% | **0.00211%** | QMC improves on legacy by 82.0% |
+| Mass | Cube | **0%** | 0.00392% | 0.18574% | 0.00674% | Same relative result as volume |
+| Mass | Sphere | 0.01170% | 0.00679% | 0.24152% | **0.00211%** | Same relative result as volume |
+| Surface area | Cube | unavailable | **0.00468%** | 0.12256% | 0.03220% | New capability; rotated wins |
+| Surface area | Sphere | unavailable | 0.01267% | 0.39150% | **0.00200%** | New capability; QMC wins |
+| Centroid | Cube | 3.6079 mm | **0.0102 mm** | 1.1650 mm | 0.1952 mm | Rotated improves on legacy by 99.7% |
+| Centroid | Sphere | **approximately 0** | 0.0261 mm | 3.0817 mm | 0.1139 mm | Axial symmetry gives legacy exact cancellation |
+| Inertia | Cube | 0.00276% | **0.00224%** | 0.32629% | 0.02884% | Rotated improves on legacy by 18.8% |
+| Inertia | Sphere | 0.01534% | **0.01222%** | 0.42519% | 0.01539% | QMC is effectively tied; rotated wins |
+
+Selected relative changes make the tradeoff clearer:
+
+| Comparison at comparable effort | Change relative to legacy error |
+|---|---:|
+| Sphere volume, rotated | 41.9% lower |
+| Sphere volume, QMC | 82.0% lower |
+| Sphere volume, random | 1,964% higher |
+| Cube centroid, rotated | 99.7% lower |
+| Cube centroid, QMC | 94.6% lower |
+| Cube inertia, QMC | 944% higher |
+| Cube inertia, random | 11,713% higher |
+| Sphere inertia, rotated | 20.4% lower |
+
+The apparent contradictions are real sampling effects.  An aligned grid can
+integrate an aligned box volume exactly while placing its estimated centroid
+at cell centers.  Symmetric axial views can cancel sphere centroid error.  That
+does not make axis-aligned sampling uniformly more accurate on arbitrary
+geometry.
+
+## Convergence rates
+
+The fitted relation is `error ~ rays^slope`; a more negative slope indicates
+faster observed error reduction.  Fits span the measured sampling ladders and
+should not be extrapolated as guarantees.
+
+| Method | Volume slope | Area slope | Centroid slope | Inertia slope |
+|---|---:|---:|---:|---:|
+| Legacy axis grid | -1.025 on sphere; cube exact | unavailable | -0.500 on cube; sphere cancels | -1.000 to -1.037 |
+| Rotated grid | -0.650 to -0.918 | -0.673 to -0.681 | -0.934 to -1.118 | -0.863 to -0.898 |
+| Crofton random | -0.460 to -0.462 | -0.480 to -0.487 | -0.487 to -0.531 | -0.480 to -0.546 |
+| Crofton QMC | -0.576 to -0.696 | -0.573 to -0.753 | -0.786 to -0.844 | -0.757 to -0.834 |
+
+Random Crofton behaves close to inverse square root.  QMC improves every
+measured random-Crofton slope.  Grid convergence is faster on these smooth
+fixtures, but its result can plateau around an orientation-specific bias.
+
+The next table reports the first tested sampling rung meeting each threshold
+on **both** analytic fixtures.  A dash means the threshold was not reached by
+the largest tested rung.
+
+| Method | Volume 1% / 0.1% / 0.01% | Area 1% / 0.1% / 0.01% | Centroid 1 mm / 0.1 mm | Inertia 1% / 0.1% / 0.01% |
+|---|---:|---:|---:|---:|
+| Legacy axis grid | 1,200 / 19,200 / 307,200 | — | 4,915,200 / — | 1,200 / 19,200 / 307,200 |
+| Rotated grid | 2,940 / 14,373 / 239,005 | 2,940 / 59,501 / 954,145 | 2,940 / 59,501 | 2,940 / 14,373 / 239,005 |
+| Crofton random | 25,000 / 1,600,000 / — | 25,000 / 800,000 / — | 800,000 / — | 25,000 / 3,200,000 / — |
+| Crofton QMC | 12,500 / 12,500 / 100,000 | 12,500 / 50,000 / 1,600,000 | 25,000 / 400,000 | 12,500 / 25,000 / 800,000 |
+
+## Time-first and adaptive behavior on `share/db`
+
+### Ten seconds or a 1% interval
+
+The implemented pilot/forecast/production planner was run on all 138 bounded,
+non-half-containing roots from the original corpus.  The requested result was
+a pointwise 95% volume-interval half-width no larger than 1%, with a 10-second
+deadline and no ray cap.
+
+| Outcome | Result |
+|---|---:|
+| Target met | **138 / 138** |
+| Median rays | 300,000 |
+| 90th-percentile rays | 1,863,584 |
+| Maximum rays | 45,943,072 |
+| Median time | 0.11 s |
+| 90th-percentile time | 1.41 s |
+| 95th-percentile time | 2.53 s |
+| Maximum time | 8.40 s |
+| Median relative half-width | 0.368% |
+| 90th-percentile relative half-width | 0.652% |
+| Maximum relative half-width | 0.896% |
+
+This replaces the misleading earlier 100,000-total-ray result, where only
+49/138 met 1% because the cap also paid for a discarded pilot and often left
+only 25,000–75,000 production rays.  Ray count varies by more than two orders
+of magnitude across the corpus, while all wall times remained under the
+declared deadline.  That strongly supports time as the primary user-facing
+resource control.
+
+### Automatic dimensional stability with a 60-second default
+
+The dimensional stability rule used a 0.00005 mm target, scaled downward for
+very small objects, and required two stable doubling comparisons.
+
+| Outcome | Result |
+|---|---:|
+| Stopped stable before deadline | 23 / 138 |
+| Reached 60-second limit | 115 / 138 |
+| Mean rays | 104.2 million |
+| Within the independent adaptive result's 95% interval | 129 / 138 |
+| Median difference from adaptive point estimate | 0.100% |
+| 90th-percentile difference | 0.334% |
+| 95th-percentile difference | 0.454% |
+| Maximum difference | 1.204% |
+
+The independent adaptive interval is not analytic truth, so the interval count
+is an agreement check rather than a coverage claim.  The practical conclusion
+is still clear: the modeling-distance stability floor is much stricter than a
+1% answer requirement for most corpus objects.  It is valuable as an absolute
+terminal criterion, but using it alone makes the default consume the full
+minute in 83.3% of cases.
+
+### Legacy comparison over the corpus
+
+At the earlier 100,000-ray cap, 123 roots produced paired legacy and Crofton
+volumes.  Median relative difference was 0.544%, the 90th percentile was
+3.573%, and 113/123 were within 5%.  Fifteen legacy runs exceeded 120 seconds;
+all had fast Crofton answers.  Analytic checks on notable disagreements favor
+the new estimator:
+
+| Model | Analytic volume (mm^3) | Legacy error | Crofton error |
+|---|---:|---:|---:|
+| `pic.g/box.r` | 1.17368243 | -10.126% | -0.175% |
+| `cube.g/light.r` | 6.08899047 | +1.940% | -0.0146% |
+| `lgt-test.g/box.r` | 19,233.2255 | +2.496% | -0.175% |
+
+## Large complex models
+
+Volumes are in billions of mm^3.  Legacy uses its finest practical tested grid;
+Havoc, M35, and Goliath use rotated level 3 and 3.2 million Crofton rays.  The
+full Generic Twin uses rotated level 5 and 12.8 million Crofton rays.
+
+| Model | Legacy | Rotated | Random | QMC | Legacy vs new-method consensus |
+|---|---:|---:|---:|---:|---:|
+| Havoc | 2.864720 | 2.420378 | 2.422244 | 2.426672 | **+18.226%** |
+| M35 cab | 0.547483 | 0.537888 | 0.538267 | 0.537940 | +1.757% |
+| Goliath | 0.708337 | 0.709737 | 0.709080 | 0.709610 | -0.161% |
+| Generic Twin, full | 29.415300 | 29.725600 | 29.830135 | 29.801455 | -1.244% |
+
+The legacy finest-grid centroids differ from the new-method centroid mean by
+about 547 mm on Havoc, 11.1 mm on M35, 2.1 mm on Goliath, and 8.1 mm on the full
+Generic Twin.  Havoc is the clearest evidence of axis-aligned aliasing in a
+real model.
+
+Cross-method spread among the three new methods summarizes consistency for all
+numeric outputs.  Volume and mass have the same relative spread because of the
+uniform synthetic density.  Inertia spread uses tensor trace about the model
+origin; centroid reports the maximum pairwise distance.
+
+| Model | Volume/mass spread | Area spread | Centroid separation | Inertia-trace spread |
+|---|---:|---:|---:|---:|
+| Havoc | 0.260% | 0.336% | 12.31 mm | 0.203% |
+| M35 cab | 0.070% | 0.113% | 1.73 mm | 0.107% |
+| Goliath | 0.093% | 0.0097% | 0.68 mm | 0.204% |
+| Generic Twin, full | 0.351% | **6.445%** | 23.75 mm | 0.344% |
+
+The Generic Twin area disagreement is the main unresolved answer regression.
+Random and QMC agree with each other to 0.108%, while rotated grid is 6.5%
+higher.  This points to orientation/boundary-count behavior in the rotated
+area estimator rather than Crofton noise.
+
+The declared 10-second, 1%-interval mode also succeeded on all four complex
+models after the half-space correction:
+
+| Model | Volume (10^9 mm^3) | Relative 95% half-width | Rays | Internal time | Difference from deep new-method consensus |
+|---|---:|---:|---:|---:|---:|
+| Havoc | 2.430949 | 0.600% | 3,119,616 | 1.691 s | +0.324% |
+| M35 cab | 0.536708 | 0.352% | 1,360,416 | 0.407 s | -0.246% |
+| Goliath | 0.710126 | 0.221% | 300,000 | 1.221 s | +0.092% |
+| Generic Twin, full | 29.871938 | 0.463% | 3,746,112 | 1.526 s | +0.289% |
+
+Goliath illustrates why time and rays are separate concerns: only 300,000 rays
+were needed statistically, but its geometry makes each ray much more expensive
+than in the other models.
+
+## Modeled-air consistency
+
+Expected volume is 2.0 billion mm^3 with air included and 1.0 billion mm^3
+without it.  Expected included mass is 1,001,000 g and excluded mass is
+1,000,000 g.
+
+| Output | Legacy | Rotated | Random | QMC |
+|---|---:|---:|---:|---:|
+| Included volume error | **0%** | +0.00339% | +0.39436% | +0.00286% |
+| Excluded volume error | **0%** | -0.00289% | +0.02053% | +0.00081% |
+| Included mass error | **0%** | -0.00288% | +0.02127% | +0.00081% |
+| Included centroid error | 14.43 mm | **0.052 mm** | 1.01 mm | 0.184 mm |
+| Excluded centroid error | 14.43 mm | **0.052 mm** | 1.01 mm | 0.184 mm |
+
+All methods apply air inclusion consistently.  The legacy exact totals again
+coexist with a cell-center centroid bias.
+
+## Diagnostic outputs
+
+Diagnostics are sampled observations, not proofs that unobserved problems are
+absent.  “Largest” values are therefore directional extreme samples and
+generally increase as more favorable rays are found.
+
+| Check | Legacy | Rotated | Random | QMC | Assessment |
+|---|---|---|---|---|---|
+| Overlap | Found; max 8,000 mm | Found; 2,888.6 mm | Found; 9,714.2 mm | Found; 10,185.5 mm | Analytic directional maximum is 11,357.8 mm; QMC is closest |
+| Gaps/internal voids | Found | Found | Found | Found | New output distinguishes `internal-void` and both region orders |
+| Adjacent air | Legacy stand-alone run crashed without a plot prefix | Found at 0 mm | Found at 0 mm | Found at 0 mm | Existing legacy defect; experimental modes are clean |
+| Exposed air | Found; max 9,000 mm | Found; 12,122.8 mm | Found; 12,765.4 mm | Found; 12,406.1 mm | Analytic box diagonal is 14,456 mm; random is closest at this budget |
+
+Crofton has much better angular coverage than three fixed views for overlap and
+exposed-air extremes.  Rotating only one triad does not guarantee that it will
+sample the direction of a narrow or maximal defect.
+
+## Output and interface coverage
+
+| Capability | Legacy | Rotated grid | Crofton random/QMC | Status |
+|---|---|---|---|---|
+| Volume | yes | yes | yes | Comparable |
+| Mass/weight | yes | yes | yes | Comparable with density data |
+| Centroid | yes | yes | yes | Comparable |
+| Moments/products of inertia | yes, about centroid | yes, about model origin | yes, about model origin | Convention mismatch must remain explicit |
+| Surface area | no | yes | yes | Experimental improvement |
+| Overlaps | yes | yes | yes | Crofton negative result remains probabilistic |
+| Gaps/internal voids | yes | yes | yes | Experimental labeling is clearer |
+| Adjacent air | yes, but tested crash | yes | yes | Legacy defect exposed |
+| Exposed air | yes | yes | yes | Better angular coverage experimentally |
+| Modeled-air include/exclude | yes | yes | yes | Consistent |
+| Per-object/per-region results | text | JSON/text | JSON/text | Experimental output is easier to consume |
+| Bounding-box report | yes | yes | yes | Experimental bounds need no sampling rays and reject unbounded results |
+| Plot3 diagnostic output | yes | all or representative | all or representative | Experimental modes also support incremental MGED display |
+| Progress, replay, uncertainty | no | progress | progress, invalid-ray replay, intervals | Experimental improvement |
+
+## Half-space behavior
+
+Half-spaces must remain in ray tracing when they clip finite geometry; simply
+removing them changes Boolean meaning.  The corrected implementation derives
+conservative sampling bounds from prepared region trees and ignores infinite
+solids only when constructing the Crofton sampling sphere.
+
+- A cube intersected by a half-space and the complementary subtraction each
+  report the expected 500,000,000 mm^3.
+- The combined split reports the original 1,000,000,000 mm^3 in both grid and
+  Crofton tests.
+- The complete Generic Twin, whose fuel regions use half-space clipping, now
+  runs successfully and meets the adaptive 1% target in 1.526 seconds.
+- `cube.g/cloud.r.bak` is another bounded half-containing root and now runs.
+- A standalone half-space, or a selected union containing an independently
+  unbounded region, is rejected with an explicit unbounded-region diagnostic.
+- Unboundedness is now propagated through `SUBMODEL`.  This prevents the prior
+  `primitives.g/other` infinite-hit crash while preserving bounded submodels
+  whose internal half-space is clipped by finite geometry.
+
+Of the twelve original direct or indirect half-containing corpus roots, two
+are bounded and now analyzable; ten are genuinely unbounded selections and are
+rejected before sampling.
+
+## Recommendations
+
+1. **Use QMC for the experimental accuracy-oriented default.**  Keep a time
+   ceiling as the primary resource control, report progress, and always expose
+   `target_met`, `target_not_met`, or `time_limit` with the achieved interval.
+2. **Retain the geometry-scale stability rule as an absolute floor, not the
+   sole ordinary stopping target.**  It consumed the full minute on 115/138
+   cases even though the independent comparison was already subpercent for
+   almost all of them.
+3. **Keep rotated grid as a deterministic control.**  It is excellent on the
+   analytic fixtures and gives users a non-statistical alternative, but it
+   should not replace QMC until the Generic Twin area discrepancy is resolved.
+4. **Investigate rotated-grid surface-area accumulation next.**  Test multiple
+   independent triads, internal-interface counting, clipped boundary-cell
+   weights, and convergence under orientation changes.  The Crofton pair's
+   0.108% agreement gives a strong reference, though not analytic truth.
+5. **Keep ordinary random Crofton for validation and debugging.**  It is an
+   important unbiased control for QMC, but should not be the performance
+   default at the tested precision levels.
+6. **Keep representative plotting as the interactive default.**  The `all`
+   policy is valuable for exact legacy-style debugging but grows directly with
+   ray count.  Continue using bounded spatial representatives for ordinary
+   live viewing, and clearly document inertia reference conversion when
+   presenting legacy and experimental tensors side by side.
+7. **Expand interval calibration.**  Repeat the 10-second corpus with multiple
+   independent outer seeds and analytic fixtures to measure empirical coverage,
+   especially for area, centroid, inertia, near-zero results, and
+   `--accuracy-scope all`.
+8. **Retain half-space/submodel regressions.**  Test finite intersection,
+   subtraction, bounded submodels, standalone unbounded regions, and nested
+   unbounded submodels.  Do not filter half-spaces out of Boolean evaluation.
+9. **Treat diagnostic negatives conservatively.**  Report sampling evidence
+   and event-probability bounds, never “no overlap exists,” unless a separate
+   exhaustive geometric proof is implemented.
@@ -4,6 +4,7 @@ set(
  api_reference.adoc
  bool_eval_development.adoc
  getting_started_as_contributor.adoc
+  gqa_analysis.adoc
  iterate_geometry.adoc
  libged_command_authoring.adoc
  raytracing_pipeline.adoc
@@ -0,0 +1,244 @@
+= GQA Analysis Design
+
+The experimental `gqa --analyze` path separates geometric observations from
+reporting.  Ray callbacks accumulate region-attributed raw moments, boundary
+crossings, and diagnostic events; the command layer applies units, object
+rollups, and schema-versioned output.  The engine is exposed through the
+C-compatible `analyze_gqa` libanalyze API; libged supplies only its open
+database, result string, and density-loader adapter.  The reusable sampling
+interval calculation is separately exposed as
+`analyze_sampling_interval95`; paired ratio estimates use
+`analyze_sampling_ratio_interval95`.
+
+== Measurement conventions
+
+Lengths use millimetres internally, mass uses grams, surface area counts region
+partition boundaries, and inertia is about the model origin.  Modeled air is
+included by default and can be excluded explicitly.  Grid boundary cells are
+clipped to the model bounds.  These choices are written to JSON so a consumer
+does not have to infer them from command options.
+
+The optional `bounds` measurement is not sampled.  Conservative bounds are
+combined from the prepared region Boolean trees rather than directly unioning
+raw primitive boxes.  This keeps finite intersections and subtractions
+involving half-spaces usable while rejecting a result whose region-tree bound
+is unbounded.  As elsewhere in librt, subtraction does not shrink this
+conservative box.  The same bounds define the sampling domain and the
+automatic plot scale.
+
+== Grid sampling
+
+The experimental grid implementation has two modes.  The `grid` mode is an
+axis-aligned three-view control.  The `grid-rotated` mode constructs an
+orthonormal direction triad from an azimuth and elevation, projects all eight
+corners of the tight model bounding box into each perpendicular view plane,
+and traces the resulting rectangular grids.  Both modes clip boundary cells
+to their projection intervals and weight each chord by its cell's actual
+area.  Rotated-cell covariance is included when accumulating mass moments;
+otherwise an oblique cell would incorrectly be treated as having only
+axis-aligned transverse extent.
+
+The grid worker also records the clipped width and orientation of every
+sampled cell.  These are visualization metadata only: volume and moment
+accumulation continues to use the same clipped area and covariance.
+
+The rotated mode changes the two-dimensional grid phase between refinements.
+The general shifted-lattice technique is described by R. Cranley and T. N. L.
+Patterson, _Randomization of Number Theoretic Methods for Multiple
+Integration_ (1976), https://doi.org/10.1137/0713071.  GQA uses deterministic
+irrational shifts for reproducibility, so they are an anti-aliasing measure,
+not the randomized construction or an uncertainty estimate from that paper.
+
+Grid rays execute in parallel using one librt resource and one GQA accumulator
+per worker.  Region metrics, problem observations, and material errors are
+merged only after a batch completes.  This preserves the existing serial
+bookkeeping semantics without locking every partition callback.  A refinement
+is publishable only after all three views finish; rays from an incomplete pass
+at the deadline are counted and reported as discarded rather than mixed into
+the preceding complete estimate.
+
+Grid convergence uses dimensional quantities.  The change between successive
+area and volume estimates is expressed as the corresponding sphere-radius
+change, and the maximum-minus-minimum directional volume radius measures
+view-dependent aliasing.  Both must be below the requested distance for two
+successive refinements.  Stability across levels alone is insufficient:
+axis-aligned geometry can produce a nearly unchanged average while one view
+remains anomalous.  When a half-spacing pass will not fit the remaining time,
+the planner may choose an intermediate spacing that is still at least 1.05
+times finer; it stops instead of starting a refinement too small to add useful
+resolution.
+
+== Crofton sampling
+
+An oriented line needs four uniform coordinates: two for direction and two for
+the perpendicular offset.  The default QMC implementation uses a randomized
+rank-one lattice when OpenQMC is available; ordinary pseudorandom sampling is
+also available as a diagnostic cross-check.  A ray offset continues one point
+set.  A seed and stream identifier select a separately randomized point set.
+The distinction matters because consecutive points from one QMC sequence are
+dependent and cannot be treated as independent error replicates.
+
+The implementation follows Art Owen's _Error estimation for quasi-Monte Carlo_
+(2025), especially section 5:
+https://arxiv.org/abs/2501.00150v3.  Owen recommends independent randomizations,
+at least ten replications, and a Student-t interval across replicate estimates.
+GQA uses 16 replications.  Libanalyze uses the conservative 95 percent
+Student-t critical value for nine degrees of freedom for any accepted replicate
+count, avoiding intervals that are too narrow because a normal critical value
+was used for a small replication count.
+
+The reported intervals are approximate and pointwise.  Selecting
+`--accuracy-scope all` requires each selected model, unambiguous object, and
+region metric to pass, but does not claim simultaneous 95 percent family-wise
+coverage.
+
+== Adaptive accuracy
+
+Adaptive accuracy has two statistically separate stages:
+
+. A pilot set of 16 randomized replications estimates the production budget.
+  The normal deadline-driven pilot uses 100000 rays.  A root-n squared-error
+  scaling and an empirically calibrated safety factor are used only as planning
+  heuristics.
+. The pilot estimates are discarded.  Sixteen fresh streams produce the
+  published estimates, intervals, and pass/fail result.
+
+Fresh validation prevents the stopping decision and the reported interval from
+using the same favorable random fluctuation.  QMC convergence need not follow a
+stable power law, so the production interval is authoritative: an optimistic
+pilot can reduce work, but it cannot assert success.  This separation, the use
+of independently randomized replications, and the Student-t interval follow
+A. B. Owen, _Error estimation for quasi-Monte Carlo_ (2025), section 5:
+https://arxiv.org/abs/2501.00150v3.  The root-n forecast and safety factor are
+engineering policy, not claims attributed to Owen.
+
+When an adaptive target is combined with `--time` and no explicit `--rays`,
+the pilot also measures throughput.  The forecast is clamped to 80 percent of
+the ray count expected to fit in the remaining wall time, reserving headroom
+for variation and reporting.  An explicit ray count remains an additional hard
+cap.  A deadline-constrained production estimate that misses its target is
+reported with `time_limit`; it is not promoted to success.
+
+The 16 streams execute serially through one reusable Crofton worker session.
+Each stream's rays use the session's all-CPU worker pool.  Worker-local events
+are merged and delivered to libanalyze callbacks serially, preserving GQA's
+per-region bookkeeping without launching competing all-CPU replications.
+
+== Default resource and stability policy
+
+Grid and Crofton measurement-only runs use elapsed time as their default
+resource limit.  Unless the caller supplies the sampler's explicit work limit
+(`--refine` or `--rays`) or `--time`, the analysis deadline is 60000 ms.  An
+explicit ray count remains available as a Crofton hard cap for reproducibility
+and specialized testing, but ray count is not the ordinary user-facing
+resource control.
+
+For non-uncertainty measurement runs without an explicit stability target,
+the automatic equivalent-radius target is
+
+`min(0.1 * BN_TOL_DIST, 0.001 * smallest_positive_tight_model_span)`.
+
+The first term asks for stability an order of magnitude below librt's normal
+0.0005 mm modeling tolerance.  The scale-relative term begins tightening the
+target for models below 0.05 mm, before the fixed tolerance becomes large
+relative to the selected geometry.  Finite region-tree bounds are used because
+the prepared model RPP is deliberately rounded outward and raw half-space
+bounds are infinite.  This stopping rule is an
+engineering convergence heuristic and is reported as `stable`; it does not
+claim the pointwise coverage of the independent-replication intervals.
+
+For Crofton sampling, only estimates needed by the selected measurements
+participate in the test: area requests stabilize the surface-area radius,
+while volume, mass, centroid, and moment requests stabilize the volume radius.
+Stability is evaluated at 20,000 rays and then whenever the cumulative ray
+count doubles.  Two successive doubling intervals must meet the target.
+Comparing fixed 2,000-ray increments would be misleading because each
+increment becomes a vanishing fraction of the cumulative estimate and its
+apparent change eventually tends to zero even when sampling error remains
+material.
+
+Grid sampling applies the selected area or volume-radius refinement change and
+always checks directional volume-radius spread as an aliasing diagnostic.
+
+The Crofton progress callback runs serially after completed parallel ray
+batches.  GQA throttles those callbacks to approximately one report per second
+and includes elapsed time, total rays, throughput, the current volume estimate,
+and observed versus requested stability.  A callback never reads worker-local
+state while rays are in flight.
+
+== Plot publication and live views
+
+Plotting is downstream of sampling.  Grid workers append immutable sample
+records to worker-local accumulators, and Crofton visitors append records
+while librt is in its serialized delivery phase.  Libanalyze publishes only
+after a parallel batch completes.  Consequently Plot3 output and live-view
+output consume the same ordered line batches and cannot alter ray selection,
+partition bookkeeping, or numerical estimates.
+
+The `all` policy preserves every chord and diagnostic line, including rays
+from successive complete grid refinements.  It is useful for exact debugging
+but its storage and display cost grows with sample count.  The default
+`representative` policy quantizes sample centers in model space and retains a
+deterministic one-element min-hash reservoir for each category, direction
+class, region pair, and occupied cell.  The reservoir construction follows
+the mergeable fixed-size sampling idea in J. S. Vitter, _Random Sampling with
+a Reservoir_ (1985), https://doi.org/10.1145/3147.3165.  Severity wins first
+for diagnostic observations so the most significant issue in a cell remains
+visible; deterministic priority resolves equal-severity samples.
+
+Representative samples are rendered as finite-support glyphs.  Grid glyphs
+use the exact clipped pixel footprint.  Crofton glyphs use the representative
+side length corresponding to their average `pi R^2/N` transverse sampling
+area, where `R` is the bounding-sphere radius and `N` is the ray count.  The
+spatial-cell scale independently limits how many glyphs are retained.  This
+follows the visual principle of surfels, which draw
+sample support instead of infinitesimal points: H. Pfister et al., _Surfels:
+Surface Elements as Rendering Primitives_ (2000),
+https://doi.org/10.1145/344779.344936.  Bounding the number of displayed
+samples by occupied spatial cells applies the fixed-detail motivation of S.
+Rusinkiewicz and M. Levoy, _QSplat: A Multiresolution Point Rendering System
+for Large Meshes_ (2000), https://doi.org/10.1145/344779.344940, without
+claiming to implement its hierarchy.
+
+An incomplete grid refinement never replaces the last complete
+representative snapshot.  Adjacent-air samples use an interface-centered
+footprint and a bounded directional tick; the full following air chord would
+overstate the location of the transition.
+
+The libanalyze callback contains plain immutable model-coordinate lines and
+has no display-system dependency.  MGED copies each batch, synchronously
+marshals the update through its existing GUI-event bridge, rebuilds the named
+`gqa::analysis` view object on the display-owning thread, and requests a
+refresh.  This is the same ownership rule used for other asynchronous GED
+work: worker threads calculate data, while all view-list mutation remains on
+the GUI thread.
+
+Relative and absolute tolerances may be combined.  For a scalar estimate
+`theta`, the accepted half-width is the larger of
+`relative_error * abs(theta)` and the applicable absolute tolerance.  This
+lets absolute tolerances govern quantities near zero.  A centroid point
+estimate is the ratio of the replicate-mean first moment to replicate-mean
+mass.  Its interval uses paired first-moment and mass residuals in a
+first-order delta calculation, retaining their covariance; zero-mass
+centroids remain unavailable.
+
+== Diagnostic evidence
+
+Overlap, gap, and modeled-air results are sampled observations, not proofs of
+absence.  Crofton overlap-only runs first find positive-volume intersections of
+region-instance bounding boxes and divide rays evenly among those candidates.
+Raw overlap events are clipped to their candidate box and filtered by the
+requested depth tolerance.  For an ordinary randomized sequence with zero
+events, the report includes the exact one-sided 95 percent binomial upper bound
+`1 - 0.05^(1/N)` on the per-ray event probability.
+
+Invalid partitions are rejected deliberately and may be written to a separate
+replay file.  Publication uses temporary sibling files followed by rename so a
+failed write does not leave a new partial report at the requested path.
+
+== Reproducibility and schema
+
+The seed, resolved sequence, total pilot and production ray counts, unused
+budget, measurement conventions, and every interval are recorded in JSON
+schema version 2.  Each diagnostic observation retains its actual ray ID,
+origin, and direction.  Invalid-ray files can be replayed by `nirt --replay`.
@@ -37,6 +37,179 @@ For error and interface calculations, the grid is refined until an error is foun
 When multiple calculations are performed, the most extensive path is taken. For example, if weight and overlap calculations are performed, then the grid will be refined until the refinement limit is reached, or an overlap is detected.


+[[experimental_analysis]]
+== EXPERIMENTAL ANALYSIS
+
+The separate `--analyze` interface preserves the behavior of existing *gqa*
+invocations.  It reports a concise summary and can write all region results and
+individual problem observations as structured JSON.  The JSON includes
+per-object rollups when every region belongs unambiguously to one requested
+object.  This interface requires an explicit object list, including when used
+inside MGED.
+
+[source]
+----
+gqa --analyze [analysis options] model.g object [object ...]
+----
+
+`--measure` accepts a comma-separated selection from
+`volume,mass,area,centroid,moments,bounds`.  With no measurement or check selection,
+the default is `volume,area`.  A check-only selection does not require density
+data.  Centroid and moments require mass.  `--check` accepts
+`overlaps,gaps,adjacent-air,exposed-air` and defaults to no checks.
+The model totals and region measurements use mm, mm^3^, mm^2^, and grams.
+Mass needs a density table, selected with `--density file` or found using the
+normal gqa density search.
+The `bounds` measure reports the finite, conservative region-tree bounding box
+and its three face areas without firing sampling rays.  It can be combined
+with the sampled measures.  Subtraction does not shrink the reported box.  A
+bounded intersection or subtraction involving a half-space is supported; a
+selected region whose region-tree bound is unbounded is rejected.
+
+`--sampler grid`, the experimental interface's default, is the axis-aligned
+control corresponding to the traditional three-view method.  `--sampler
+grid-rotated` uses a non-axis-aligned orthonormal view triad.  Its first
+direction defaults to 35 degrees azimuth and 25 degrees elevation and may be
+set with `--azimuth degrees` and
+`--elevation degrees`; the other two directions are derived orthogonally.
+Rotated refinements use deterministic plane-grid shifts so successive passes
+do not repeatedly sample the same alignment.  These shifts reduce aliasing but
+do not make grid results statistical confidence estimates.
+
+Both grid modes use `--spacing mm`.  When it is omitted, the initial spacing
+is the larger of 50 mm and one 128th of the tight model bounding-box diagonal.
+This keeps the first pass tractable for exceptionally large models while
+preserving the familiar 50 mm start for ordinary models.  `--refine levels`
+sets the maximum number of additional grid passes and disables the default
+time control when used alone.  Each normal refinement halves the spacing;
+under a time limit, a final intermediate spacing may be selected so a complete
+three-view pass fits the remaining budget.  Boundary cells are clipped and
+weighted by their actual projected area rather than counted as full cells.
+
+`--sampler crofton` uses isotropic rays.  For a measurement-only run with
+either grid or Crofton sampling and no explicit resource control, sampling has
+a 60000 ms budget and stops sooner when its automatic equivalent-radius
+stability target is reached.  The planner may stop before the exact deadline
+when there is not enough time for another complete sampling checkpoint.
+For Crofton sampling, `--rays count` replaces that time-first default with an
+explicit hard ray cap; `--sequence random|qmc` chooses an explicit sequence.
+If omitted, the standard librt sequence selection is used.  A QMC request
+requires a build with OpenQMC.
+`--seed value` makes the selected randomized sequence reproducible.
+`--stability mm` requests equivalent-radius convergence.  For a grid, both the
+change from the preceding refinement and the spread among its three
+directional volume estimates must meet the target for two successive passes.
+For Crofton sampling, two successive doubling checkpoints must meet the
+target.  `--time ms` sets an internal wall-clock budget.  A time-controlled run
+reports progress at approximately one-second intervals.  Grid progress
+includes the level, spacing, completed rays, throughput, and directional
+volumes; Crofton progress includes the current volume and stability when
+available.
+
+The automatic stability target for either measurement sampler is the smaller
+of `BN_TOL_DIST / 10` (0.00005 mm) and 0.1 percent of the smallest positive
+tight model-bounding-box span.  The scale-dependent term prevents the fixed
+modeling tolerance from becoming coarse relative to a small valid object.
+This is a convergence heuristic, not a confidence interval.  It applies only
+to estimates needed by the requested measurements and requires two stable
+comparisons at doubling ray-count checkpoints for Crofton sampling.  Grid
+refinement applies the same two-interval policy to complete three-view passes
+and also requires its directional spread to meet the target.  The result
+distinguishes `stable`, `time_limit`, and an explicitly requested ray or
+refinement limit stopping condition.
+`--tolerance mm` ignores overlap and gap depths at or below the threshold
+(default 0).  Crofton diagnostic checks use the observed ray partition paths
+and therefore remain sampling evidence rather than exhaustive proofs.  Gaps
+between two partitions of the same region are labeled `internal-void` in the
+summary.
+When Crofton is used for an overlap-only request (that is, `--check overlaps`
+with no `--measure`), the experimental path first indexes finite region-instance
+bounding boxes and samples only their intersecting candidate boxes.  The summary
+and JSON record report the number of candidate boxes considered.  This is a
+broad-phase acceleration and does not establish that two solids actually
+intersect; the targeted rays still provide the overlap evidence.
+
+For targeted overlap checks, `--stability mm` compares successive equivalent-radius
+estimates of each candidate box, using the same rays that test the geometry.  The
+known box volume and surface area provide calibration checks, and the sampled rays
+must cover all three dominant direction classes.  A candidate stops only when the
+radius change and these calibration checks meet the requested target, subject to
+the ray budget.  JSON records the calibration errors and whether the target was
+reached in `box_stability_reached` and `box_calibration_reached`.
+
+For randomized sequences, a candidate with no observed overlap event also reports
+a 95 percent upper bound on the per-ray event probability.  This quantifies the
+strength of a sampling-negative result; it is not an exhaustive proof of absence.
+
+The measurements are
+sampled estimates, and diagnostic checks can miss problems that no ray hits.
+
+`--plot-prefix prefix` writes the selected observations as binary Plot3 files
+named `prefixvolume.plot3`, `prefixgaps.plot3`, `prefixoverlaps.plot3`,
+`prefixadj_air.plot3`, and `prefixexp_air.plot3`.  Only files corresponding to
+requested measurements or checks are created.  `--plot-mode all` writes every
+sampled chord, matching the legacy plotting intent and retaining all grid
+refinement passes.  The default, `--plot-mode representative`, keeps one
+deterministic sample per occupied spatial cell, direction class, category, and
+region pair.  It draws the actual projected grid-cell footprint; a Crofton
+glyph represents its average `pi R^2/N` transverse sampling area.  This bounded-detail
+mode preserves visible coverage as refinement proceeds instead of making the
+plot grow in direct proportion to the ray count.  `--plot-resolution mm`
+overrides the automatic representative cell size of one 128th of the model
+diagonal.  Plot selection never changes the rays or numerical estimates.
+
+Inside MGED, `--plot-view` displays the same output in the active view as
+`gqa::analysis` and refreshes it incrementally at approximately one-second
+intervals.  Representative mode replaces the overlay with its latest bounded
+snapshot; all mode appends new rays and can become expensive.  Adjacent-air
+snapshots mark the sampled interface with a footprint and a short directional
+tick, avoiding an invisible zero-length line without implying that the whole
+air chord is defective.  `--plot-view` is unavailable to standalone *gqa*.
+
+`--air include|exclude` controls whether modeled air contributes to
+measurement totals; the default is `include`.  Diagnostic air checks still
+trace modeled air when it is excluded from measurements.  Surface area counts
+region partition boundaries, including internal interfaces, and moments of
+inertia are reported about the model origin.
+
+`--uncertainty` divides the Crofton ray budget among 16 independently
+randomized replications and reports approximate pointwise 95 percent
+Student-t intervals.  It applies only to measurement-only Crofton runs.
+`--relative-error fraction` requests adaptive accuracy and
+`--absolute-error measure=value[,measure=value...]` sets absolute targets for
+`volume` (mm^3^), `mass` (g), `area` (mm^2^), `centroid` (mm), or
+`moments` (g mm^2^).  When both forms apply, the larger tolerance is used.
+`--accuracy-scope model|all` selects model totals only or model, object, and
+region results.  The `all` intervals remain pointwise rather than a
+simultaneous family-wise confidence statement.
+
+Adaptive accuracy uses a pilot stage to forecast a production ray count, then
+discards the pilot estimates and evaluates the target with fresh independent
+production streams.  With `--time ms` and no explicit `--rays`, the deadline
+is the resource limit: the pilot measures ray throughput and the production
+forecast is restricted to fit the remaining time with scheduling headroom.
+When both are specified, the ray count is also a hard cap.  When neither is
+specified, adaptive accuracy uses the default 60000 ms deadline.  A bare
+`--uncertainty` request without an accuracy target retains the 100000-ray
+budget because it has no accuracy condition with which to plan a deadline.
+
+The result explicitly reports `target_met`, `target_not_met`, `time_limit`, or
+`insufficient_evidence`.  `time_limit` means the deadline constrained a
+production result that did not meet the requested interval width; the estimate
+and its achieved interval are still reported.  The pilot, forecast, production
+count and timing, planning safety factor, deadline state, and actual interval
+are recorded in JSON.  The planning heuristic and safety factor do not alter
+the published interval: only the fresh production replications determine it.
+
+`--json file` writes a schema-version 2 JSON record with model and region
+measurements and every observed problem, including full region paths and ray
+information.
+`--invalid-rays file` writes a separate replay-oriented JSON file for rejected
+Crofton partitions.  Each record includes the ray ID, origin, direction, hit
+distances, thickness, region path, and rejection reason.
+Without `--json`, the command prints only the summary.  Reported surface area
+counts region partition boundaries, including internal interfaces.
+
 [[volume_and_weight_calculation]]
 == VOLUME AND WEIGHT CALCULATION

@@ -19,6 +19,20 @@ nirt [options] model.g [objects]...

 *nirt* (Natalie's Interactive Ray-Tracer) traces rays through the specified _objects_ in the database _model.g_, using *librt*(3), and produces an ASCII report of what each ray encounters. Commands are read from standard input. By default the user interacts with *nirt*, repeatedly specifying ray origination points and directions and controlling the format and destination of the reports. Locations may be input and output in either model coordinates (_x_, _y_, _z_) or view (a.k.a. grid-plane) coordinates (_h_, _v_, _d_), and direction may be expressed either as a vector in model coordinates or as azimuth and elevation angles.

+Replay invalid-ray records generated by analysis tools with `--replay` and
+`--ray-id`.  The selected record's origin and direction are restored at JSON
+precision, and one shot is performed before *nirt* exits:
+
+[source]
+----
+nirt --replay invalid_rays.json --ray-id 18432 model.g object
+----
+
+Replay files use schema version 1 or 2 and contain a `rays` array.  Each record
+has a unique nonnegative integer `ray_id`, a finite numeric `origin_mm`
+three-vector, and a finite nonzero `direction` three-vector.  Both options are
+required.  The `--replay=FILE` and `--ray-id=N` forms are also accepted.
+
 When standard input is not a terminal (for example, when it is redirected from a file or a pipe), *nirt* defaults to silent mode, making it convenient for scripted ray queries.

 [[options]]
@@ -38,6 +38,108 @@ region.  Overlaps and other geometry errors invalidate reported volume and
 mass properties.


+[[experimental_analysis]]
+== EXPERIMENTAL ANALYSIS
+
+The separate `--analyze` interface preserves the behavior of existing *gqa*
+invocations.  It reports a concise summary and can write all region results and
+individual problem observations as structured JSON.  The JSON includes
+per-object rollups when every region belongs unambiguously to one requested
+object.  This interface requires an explicit object list, including when used
+inside MGED.
+
+[source]
+----
+gqa --analyze [analysis options] object [object ...]
+----
+
+`--measure` accepts a comma-separated selection from
+`volume,mass,area,centroid,moments,bounds`.  With no measurement or check selection,
+the default is `volume,area`.  A check-only selection does not require density
+data.  Centroid and moments require mass.  `--check` accepts
+`overlaps,gaps,adjacent-air,exposed-air` and defaults to no checks.
+The model totals and region measurements use mm, mm^3^, mm^2^, and grams.
+Mass needs a density table, selected with `--density file` or found using the
+normal gqa density search.
+The `bounds` measure reports the finite, conservative region-tree bounding box
+and its three face areas without firing rays.  Subtraction does not shrink the
+reported box.  Bounded results involving half-spaces are supported; regions
+whose region-tree bounds are unbounded are rejected.
+
+`--sampler grid`, the experimental interface's default, is the axis-aligned
+three-view control.  `--sampler grid-rotated` uses a non-axis-aligned
+orthonormal view triad and deterministically shifts the plane grids between
+refinements to reduce alignment aliasing.  Its first direction defaults to 35
+degrees azimuth and 25 degrees elevation and may be changed with `--azimuth`
+and `--elevation`.
+Both grid modes use `--spacing mm`; when omitted, the initial spacing is the
+larger of 50 mm and one 128th of the tight model diagonal.  `--refine levels`
+sets a maximum number of additional passes.  `--sampler crofton` uses
+isotropic rays.  Measurement-only runs without an explicit resource control
+use a 60000 ms budget and stop sooner when the automatic stability target is
+reached.  The planner may stop before the exact deadline when there is not
+enough time for another complete sampling checkpoint.  `--rays count`
+supplies a Crofton hard ray cap;
+`--sequence random|qmc` chooses an explicit sequence.  A QMC request requires
+a build with OpenQMC.  `--seed value` makes randomized sampling reproducible.
+`--stability mm` requests equivalent-radius convergence.  A grid must meet the
+target for both refinement change and three-view directional spread twice in
+succession.  `--time ms` sets an internal wall-clock budget.
+The automatic target is the smaller of `BN_TOL_DIST / 10` (0.00005 mm)
+and 0.1 percent of the smallest positive tight model span, allowing valid small
+objects to use a proportionally smaller target.  Only estimates needed by the
+requested measurements participate in Crofton convergence, and two
+comparisons at doubling ray-count checkpoints must be stable.  Time-controlled
+runs report progress at approximately one-second intervals.  The summary and
+JSON identify `stable`, `time_limit`, or an explicitly requested ray or
+refinement limit as appropriate.
+`--tolerance mm` ignores overlap and gap depths at or below the threshold
+(default 0).  Gaps between two
+partitions of the same region are labeled `internal-void` in the summary.
+The measurements are
+sampled estimates, and diagnostic checks can miss problems that no ray hits.
+
+`--plot-prefix prefix` writes selected results to the traditional
+`volume.plot3`, `gaps.plot3`, `overlaps.plot3`, `adj_air.plot3`, and
+`exp_air.plot3` suffixes.  `--plot-mode all` records every sampled chord and
+all grid refinements.  The default `representative` mode keeps one
+deterministic sample per occupied spatial cell, direction, category, and
+region pair, and draws finite sample footprints so the plot remains visually
+useful as refinement increases.  `--plot-resolution mm` overrides the default
+cell size of one 128th of the model diagonal.  Plotting does not change the
+sampling calculation.
+
+In MGED, `--plot-view` incrementally refreshes a `gqa::analysis` overlay.  The
+representative overlay is replaced by each bounded snapshot; all mode appends
+rays and can be expensive.  Adjacent-air snapshots show an interface footprint
+and a short directional tick.  Standalone *gqa* does not provide a live view.
+
+`--air include|exclude` controls whether modeled air contributes to
+measurements; inclusion is the default.  Surface area counts region partition
+boundaries and inertia is about the model origin.
+
+`--uncertainty` reports approximate pointwise 95 percent intervals from 16
+independently randomized Crofton replications.
+`--relative-error fraction` and
+`--absolute-error measure=value[,measure=value...]` request adaptive accuracy;
+`--accuracy-scope model|all` controls which reported measurements must meet
+the target.  A pilot estimates the required production budget, but only fresh
+independent production streams determine the reported accuracy status.
+With `--time ms` and no explicit `--rays`, measured throughput limits the
+production plan to the remaining time.  If both are specified, the ray count
+is also a hard cap.  With neither specified, the default deadline is 60000 ms.
+Results distinguish `target_met`, `target_not_met`,
+`time_limit`, and `insufficient_evidence`, and JSON records the pilot,
+forecast, production timing, deadline state, and achieved interval.
+
+`--json file` writes a versioned JSON record with model and region measurements
+and every observed problem, including full region paths and ray information.
+`--invalid-rays file` writes a separate replay-oriented JSON file for rejected
+Crofton partitions.  Each record includes the ray ID, origin, direction, hit
+distances, thickness, region path, and rejection reason.
+Without `--json`, the command prints only the summary.  Reported surface area
+counts region partition boundaries, including internal interfaces.
+
 [[how_gqa_works]]
 == HOW GQA WORKS

@@ -35,14 +35,16 @@
 #include "analyze/debug.h"
 #include "analyze/diff.h"
 #include "analyze/density.h"
+#include "analyze/gqa.h"
 #include "analyze/grid.h"
 #include "analyze/heal.h"
 #include "analyze/info.h"
+#include "analyze/nirt.h"
 #include "analyze/pnts.h"
 #include "analyze/polygonize.h"
-#include "analyze/nirt.h"
-#include "analyze/worker.h"
+#include "analyze/sampling.h"
 #include "analyze/voxelize.h"
+#include "analyze/worker.h"

 #endif /* ANALYZE_H */

@@ -5,12 +5,14 @@ set(
  defines.h
  density.h
  diff.h
+  gqa.h
  grid.h
  heal.h
  info.h
  nirt.h
  pnts.h
  polygonize.h
+  sampling.h
  voxelize.h
  worker.h
 )
@@ -0,0 +1,153 @@
+/*                         G Q A . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @addtogroup libanalyze */
+/** @{ */
+/** @file analyze/gqa.h */
+
+#ifndef ANALYZE_GQA_H
+#define ANALYZE_GQA_H
+
+#include <stddef.h>
+
+#include "common.h"
+#include "analyze/defines.h"
+
+__BEGIN_DECLS
+
+struct analyze_densities;
+struct bu_vls;
+struct db_i;
+
+/**
+ * Load the density table selected by @p filename, or the caller's normal
+ * default density source when it is NULL.  On success, ownership of both
+ * returned allocations passes to analyze_gqa(); @p source must be allocated
+ * with bu_malloc() and @p densities must be releasable by
+ * analyze_densities_destroy().
+ */
+typedef int (*analyze_gqa_density_loader)(
+    struct analyze_densities **densities,
+    char **source,
+    const char *filename,
+    void *data);
+
+/** Receive an informational progress line during a long analysis. */
+typedef void (*analyze_gqa_progress_handler)(const char *message, void *data);
+
+/** Categories used to separate quantitative-analysis plot output. */
+enum analyze_gqa_plot_category {
+    ANALYZE_GQA_PLOT_VOLUME = 0,
+    ANALYZE_GQA_PLOT_GAP,
+    ANALYZE_GQA_PLOT_OVERLAP,
+    ANALYZE_GQA_PLOT_ADJACENT_AIR,
+    ANALYZE_GQA_PLOT_EXPOSED_AIR,
+    ANALYZE_GQA_PLOT_CATEGORY_COUNT
+};
+
+/** Whether a plot batch extends or replaces the preceding visualization. */
+enum analyze_gqa_plot_action {
+    ANALYZE_GQA_PLOT_APPEND = 0,
+    ANALYZE_GQA_PLOT_REPLACE
+};
+
+/** One colored line in model coordinates. */
+struct analyze_gqa_plot_line {
+    int category;
+    unsigned char color[3];
+    double start[3];
+    double end[3];
+};
+
+/**
+ * Immutable plot update.  Storage is valid only for the duration of the
+ * callback.  An empty replacement clears a previous visualization.
+ */
+struct analyze_gqa_plot_batch {
+    int action;
+    const struct analyze_gqa_plot_line *lines;
+    size_t line_count;
+};
+
+/** Receive a serialized plot update after a completed ray batch. */
+typedef void (*analyze_gqa_plot_handler)(
+    const struct analyze_gqa_plot_batch *batch,
+    void *data);
+
+/**
+ * Invocation-local services needed by the quantitative analysis engine.
+ *
+ * The database and result string remain caller-owned.  A density loader is
+ * required only for analyses that request mass-dependent measurements.  The
+ * progress and plot handlers are optional.  Plot handlers run serially after
+ * completed ray batches and must not call back into the active analysis.
+ */
+struct analyze_gqa_context {
+    struct db_i *dbip;
+    struct bu_vls *result;
+    analyze_gqa_density_loader load_densities;
+    void *density_data;
+    analyze_gqa_progress_handler report_progress;
+    void *progress_data;
+    analyze_gqa_plot_handler report_plot;
+    void *plot_data;
+};
+
+/**
+ * Locate the model filename in a standalone @c gqa @c --analyze command line.
+ * This lets a caller open the database before invoking analyze_gqa() without
+ * duplicating the engine's option-arity rules.
+ *
+ * @return the zero-based argument index, or -1 for malformed input.
+ */
+ANALYZE_EXPORT extern int analyze_gqa_model_argument(
+    int argc,
+    const char *argv[]);
+
+/**
+ * Run the versioned GQA analysis interface.
+ *
+ * @p argv uses the command form accepted by @c gqa @c --analyze; argv[0] is
+ * the command name and argv[1] is @c --analyze.  The selected objects follow
+ * all options.  Results and diagnostics replace context->result.
+ *
+ * The function does not retain any context or geometry pointers after return.
+ *
+ * @return ANALYZE_OK on success; ANALYZE_ERROR on invalid input or failure.
+ */
+ANALYZE_EXPORT extern int analyze_gqa(
+    const struct analyze_gqa_context *context,
+    int argc,
+    const char *argv[]);
+
+__END_DECLS
+
+#endif /* ANALYZE_GQA_H */
+
+/** @} */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -0,0 +1,80 @@
+/*                     S A M P L I N G . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+/** @addtogroup libanalyze */
+/** @{ */
+/** @file analyze/sampling.h */
+
+#ifndef ANALYZE_SAMPLING_H
+#define ANALYZE_SAMPLING_H
+
+#include "common.h"
+#include "analyze/defines.h"
+
+__BEGIN_DECLS
+
+/** Minimum independent replications accepted by the interval estimator. */
+#define ANALYZE_SAMPLING_MIN_REPLICATES 10u
+
+struct analyze_sampling_interval {
+    double estimate;   /**< arithmetic mean of the replicate estimates */
+    double half_width; /**< approximate two-sided 95% half-width */
+};
+
+/**
+ * Estimate a pointwise two-sided 95% sampling interval from independent
+ * randomized replicate estimates.  At least ANALYZE_SAMPLING_MIN_REPLICATES
+ * finite values are required.  The implementation uses the conservative
+ * 95% Student-t critical value for nine degrees of freedom for every accepted
+ * replicate count.
+ *
+ * This API is suitable for ordinary Monte Carlo or independently randomized
+ * QMC replications.  Consecutive points from one deterministic or randomized
+ * QMC sequence are not independent replications.
+ *
+ * @return ANALYZE_OK on success; ANALYZE_ERROR for invalid input.
+ */
+ANALYZE_EXPORT extern int analyze_sampling_interval95(
+    struct analyze_sampling_interval *interval,
+    const double *replicate_estimates,
+    size_t replicate_count);
+
+/**
+ * Estimate a pointwise two-sided 95% interval for a ratio of replicate means.
+ * Numerator and denominator values must be paired independent replications.
+ * The point estimate is mean(numerator) / mean(denominator), and the standard
+ * error uses their paired first-order delta-method residuals.  This retains
+ * numerator/denominator covariance and is appropriate for derived quantities
+ * such as a mass centroid.
+ *
+ * @return ANALYZE_OK on success; ANALYZE_ERROR for invalid input or a zero
+ * denominator mean.
+ */
+ANALYZE_EXPORT extern int analyze_sampling_ratio_interval95(
+    struct analyze_sampling_interval *interval,
+    const double *replicate_numerators,
+    const double *replicate_denominators,
+    size_t replicate_count);
+
+__END_DECLS
+
+#endif /* ANALYZE_SAMPLING_H */
+
+/** @} */
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -82,6 +82,10 @@ typedef int (*ged_func_ptr)(struct ged *, int, const char *[]);
 /* Callback related definitions */
 typedef void (*ged_io_func_t)(void *, int);
 typedef void (*ged_refresh_func_t)(void *);
+typedef void (*ged_main_thread_callback_t)(void *);
+/** Run callback synchronously on the display-owning thread. */
+typedef void (*ged_run_on_main_thread_func_t)(void *,
+	ged_main_thread_callback_t, void *);
 typedef void (*ged_create_vlist_solid_func_t)(void *, struct bv_scene_obj *);
 typedef void (*ged_create_vlist_display_list_func_t)(void *, struct display_list *);
 typedef void (*ged_destroy_vlist_func_t)(void *, unsigned int, int);
@@ -235,6 +239,8 @@ struct ged {
    struct ged_callback_state *ged_cbs;
    void (*ged_refresh_handler)(void *);	/**< @brief  function for handling refresh requests */
    void *ged_refresh_clientdata;	/**< @brief  client data passed to refresh handler */
+    ged_run_on_main_thread_func_t ged_run_on_main_thread; /**< @brief run display-affine work on the owning thread */
+    void *ged_main_thread_clientdata; /**< @brief client data passed to the main-thread handler */
    void (*ged_output_handler)(struct ged *, char *);	/**< @brief  function for handling output */
    void (*ged_create_vlist_scene_obj_callback)(void *, struct bv_scene_obj *);	/**< @brief  function to call after creating a vlist to create display list for solid */
    void (*ged_create_vlist_display_list_callback)(void *, struct display_list *);	/**< @brief  function to call after all vlist created that loops through creating display list for each solid  */
@@ -196,12 +196,12 @@ RT_EXPORT extern int rt_obj_prep_serialize(struct soltab *stp, const struct rt_d
 /**
 * Stopping-criteria parameters for rt_crofton_shoot() and rt_crofton_sample().
 *
- * All fields default to zero.  Behaviour when all three are zero (or the
- * pointer is NULL) is identical to the historical default: fire 2 000 rays
- * and stop when two successive 1 %-threshold iterations agree.
+ * All stopping fields default to zero.  Behaviour when all three are zero
+ * (or the pointer is NULL) is identical to the historical default: fire
+ * 2 000 rays and stop when two successive 1 %-threshold iterations agree.
 *
- * When one or more non-zero fields are provided, sampling continues until
- * the FIRST criterion that is met:
+ * When one or more non-zero stopping fields are provided, sampling continues
+ * until the FIRST criterion that is met:
 *
 *   n_rays > 0        Stop once this many rays have been fired in total.
 *
@@ -211,30 +211,65 @@ RT_EXPORT extern int rt_obj_prep_serialize(struct soltab *stp, const struct rt_d
 *                     iterations: r_sa = sqrt(SA / (4*pi)) for surface area
 *                     and r_v = cbrt(3*V / (4*pi)) for volume.  Sampling
 *                     stops when both |Δr_sa| and |Δr_v| (for whichever
- *                     outputs are requested) remain < stability_mm for
+ *                     outputs are requested) remain <= stability_mm for
 *                     consecutive stable windows.
 *
 *   time_ms > 0       Stop once this many wall-clock milliseconds have
 *                     elapsed, returning the best estimate accumulated so far.
 *
- * When multiple fields are non-zero the first criterion to fire wins,
- * giving callers fine control over the accuracy / speed trade-off.
+ * When multiple stopping fields are non-zero, the first criterion to fire
+ * wins, giving callers fine control over the accuracy / speed trade-off.
 *
+ * @c stability_metrics selects which estimates must meet @c stability_mm.
+ * The default derives them from the requested output arguments; visitor APIs
+ * return both statistics and therefore consider both unless told otherwise.
+ *
+ * The optional progress callback runs serially after completed ray batches.
+ * It is informational and must not call back into the active Crofton
+ * session.
 */
+typedef void (*rt_crofton_progress_fn)(size_t ray_count,
+                                      size_t crossing_count,
+                                      double surface_area,
+                                      double volume,
+                                      double stability_mm,
+                                      int stability_evaluated,
+                                      double elapsed_ms,
+                                      void *data);
+
+/** Estimates considered by stability-based stopping. */
+enum rt_crofton_stability_metric {
+    RT_CROFTON_STABILITY_DEFAULT = 0,
+    RT_CROFTON_STABILITY_SURFACE_AREA = 1,
+    RT_CROFTON_STABILITY_VOLUME = 2,
+    RT_CROFTON_STABILITY_ALL =
+        RT_CROFTON_STABILITY_SURFACE_AREA | RT_CROFTON_STABILITY_VOLUME
+};
+
 struct rt_crofton_params {
    size_t n_rays;       /**< max total rays; 0 = no limit via this criterion */
    double stability_mm; /**< equivalent-radius stability target (mm); 0 = disabled */
    double time_ms;      /**< wall-clock time budget (ms); 0 = disabled */
+    unsigned int stability_metrics; /**< rt_crofton_stability_metric bitmask */
+    rt_crofton_progress_fn progress; /**< optional serial progress callback */
+    void *progress_data; /**< caller data for progress */
 };

 /** One solid partition observed by a Crofton ray. */
 struct rt_crofton_segment {
+    point_t ray_origin;       /**< sampled line origin */
+    vect_t ray_direction;     /**< sampled unit line direction */
    point_t in_point;
    vect_t in_normal;
    point_t out_point;
    vect_t out_normal;
+    fastf_t in_distance;      /**< distance from ray_origin to in_point */
+    fastf_t out_distance;     /**< distance from ray_origin to out_point */
    fastf_t thickness;
    size_t ray_id;
+    const struct region *region; /**< owning region, valid while rtip lives */
+    const struct soltab *in_solid;  /**< solid producing the entry hit */
+    const struct soltab *out_solid; /**< solid producing the exit hit */
 };

 /**
@@ -250,9 +285,83 @@ struct rt_crofton_result {
    double volume;
 };

-/** Initialize an rt_crofton_result before first use. */
 #define RT_CROFTON_RESULT_INIT { NULL, 0, 0, 0, 0.0, 0.0 }

+/** Explicit sequence selection for experimental Crofton callers. */
+enum rt_crofton_sequence {
+    RT_CROFTON_SEQUENCE_DEFAULT = 0,
+    RT_CROFTON_SEQUENCE_RANDOM,
+    RT_CROFTON_SEQUENCE_QMC
+};
+
+/** Criterion that ended a Crofton sampling call. */
+enum rt_crofton_stop_reason {
+    RT_CROFTON_STOP_NONE = 0,
+    RT_CROFTON_STOP_RAYS,
+    RT_CROFTON_STOP_STABILITY,
+    RT_CROFTON_STOP_TIME
+};
+
+/** Estimates from a Crofton run that does not retain every segment. */
+struct rt_crofton_stats {
+    size_t ray_count;
+    size_t crossing_count;
+    double surface_area;
+    double volume;
+    double stability_mm; /**< last observed equivalent-radius change */
+    int stability_evaluated; /**< non-zero when stability_mm was observed */
+    size_t invalid_partition_count; /**< rejected non-finite partitions */
+    enum rt_crofton_stop_reason stop_reason; /**< criterion that ended sampling */
+};
+
+/** Opaque reusable worker session for serial Crofton sampling calls. */
+struct rt_crofton_session;
+
+/** Called serially after each ray batch; the segment is borrowed. */
+typedef void (*rt_crofton_segment_fn)(
+    const struct rt_crofton_segment *segment, void *data);
+
+/** A complete sampled ray and its valid positive-length partitions. */
+struct rt_crofton_ray {
+    size_t ray_id;
+    point_t origin;
+    vect_t direction;
+    const struct rt_crofton_segment *segments;
+    size_t segment_count;
+};
+
+/** Called serially after each ray batch; the ray and segments are borrowed. */
+typedef void (*rt_crofton_ray_fn)(const struct rt_crofton_ray *ray, void *data);
+
+/** A raw overlap event observed before boolean overlap resolution. */
+struct rt_crofton_overlap {
+    size_t ray_id;
+    point_t origin;
+    vect_t direction;
+    point_t in_point;
+    point_t out_point;
+    fastf_t depth;
+    const struct region *first_region;
+    const struct region *second_region;
+};
+
+/** Called serially after each ray batch for each raw overlap event. */
+typedef void (*rt_crofton_overlap_fn)(const struct rt_crofton_overlap *overlap,
+                                      void *data);
+
+struct rt_crofton_invalid_ray {
+    size_t ray_id;
+    point_t origin;
+    vect_t direction;
+    fastf_t in_distance;
+    fastf_t out_distance;
+    fastf_t thickness;
+    const struct region *region;
+    const char *reason;
+};
+
+typedef void (*rt_crofton_invalid_fn)(const struct rt_crofton_invalid_ray *ray,
+                                      void *data);

 /**
 * Run the Cauchy-Crofton ray-sampling estimator on an already-prepared
@@ -284,6 +393,26 @@ struct rt_crofton_result {
 *          bound fit failure.  A return
 *          value of 0 means no geometry was intersected by the sampler.
 *
+ * @section crofton_sampler Sampling method
+ *
+ * When BRL-CAD is built with BRLCAD_ENABLE_QMC, Crofton sampling defaults to
+ * OpenQMC's randomized rank-one lattice.  It maps four uniformly distributed
+ * coordinates directly to the two direction and two perpendicular-offset
+ * coordinates of an oriented line.  Directions are isotropic, and the offset
+ * radius uses R*sqrt(u) so lines are uniform by area across the bounding disk,
+ * not concentrated toward its center.  This is the same invariant line
+ * measure used by the independent-random-endpoint construction, so the
+ * Crofton area and volume formulas and stopping criteria do not change.
+ *
+ * The lattice improves finite-sample coverage but its samples are correlated,
+ * unlike ordinary Monte Carlo samples.  OpenQMC shuffles indices and applies
+ * toroidal shifts to mitigate alignment artifacts; a fixed pattern schedule
+ * keeps runs reproducible.  Set the LIBRT_CROFTON_USE_PRNG environment
+ * variable to a value accepted by bu_str_true() to force the original
+ * independent-PRNG generator.  That path is useful for compatibility and for
+ * checking whether an unexpected result is specific to lattice correlation.
+ *
+ *
 * @section crofton_near_tol Near-tolerance sliver geometry: CSG vs BoT divergence
 *
 * The Cauchy-Crofton formula is mathematically exact for any well-defined
@@ -379,6 +508,136 @@ RT_EXPORT extern int rt_crofton_collect(struct rt_crofton_result *result,
 					const fastf_t *bbox_min,
 					const fastf_t *bbox_max);

+/** Whether the librt build supports the QMC Crofton sequence. */
+RT_EXPORT extern int rt_crofton_qmc_available(void);
+
+/** Resolve DEFAULT to the sequence selected by this process. */
+RT_EXPORT extern enum rt_crofton_sequence rt_crofton_resolve_sequence(
+                                            enum rt_crofton_sequence sequence);
+
+/**
+ * Visit region-attributed segments without retaining the full sample set.
+ * The callback runs serially between bounded ray batches.  It receives every
+ * positive-length partition; a normal may be zero when the raytracer cannot
+ * calculate one.  Borrowed pointers are valid only during the callback,
+ * except for the region pointer, which remains valid while @p rtip is alive.
+ */
+RT_EXPORT extern int rt_crofton_visit(struct rt_crofton_stats *stats,
+                                     struct rt_i *rtip,
+                                     const struct rt_crofton_params *params,
+                                     size_t ray_offset,
+                                     const fastf_t *bbox_min,
+                                     const fastf_t *bbox_max,
+                                     enum rt_crofton_sequence sequence,
+                                     rt_crofton_segment_fn visitor,
+                                     rt_crofton_invalid_fn invalid_visitor,
+                                     void *data);
+
+/**
+ * Seeded form of rt_crofton_visit().  A zero seed selects the library default;
+ * distinct stream identifiers select reproducible randomized replications.
+ * Independent replications, rather than offsets within one point set, are
+ * required when estimating randomized-QMC sampling error.
+ */
+RT_EXPORT extern int rt_crofton_visit_seeded(
+                                     struct rt_crofton_stats *stats,
+                                     struct rt_i *rtip,
+                                     const struct rt_crofton_params *params,
+                                     size_t ray_offset,
+                                     const fastf_t *bbox_min,
+                                     const fastf_t *bbox_max,
+                                     enum rt_crofton_sequence sequence,
+                                     uint64_t seed,
+                                     uint64_t stream_id,
+                                     rt_crofton_segment_fn visitor,
+                                     rt_crofton_invalid_fn invalid_visitor,
+                                     void *data);
+
+/**
+ * Create a reusable Crofton worker session for one prepared raytrace instance.
+ *
+ * A session retains its parallel worker pool and per-CPU raytrace resources
+ * across calls.  Calls using a session must be serialized, and @p rtip must
+ * remain alive and unchanged until rt_crofton_session_destroy() returns.
+ * This facility avoids repeatedly starting an all-CPU worker pool when a
+ * caller evaluates several independent randomized replications.
+ */
+RT_EXPORT extern struct rt_crofton_session *rt_crofton_session_create(
+                                     struct rt_i *rtip);
+
+/** Release a reusable Crofton worker session and its raytrace resources. */
+RT_EXPORT extern void rt_crofton_session_destroy(
+                                     struct rt_crofton_session *session);
+
+/**
+ * Session-backed form of rt_crofton_visit_seeded().
+ *
+ * The seed and stream semantics are identical to
+ * rt_crofton_visit_seeded().  Keeping replications on distinct stream
+ * identifiers is required for randomized-QMC error estimation; reusing the
+ * worker session does not combine or otherwise correlate their samples.
+ */
+RT_EXPORT extern int rt_crofton_session_visit_seeded(
+                                     struct rt_crofton_session *session,
+                                     struct rt_crofton_stats *stats,
+                                     const struct rt_crofton_params *params,
+                                     size_t ray_offset,
+                                     const fastf_t *bbox_min,
+                                     const fastf_t *bbox_max,
+                                     enum rt_crofton_sequence sequence,
+                                     uint64_t seed,
+                                     uint64_t stream_id,
+                                     rt_crofton_segment_fn visitor,
+                                     rt_crofton_invalid_fn invalid_visitor,
+                                     void *data);
+
+/**
+ * Visit complete rays, including rays with no valid partitions.  Segment
+ * pointers are borrowed and valid only during the callback.  The ray record
+ * contains only positive-length partitions; rejected partitions are reported
+ * through @p invalid_visitor when supplied.
+ */
+RT_EXPORT extern int rt_crofton_visit_rays(struct rt_crofton_stats *stats,
+                                           struct rt_i *rtip,
+                                           const struct rt_crofton_params *params,
+                                           size_t ray_offset,
+                                           const fastf_t *bbox_min,
+                                           const fastf_t *bbox_max,
+                                           enum rt_crofton_sequence sequence,
+                                           rt_crofton_ray_fn ray_visitor,
+                                           rt_crofton_invalid_fn invalid_visitor,
+                                           void *data);
+
+/** Extended ray visitor with raw overlap events. */
+RT_EXPORT extern int rt_crofton_visit_rays_ex(struct rt_crofton_stats *stats,
+                                              struct rt_i *rtip,
+                                              const struct rt_crofton_params *params,
+                                              size_t ray_offset,
+                                              const fastf_t *bbox_min,
+                                              const fastf_t *bbox_max,
+                                              enum rt_crofton_sequence sequence,
+                                              rt_crofton_ray_fn ray_visitor,
+                                              rt_crofton_overlap_fn overlap_visitor,
+                                              rt_crofton_invalid_fn invalid_visitor,
+                                              void *data);
+
+/** Seeded form of rt_crofton_visit_rays_ex(). */
+RT_EXPORT extern int rt_crofton_visit_rays_seeded_ex(
+                                              struct rt_crofton_stats *stats,
+                                              struct rt_i *rtip,
+                                              const struct rt_crofton_params *params,
+                                              size_t ray_offset,
+                                              const fastf_t *bbox_min,
+                                              const fastf_t *bbox_max,
+                                              enum rt_crofton_sequence sequence,
+                                              uint64_t seed,
+                                              uint64_t stream_id,
+                                              rt_crofton_ray_fn ray_visitor,
+                                              rt_crofton_overlap_fn overlap_visitor,
+                                              rt_crofton_invalid_fn invalid_visitor,
+                                              void *data);
+
+
 /** Release storage owned by a structured Crofton result. */
 RT_EXPORT extern void rt_crofton_result_free(
 					struct rt_crofton_result *result);
@@ -387,8 +646,8 @@ RT_EXPORT extern void rt_crofton_result_free(
 /**
 * Cauchy-Crofton surface-area and/or volume estimator for a primitive.
 *
- * Creates a temporary in-memory raytrace of @p ip and fires random chord
- * rays until the stopping criteria in @p params are satisfied (or all
+ * Creates a temporary in-memory raytrace of @p ip and samples Crofton lines
+ * until the stopping criteria in @p params are satisfied (or all
 * criteria are zero / params is NULL, in which case the 2 000-ray default
 * is used).  Stores the estimated surface area in @p *area (if non-NULL)
 * and the estimated volume in @p *vol (if non-NULL).
@@ -258,6 +258,7 @@ function(BRLCAD_Summary)
  ext_report("Netpbm" NETPBM_STATUS)
  ext_report("OpenCV" OPENCV_STATUS)
  ext_report("OpenNURBS" OPENNURBS_STATUS)
+  ext_report("OpenQMC" OPENQMC_STATUS REQUIRED_VARS "BRLCAD_ENABLE_QMC")
  ext_report("OSMesa" OSMESA_STATUS)
  ext_report("Portable Network Graphics" PNG_STATUS)

@@ -308,6 +308,10 @@ if(BRLCAD_ENABLE_OPENCL)
  endif(NOT OPENCL_FOUND)
 endif(BRLCAD_ENABLE_OPENCL)

+# Enable quasi-Monte Carlo sampling where supported
+option(BRLCAD_ENABLE_QMC "Enable quasi-Monte Carlo sampling" OFF)
+mark_as_advanced(BRLCAD_ENABLE_QMC)
+
 # Enable features requiring OpenVDB
 option(BRLCAD_ENABLE_OPENVDB "Enable features requiring OpenVDB" OFF)
 mark_as_advanced(BRLCAD_ENABLE_OPENVDB)
@@ -300,7 +300,9 @@ crofton_estimate(const char *gfile, const char *obj_name,
    rt_prep_parallel(rtip, 1);

    /* REGRESS_CROFTON_SAMPLES rays, REGRESS_CROFTON_PCT% convergence threshold */
-    struct rt_crofton_params p = { REGRESS_CROFTON_SAMPLES, REGRESS_CROFTON_STABLE_TARGET, 0.0};
+    struct rt_crofton_params p = {
+        REGRESS_CROFTON_SAMPLES, REGRESS_CROFTON_STABLE_TARGET, 0.0,
+        RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
    int ret = rt_crofton_shoot(sa, vol, NULL, NULL, NULL, NULL, NULL,
 			       rtip, &p, NULL, NULL);
    rt_i_destroy(rtip);
@@ -4,10 +4,10 @@ if(TARGET gqa AND TARGET bview_plot3 AND TARGET plot3stat)
  brlcad_regression_test(regress-gqa_plot "bview_plot3;gqa;plot3stat" TEST_SCRIPT "${TSCRIPT}" EXEC gqa VEXEC bview_plot3)
 endif(TARGET gqa AND TARGET bview_plot3 AND TARGET plot3stat)

-if(SH_EXEC AND TARGET mged AND TARGET gqa)
+if(SH_EXEC AND TARGET mged AND TARGET gqa AND TARGET nirt)
  brlcad_add_test(NAME regress-gqa COMMAND ${SH_EXEC} "${CMAKE_CURRENT_SOURCE_DIR}/gqa.sh" ${CMAKE_SOURCE_DIR})
-  brlcad_regression_test(regress-gqa "gqa;mged" TEST_DEFINED)
-endif(SH_EXEC AND TARGET mged AND TARGET gqa)
+  brlcad_regression_test(regress-gqa "gqa;mged;nirt" TEST_DEFINED)
+endif(SH_EXEC AND TARGET mged AND TARGET gqa AND TARGET nirt)

 # list of temporary files
 set(
@@ -18,6 +18,53 @@ set(
  gqa.g
  gqa.gaps.plot3
  gqa.log
+  gqa.matrix.g
+  gqa.matrix.mged
+  gqa.modern.options.out
+  gqa.modern.invalid_value.out
+  gqa.modern.unknown_option.out
+  gqa.modern.grid.out
+  gqa.modern.grid.json
+  gqa.modern.collision.out
+  gqa.modern.alias_collision.out
+  gqa.alias.json
+  gqa.modern.grid_edge.out
+  gqa.modern.rotated.out
+  gqa.modern.rotated.json
+  gqa.modern.grid_defaults.out
+  gqa.modern.grid_defaults.json
+  gqa.modern.axis_angle.out
+  gqa.modern.refine.out
+  gqa.modern.refine.json
+  gqa.modern.air_excluded.out
+  gqa.modern.matrix.out
+  gqa.modern.random.out
+  gqa.modern.uncertainty.out
+  gqa.modern.uncertainty.json
+  gqa.modern.centroid.out
+  gqa.modern.centroid.json
+  gqa.modern.accuracy.out
+  gqa.modern.accuracy.json
+  gqa.modern.absolute.out
+  gqa.modern.absolute.json
+  gqa.modern.seed_a.out
+  gqa.modern.seed_b.out
+  gqa.modern.issues.out
+  gqa.modern.issues.json
+  gqa.modern.crofton.out
+  gqa.modern.crofton.json
+  gqa.modern.tolerance.out
+  gqa.modern.stability_cap.out
+  gqa.modern.defaults.out
+  gqa.modern.defaults.json
+  gqa.modern.deadline.out
+  gqa.modern.deadline.json
+  gqa.modern.time_limit.out
+  gqa.modern.time_limit.json
+  gqa.modern.gaps.out
+  gqa.modern.gaps.json
+  gqa.modern.air.out
+  gqa.modern.air.json
  gqa.mass_air.am.out
  gqa.mass_air.aw.out
  gqa.mass_mix.am.out
@@ -25,6 +72,8 @@ set(
  gqa.mass_mix.u0.out
  gqa.mged
  gqa.overlaps.plot3
+  gqa.replay.json
+  gqa.replay.out
  gqa.sphere_adaptive.out
  gqa.volume.plot3
  ovlp_overlaps.plot3
@@ -112,6 +112,11 @@ if test ! -f "$GQABIN" ; then
    log "Unable to find gqa, aborting"
    exit 1
 fi
+NIRT="`ensearch nirt`"
+if test ! -f "$NIRT" ; then
+    log "Unable to find nirt, aborting"
+    exit 1
+fi

 rm -f density_table.txt
 echo "5 1 stuff" > density_table.txt
@@ -169,6 +174,23 @@ in mass_solid.s rpp 0 1 0 1 0 1
 r mass_solid.r u mass_solid.s
 adjust mass_solid.r GIFTmater 5

+# Half-spaces may split finite geometry but do not independently define a
+# measurable volume.
+in split.s half 1 0 0 0.5
+r split_intersection.r u mass_solid.s + split.s
+r split_subtraction.r u mass_solid.s - split.s
+g split.g split_intersection.r split_subtraction.r
+r unbounded.r u split.s
+in bounded_submodel.s submodel split_intersection.r 0 ""
+r bounded_submodel.r u bounded_submodel.s
+in unbounded_submodel.s submodel unbounded.r 0 ""
+r unbounded_submodel.r u unbounded_submodel.s
+
+# A 0.01 mm solid verifies that automatic stability scales below
+# BN_TOL_DIST rather than treating sub-tolerance geometry as ordinary size.
+in tiny.s rpp 0 0.00001 0 0.00001 0 0.00001
+r tiny.r u tiny.s
+
 in mass_air.s rpp 1 2 0 1 0 1
 r mass_air.r u mass_air.s
 adjust mass_air.r air 2
@@ -191,6 +213,19 @@ $MGED -c gqa.g <<EOF >> $LOGFILE 2>&1
 `cat gqa.mged`
 EOF

+cat > gqa.matrix.mged <<EOF
+units m
+in box.s rpp 0 1 0 1 0 1
+r left.r u box.s
+r right.r u box.s
+comb moved.c u right.r
+arced moved.c/right.r matrix rarc xlate 0.5 0 0
+g model left.r moved.c
+q
+EOF
+rm -f gqa.matrix.g
+$MGED -c gqa.matrix.g < gqa.matrix.mged >> $LOGFILE 2>&1
+
 GQA="$GQABIN -u m,m^3,kg -g 250mm-50mm -p gqa."
 STATUS=0

@@ -259,6 +294,361 @@ assert_close 1000 "$MASS_MIX_U0" 0.001 "-U 0 excludes modeled air from weight"
 assert_close 4.18879 "$SPHERE_VOLUME" 0.01 "adaptive sphere volume refines beyond symmetric first pass"


+# The experimental interface must remain separate from the legacy options.
+run_capture gqa.modern.grid.out $GQABIN --analyze --measure volume,mass,area --spacing 50 --refine 0 --json gqa.modern.grid.json gqa.g mass_solid.r
+MODERN_GRID_VOLUME="`extract_last_number 'Volume:' gqa.modern.grid.out`"
+MODERN_GRID_MASS="`extract_last_number 'Mass:' gqa.modern.grid.out`"
+assert_close 1000000000 "$MODERN_GRID_VOLUME" 1000 "experimental grid cube volume"
+assert_close 1000000 "$MODERN_GRID_MASS" 1 "experimental grid cube mass"
+if test ! -s gqa.modern.grid.json || ! grep -q '"schema_version": 2' gqa.modern.grid.json || ! grep -q '"path": "/mass_solid.r"' gqa.modern.grid.json ; then
+    log "FAIL: experimental JSON region record"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.bounds.out $GQABIN --analyze --measure bounds --json gqa.modern.bounds.json gqa.g mass_solid.r
+if ! grep -q '^gqa analysis (grid), 0 rays$' gqa.modern.bounds.out ||
+   ! grep -q '^Bounding box: 0 0 0  1000 1000 1000 mm$' gqa.modern.bounds.out ||
+   ! grep -q '"face_areas_mm2":' gqa.modern.bounds.json ; then
+    log "FAIL: experimental Boolean bounds report"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.plot_all.out $GQABIN --analyze --measure volume --spacing 50 --refine 1 --plot-prefix gqa.modern.plot.all. --plot-mode all gqa.g mass_solid.r
+run_capture gqa.modern.plot_representative.out $GQABIN --analyze --measure volume --spacing 50 --refine 1 --plot-prefix gqa.modern.plot.representative. --plot-mode representative --plot-resolution 250 gqa.g mass_solid.r
+ALL_PLOT_SIZE="`wc -c < gqa.modern.plot.all.volume.plot3 | tr -d ' '`"
+REPRESENTATIVE_PLOT_SIZE="`wc -c < gqa.modern.plot.representative.volume.plot3 | tr -d ' '`"
+if ! cmp -s gqa.modern.plot_all.out gqa.modern.plot_representative.out ||
+   test ! -s gqa.modern.plot.all.volume.plot3 ||
+   test ! -s gqa.modern.plot.representative.volume.plot3 ||
+   test "$REPRESENTATIVE_PLOT_SIZE" -ge "$ALL_PLOT_SIZE" ; then
+    log "FAIL: experimental all and representative volume plots"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.plot_overlap.out $GQABIN --analyze --check overlaps --spacing 250 --refine 0 --plot-prefix gqa.modern.plot.overlap. gqa.g overlaps
+run_capture gqa.modern.plot_gap.out $GQABIN --analyze --check gaps --spacing 50 --refine 0 --plot-prefix gqa.modern.plot.gap. gqa.g gap.g
+run_capture gqa.modern.plot_adj_air.out $GQABIN --analyze --check adjacent-air --spacing 250 --refine 0 --plot-prefix gqa.modern.plot.air. gqa.g adj_air.g
+run_capture gqa.modern.plot_exp_air.out $GQABIN --analyze --check exposed-air --spacing 250 --refine 0 --plot-prefix gqa.modern.plot.air. gqa.g exposed_air.g
+if test ! -s gqa.modern.plot.overlap.overlaps.plot3 ||
+   test ! -s gqa.modern.plot.gap.gaps.plot3 ||
+   test ! -s gqa.modern.plot.air.adj_air.plot3 ||
+   test ! -s gqa.modern.plot.air.exp_air.plot3 ; then
+    log "FAIL: experimental diagnostic Plot3 output"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+if $GQABIN --analyze --measure volume --spacing 50 --refine 0 --plot-view gqa.g mass_solid.r > gqa.modern.plot_view.out 2>&1 ; then
+    log "FAIL: standalone gqa accepted --plot-view"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+$MGED -c gqa.g > gqa.modern.mged_plot_view.out 2>&1 <<EOF
+gqa --analyze --measure volume --spacing 100 --refine 0 --plot-view mass_solid.r
+view objs
+quit
+EOF
+cat gqa.modern.mged_plot_view.out >> "$LOGFILE"
+if ! grep -q '^Volume: 1e+09 mm\^3$' gqa.modern.mged_plot_view.out ||
+   ! grep -q '^gqa::analysis$' gqa.modern.mged_plot_view.out ; then
+    log "FAIL: MGED live analysis plot"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+GQA_DB_CHECKSUM_BEFORE="`cksum gqa.g`"
+if $GQABIN --analyze --measure volume --json gqa.g gqa.g mass_solid.r > gqa.modern.collision.out 2>&1 ; then
+    log "FAIL: analysis accepted the database as its JSON output"
+    STATUS="`expr $STATUS + 1`"
+fi
+GQA_DB_CHECKSUM_AFTER="`cksum gqa.g`"
+if test "x$GQA_DB_CHECKSUM_BEFORE" != "x$GQA_DB_CHECKSUM_AFTER" ; then
+    log "FAIL: rejected output collision modified the database"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+if $GQABIN --analyze --measure volume --json gqa.alias.json --invalid-rays ./gqa.alias.json gqa.g mass_solid.r > gqa.modern.alias_collision.out 2>&1 ; then
+    log "FAIL: analysis accepted aliased JSON and invalid-ray output paths"
+    STATUS="`expr $STATUS + 1`"
+fi
+if test -e gqa.alias.json ; then
+    log "FAIL: rejected aliased output paths created an output file"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+# The shared bu_opt parser must recognize assignment syntax and options on
+# either side of positional database and object arguments.
+run_capture gqa.modern.options.out $GQABIN --analyze --measure=volume gqa.g --sampler=crofton --sequence=random mass_solid.r --rays=1000
+if ! grep -q '^gqa analysis (crofton, random), 1000 rays$' gqa.modern.options.out ||
+   ! grep -q '^Stopping condition: ray_limit$' gqa.modern.options.out ; then
+    log "FAIL: experimental bu_opt syntax and positional argument handling"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+# Exercise QMC when this build provides it.  A failed probe is a valid skip
+# only when gqa explicitly reports that QMC support is unavailable.
+log "... checking for QMC support"
+QMC_STATUS=0
+$GQABIN --analyze --measure=volume --sampler=crofton --sequence=qmc --rays=1000 gqa.g mass_solid.r > gqa.modern.qmc.out 2>&1 || QMC_STATUS=$?
+cat gqa.modern.qmc.out >> "$LOGFILE"
+if test $QMC_STATUS -eq 0 ; then
+    if ! grep -q '^gqa analysis (crofton, qmc), 1000 rays$' gqa.modern.qmc.out ||
+       ! grep -q '^Stopping condition: ray_limit$' gqa.modern.qmc.out ; then
+	log "FAIL: experimental QMC analysis"
+	STATUS="`expr $STATUS + 1`"
+    fi
+else
+    if ! grep -q '^QMC is unavailable in this build\.$' gqa.modern.qmc.out ; then
+	log "FAIL: experimental QMC availability check"
+	STATUS="`expr $STATUS + 1`"
+    fi
+fi
+
+if $GQABIN --analyze --measure=invalid gqa.g mass_solid.r > gqa.modern.invalid_value.out 2>&1 ; then
+    log "FAIL: analysis accepted an invalid bu_opt value"
+    STATUS="`expr $STATUS + 1`"
+fi
+if $GQABIN --analyze --unknown-option gqa.g mass_solid.r > gqa.modern.unknown_option.out 2>&1 ; then
+    log "FAIL: analysis accepted an unknown bu_opt option"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.grid_edge.out $GQABIN --analyze --measure volume,area --spacing 600 --refine 0 gqa.g mass_solid.r
+MODERN_GRID_EDGE_VOLUME="`extract_last_number 'Volume:' gqa.modern.grid_edge.out`"
+MODERN_GRID_EDGE_AREA="`extract_last_number 'Area:' gqa.modern.grid_edge.out`"
+assert_close 1000000000 "$MODERN_GRID_EDGE_VOLUME" 1000 "experimental clipped grid-cell volume"
+assert_close 6000000 "$MODERN_GRID_EDGE_AREA" 1 "experimental clipped grid-cell area"
+
+run_capture gqa.modern.half_intersection.out $GQABIN --analyze --measure volume --spacing 50 --refine 0 gqa.g split_intersection.r
+run_capture gqa.modern.half_subtraction.out $GQABIN --analyze --measure volume --spacing 50 --refine 0 gqa.g split_subtraction.r
+run_capture gqa.modern.half_split.out $GQABIN --analyze --measure volume --spacing 50 --refine 0 gqa.g split.g
+HALF_INTERSECTION_VOLUME="`extract_last_number 'Volume:' gqa.modern.half_intersection.out`"
+HALF_SUBTRACTION_VOLUME="`extract_last_number 'Volume:' gqa.modern.half_subtraction.out`"
+HALF_SPLIT_VOLUME="`extract_last_number 'Volume:' gqa.modern.half_split.out`"
+assert_close 500000000 "$HALF_INTERSECTION_VOLUME" 1000 "experimental half-space intersection volume"
+assert_close 500000000 "$HALF_SUBTRACTION_VOLUME" 1000 "experimental half-space subtraction volume"
+assert_close 1000000000 "$HALF_SPLIT_VOLUME" 1000 "experimental half-space split volume"
+
+run_capture gqa.modern.half_crofton.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --seed 42 --rays 20000 gqa.g split.g
+HALF_CROFTON_VOLUME="`extract_last_number 'Volume:' gqa.modern.half_crofton.out`"
+assert_close 1000000000 "$HALF_CROFTON_VOLUME" 50000000 "Crofton half-space split volume"
+
+run_capture gqa.modern.half_submodel.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --seed 42 --rays 20000 gqa.g bounded_submodel.r
+HALF_SUBMODEL_VOLUME="`extract_last_number 'Volume:' gqa.modern.half_submodel.out`"
+assert_close 500000000 "$HALF_SUBMODEL_VOLUME" 25000000 "Crofton submodel half-space clipping volume"
+
+if $GQABIN --analyze --measure volume --spacing 50 --refine 0 gqa.g unbounded.r > gqa.modern.unbounded.out 2>&1 ; then
+    log "FAIL: experimental analysis accepted an unbounded region"
+    STATUS="`expr $STATUS + 1`"
+fi
+if ! grep -q '^Cannot analyze unbounded region /unbounded.r\.$' gqa.modern.unbounded.out ; then
+    log "FAIL: experimental unbounded-region diagnostic"
+    STATUS="`expr $STATUS + 1`"
+fi
+if $GQABIN --analyze --measure volume --sampler crofton --rays 100 gqa.g unbounded_submodel.r > gqa.modern.unbounded_submodel.out 2>&1 ; then
+    log "FAIL: experimental analysis accepted an unbounded submodel"
+    STATUS="`expr $STATUS + 1`"
+fi
+if ! grep -q '^Cannot analyze unbounded region /unbounded_submodel.r\.$' gqa.modern.unbounded_submodel.out ; then
+    log "FAIL: experimental unbounded-submodel diagnostic"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.refine.out $GQABIN --analyze --measure volume,area --spacing 600 --refine 1 --json gqa.modern.refine.json gqa.g mass_solid.r
+if ! grep -q '^Final refinement volume change:' gqa.modern.refine.out ||
+   ! grep -q '"level": 1' gqa.modern.refine.json ; then
+    log "FAIL: experimental grid refinement evidence"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.rotated.out $GQABIN --analyze --measure volume,mass,area,centroid,moments --density density_table.txt --sampler grid-rotated --spacing 50 --refine 2 --azimuth 35 --elevation 25 --json gqa.modern.rotated.json gqa.g mass_solid.r
+MODERN_ROTATED_VOLUME="`extract_last_number 'Volume:' gqa.modern.rotated.out`"
+MODERN_ROTATED_MASS="`extract_last_number 'Mass:' gqa.modern.rotated.out`"
+MODERN_ROTATED_AREA="`extract_last_number 'Area:' gqa.modern.rotated.out`"
+assert_close 1000000000 "$MODERN_ROTATED_VOLUME" 1000000 "experimental rotated-grid cube volume"
+assert_close 1000000 "$MODERN_ROTATED_MASS" 1000 "experimental rotated-grid cube mass"
+assert_close 6000000 "$MODERN_ROTATED_AREA" 10000 "experimental rotated-grid cube area"
+if ! grep -q '^gqa analysis (grid-rotated),' gqa.modern.rotated.out ||
+   ! grep -q '^Centroid:' gqa.modern.rotated.out ||
+   ! grep -q '^Inertia tensor' gqa.modern.rotated.out ||
+   ! grep -q '^Grid directional spread:' gqa.modern.rotated.out ||
+   ! grep -q '"sampler": "grid-rotated"' gqa.modern.rotated.json ||
+   ! grep -q '"centroid_mm":' gqa.modern.rotated.json ||
+   ! grep -q '"inertia_tensor_g_mm2":' gqa.modern.rotated.json ||
+   ! grep -q '"grid_view_directions":' gqa.modern.rotated.json ||
+   ! grep -q '"phase_u":' gqa.modern.rotated.json ; then
+    log "FAIL: experimental rotated-grid evidence"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.grid_defaults.out $GQABIN --analyze --measure volume --json gqa.modern.grid_defaults.json gqa.g mass_solid.r
+if ! grep -q '^Stopping condition: stable$' gqa.modern.grid_defaults.out ||
+   ! grep -q '"sampler": "grid"' gqa.modern.grid_defaults.json ||
+   ! grep -q '"spacing_defaulted": true' gqa.modern.grid_defaults.json ||
+   ! grep -q '"stability_defaulted": true' gqa.modern.grid_defaults.json ||
+   ! grep -q '"time_defaulted": true' gqa.modern.grid_defaults.json ||
+   ! grep -q '"time_ms": 60000.0' gqa.modern.grid_defaults.json ; then
+    log "FAIL: grid time-first automatic stability defaults"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+if $GQABIN --analyze --measure volume --sampler grid --azimuth 35 --refine 0 gqa.g mass_solid.r > gqa.modern.axis_angle.out 2>&1 ; then
+    log "FAIL: axis-aligned grid accepted a rotated-grid orientation"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.air_excluded.out $GQABIN --analyze --measure volume,mass --air exclude --spacing 50 --refine 0 gqa.g mass_air_mix.g
+MODERN_SOLID_VOLUME="`extract_last_number 'Volume:' gqa.modern.air_excluded.out`"
+MODERN_SOLID_MASS="`extract_last_number 'Mass:' gqa.modern.air_excluded.out`"
+assert_close 1000000000 "$MODERN_SOLID_VOLUME" 1000 "experimental air-excluded volume"
+assert_close 1000000 "$MODERN_SOLID_MASS" 1 "experimental air-excluded mass"
+
+run_capture gqa.modern.random.out $GQABIN --analyze --measure volume,area --sampler crofton --sequence random --rays 20000 gqa.g mass_solid.r
+MODERN_RANDOM_VOLUME="`extract_last_number 'Volume:' gqa.modern.random.out`"
+assert_close 1000000000 "$MODERN_RANDOM_VOLUME" 50000000 "experimental random Crofton cube volume"
+
+run_capture gqa.modern.uncertainty.out $GQABIN --analyze --measure volume,area --sampler crofton --sequence random --rays 3200 --uncertainty --accuracy-scope all --json gqa.modern.uncertainty.json gqa.g mass_solid.r
+if ! grep -q '^Approximate 95% sampling intervals (16 independent replicates):' gqa.modern.uncertainty.out ||
+   ! grep -q '"replicates": 16' gqa.modern.uncertainty.json ||
+   ! grep -q '"metric": "object:mass_solid.r.volume_mm3"' gqa.modern.uncertainty.json ||
+   ! grep -q '"metric": "region:/mass_solid.r#' gqa.modern.uncertainty.json ||
+   ! grep -q '"half_width":' gqa.modern.uncertainty.json ; then
+    log "FAIL: Crofton independent-replicate uncertainty"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.centroid.out $GQABIN --analyze --measure mass,centroid --density density_table.txt --sampler crofton --sequence random --rays 3200 --uncertainty --json gqa.modern.centroid.json gqa.g mass_solid.r
+if ! grep -q '^Centroid:' gqa.modern.centroid.out ||
+   ! grep -q '"metric": "model.centroid_x_mm"' gqa.modern.centroid.json ||
+   ! grep -q '"centroid_mm":' gqa.modern.centroid.json ; then
+    log "FAIL: Crofton paired-ratio centroid uncertainty"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.absolute.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --rays 6400 --absolute-error volume=100000000 --json gqa.modern.absolute.json gqa.g mass_solid.r
+if ! grep -q '^Accuracy target:' gqa.modern.absolute.out ||
+   ! grep -q '"volume": 100000000' gqa.modern.absolute.json ||
+   ! grep -q '"target_tolerance": 100000000' gqa.modern.absolute.json ; then
+    log "FAIL: Crofton absolute accuracy target"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.seed_a.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --seed 42 --rays 1000 gqa.g mass_solid.r
+run_capture gqa.modern.seed_b.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --seed 42 --rays 1000 gqa.g mass_solid.r
+if ! cmp -s gqa.modern.seed_a.out gqa.modern.seed_b.out ; then
+    log "FAIL: seeded Crofton sampling is not reproducible"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.defaults.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --seed 42 --json gqa.modern.defaults.json gqa.g tiny.r
+if ! grep -q '^Stopping condition: stable' gqa.modern.defaults.out ||
+   ! grep -q '"rays": null' gqa.modern.defaults.json ||
+   ! grep -q '"stability_defaulted": true' gqa.modern.defaults.json ||
+   ! grep -q '"stability_mm": 1e-05' gqa.modern.defaults.json ||
+   ! grep -q '"time_defaulted": true' gqa.modern.defaults.json ||
+   ! grep -q '"time_ms": 60000.0' gqa.modern.defaults.json ; then
+    log "FAIL: Crofton time-first automatic stability defaults"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.accuracy.out $GQABIN --analyze --measure volume,area --sampler crofton --sequence random --rays 6400 --uncertainty --relative-error 0.1 --json gqa.modern.accuracy.json gqa.g mass_solid.r
+if ! grep -q '^Accuracy target:' gqa.modern.accuracy.out ||
+   ! grep -q '"pilot_rays":' gqa.modern.accuracy.json ||
+   ! grep -q '"production_rays":' gqa.modern.accuracy.json ; then
+    log "FAIL: Crofton pilot and production accuracy evaluation"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.deadline.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --uncertainty --relative-error 0.1 --time 1000 --json gqa.modern.deadline.json gqa.g mass_solid.r
+if ! grep -q '^Accuracy target: target_met' gqa.modern.deadline.out ||
+   ! grep -q '^Adaptive plan:' gqa.modern.deadline.out ||
+   ! grep -q '^Deadline:' gqa.modern.deadline.out ||
+   ! grep -q '"rays": null' gqa.modern.deadline.json ||
+   ! grep -q '"forecast_production_rays":' gqa.modern.deadline.json ||
+   ! grep -q '"planning_safety_factor": 2' gqa.modern.deadline.json ||
+   ! grep -q '"deadline_hit": false' gqa.modern.deadline.json ; then
+    log "FAIL: Crofton deadline-driven adaptive accuracy"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.time_limit.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --uncertainty --relative-error 0.000000001 --time 500 --json gqa.modern.time_limit.json gqa.g mass_solid.r
+if ! grep -q '^Accuracy target: time_limit' gqa.modern.time_limit.out ||
+   ! grep -q '"accuracy_status": "time_limit"' gqa.modern.time_limit.json ||
+   ! grep -q '"deadline_limited": true' gqa.modern.time_limit.json ; then
+    log "FAIL: Crofton adaptive deadline status"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.issues.out $GQABIN --analyze --measure volume --check overlaps --spacing 250 --json gqa.modern.issues.json gqa.g overlaps
+if ! grep -q '^overlap: ' gqa.modern.issues.out ||
+   ! grep -q '"type": "overlap"' gqa.modern.issues.json ||
+   ! grep -q '"start_point_mm"' gqa.modern.issues.json ||
+   ! grep -q '"end_point_mm"' gqa.modern.issues.json ; then
+    log "FAIL: experimental overlap report"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.matrix.out $GQABIN --analyze --measure volume --check overlaps --spacing 100 gqa.matrix.g model
+if ! grep -q '/model/moved.c/right.r' gqa.modern.matrix.out ||
+   grep -q '^Area:' gqa.modern.matrix.out ; then
+    log "FAIL: transformed overlap path must remain complete"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+# Crofton diagnostics use complete ray partitions and must report the same
+# intentionally constructed overlap, gap, and air cases as the grid path.
+run_capture gqa.modern.crofton.out $GQABIN --analyze --sampler crofton --sequence random --rays 12000 --check overlaps --json gqa.modern.crofton.json gqa.g overlaps
+if ! grep -q 'Overlap candidate boxes: ' gqa.modern.crofton.out ||
+   ! grep -q '^overlap: ' gqa.modern.crofton.out ||
+   ! grep -q '"targeted_overlap": true' gqa.modern.crofton.json ||
+   ! grep -q '"event_probability_upper_bound_95"' gqa.modern.crofton.json ||
+   ! grep -q '"type": "overlap"' gqa.modern.crofton.json ; then
+    log "FAIL: Crofton targeted overlap report"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+CROFTON_EVENT_COUNT="`grep -c '"ray_id":' gqa.modern.crofton.json`"
+CROFTON_RAY_COUNT="`grep '"ray_id":' gqa.modern.crofton.json | sort -u | wc -l | tr -d ' '`"
+if test "$CROFTON_EVENT_COUNT" -ne "$CROFTON_RAY_COUNT" ||
+   test "$CROFTON_RAY_COUNT" -le 1 ; then
+    log "FAIL: Crofton overlap events must be delivered once with their own ray IDs"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.tolerance.out $GQABIN --analyze --sampler crofton --sequence random --rays 1000 --check overlaps --tolerance 1000000000 gqa.g overlaps
+if grep -q '^overlap: ' gqa.modern.tolerance.out ; then
+    log "FAIL: Crofton overlap tolerance"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.stability_cap.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --rays 100 --stability 1 gqa.g mass_solid.r
+if ! grep -q '), 100 rays$' gqa.modern.stability_cap.out ; then
+    log "FAIL: Crofton stability must honor the ray limit"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.gaps.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --rays 5000 --check gaps --json gqa.modern.gaps.json gqa.g gap.g
+if ! grep -q '^gap: ' gqa.modern.gaps.out ||
+   ! grep -q '"type": "gap"' gqa.modern.gaps.json ; then
+    log "FAIL: Crofton gap report"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+run_capture gqa.modern.air.out $GQABIN --analyze --measure volume --sampler crofton --sequence random --rays 5000 --check adjacent-air,exposed-air --json gqa.modern.air.json gqa.g adj_air.g exposed_air.g
+if ! grep -q '"sampler": "crofton"' gqa.modern.air.json ||
+   ! grep -q '"ray_count": 5000' gqa.modern.air.json ; then
+    log "FAIL: Crofton air diagnostic run"
+    STATUS="`expr $STATUS + 1`"
+fi
+
+cat > gqa.replay.json <<EOF
+{"schema_version":2,"model":["mass_solid.r"],"rays":[{"ray_id":7,"origin_mm":[1500,500,500],"direction":[-1,0,0]}]}
+EOF
+run_capture gqa.replay.out $NIRT --replay gqa.replay.json --ray-id 7 gqa.g mass_solid.r
+if ! grep -q 'Replaying ray 7' gqa.replay.out ||
+   ! grep -q 'mass_solid.r' gqa.replay.out ; then
+    log "FAIL: nirt replay of analysis ray record"
+    STATUS="`expr $STATUS + 1`"
+fi
+
 if [ $STATUS = 0 ] ; then
    log "-> gqa.sh succeeded"
 else
@@ -9,7 +9,7 @@ if (TARGET ged-check)
 endif()
 brlcad_addexec(gex "gex.cpp" "libbu;librt")
 brlcad_addexec(glint glint.cpp "librt;${M_LIBRARY}")
-brlcad_addexec(gqa gqa.c "libged;librt")
+brlcad_addexec(gqa gqa.c "libged;libanalyze;librt")

 brlcad_addexec(gtransfer gtransfer.c "librt;libpkg" NO_INSTALL)

@@ -21,8 +21,6 @@
 *
 * Perform quantitative analysis checks on geometry.
 *
- * XXX need to look at gap computation
- *
 * plot the points where overlaps start/stop
 *
 * Designed to be a framework for 3d sampling of the geometry volume.
@@ -30,6 +28,7 @@

 #include "common.h"

+#include "analyze/gqa.h"
 #include "bu/app.h"
 #include "bu/cmd.h"
 #include "bu/getopt.h"
@@ -51,6 +50,30 @@ main(int argc, char *argv[])
    bu_opterr = 0;
    bu_optind = 1;

+    if (argc > 1 && bu_strcmp(argv[1], "--analyze") == 0) {
+        int model_index = analyze_gqa_model_argument(argc,
+            (const char **)argv);
+        if (model_index < 0 || model_index + 1 >= argc)
+            bu_exit(1, "Usage: %s --analyze [options] model object [objects...]\n", argv[0]);
+        av = (const char **)bu_calloc(argc, sizeof(char *), "analysis argv");
+        for (i = j = 0; i < argc; i++) {
+            if (i != model_index)
+                av[j++] = argv[i];
+        }
+        av[0] = "gqa";
+        gedp = ged_open("db", argv[model_index], 1);
+        if (gedp == GED_NULL) {
+            bu_free((void *)av, "analysis argv");
+            bu_exit(1, "Cannot open %s\n", argv[model_index]);
+        }
+        c = ged_exec_gqa(gedp, j, av);
+        if (bu_vls_strlen(gedp->ged_result_str) > 0)
+            bu_log("%s", bu_vls_addr(gedp->ged_result_str));
+        ged_close(gedp);
+        bu_free((void *)av, "analysis argv");
+        return c == BRLCAD_OK ? 0 : 1;
+    }
+
    /* Get past command line options. */
    while ((c = bu_getopt(argc, argv, "A:a:de:f:g:Gn:N:p:P:qrS:t:U:u:vV:W:h?")) != -1) {
 	switch (c) {
@@ -32,6 +32,7 @@ set(
  contour.cpp
  density.cpp
  globals.c
+  gqa.cpp
  inside.c
  mass.c
  moments.c
@@ -42,6 +43,7 @@ set(
  overlaps.c
  polygonizer.c
  raydiff.c
+  sampling.c
  #find_subtracted_shapes.cpp
  surf_area.c
  util.cpp
@@ -0,0 +1,108 @@
+/*                     S A M P L I N G . C
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+
+#include "common.h"
+
+#include <float.h>
+#include <math.h>
+
+#include "analyze/sampling.h"
+
+/*
+ * Owen recommends at least ten independent randomized-QMC replications and
+ * an ordinary Student-t interval across their estimates.  See A. B. Owen,
+ * "Error estimation for quasi-Monte Carlo" (2025), section 5,
+ * https://arxiv.org/abs/2501.00150v3.
+ *
+ * Using t(0.975, 9) for all R >= 10 is conservative relative to the smaller
+ * critical value for R - 1 degrees of freedom and avoids making the public
+ * API depend on a special-functions package.
+ */
+#define ANALYZE_STUDENT_T_95_DF9 2.2621571627409915
+
+int
+analyze_sampling_interval95(struct analyze_sampling_interval *interval,
+    const double *replicate_estimates, size_t replicate_count)
+{
+    if (!interval || !replicate_estimates ||
+        replicate_count < ANALYZE_SAMPLING_MIN_REPLICATES)
+	return ANALYZE_ERROR;
+
+    double sum = 0.0;
+    for (size_t i = 0; i < replicate_count; i++) {
+	if (!isfinite(replicate_estimates[i])) return ANALYZE_ERROR;
+	sum += replicate_estimates[i];
+    }
+
+    interval->estimate = sum / (double)replicate_count;
+    double squared_deviation = 0.0;
+    for (size_t i = 0; i < replicate_count; i++) {
+	const double deviation = replicate_estimates[i] - interval->estimate;
+	squared_deviation += deviation * deviation;
+    }
+    const double variance = squared_deviation / (double)(replicate_count - 1);
+    interval->half_width = ANALYZE_STUDENT_T_95_DF9 *
+	sqrt(variance / (double)replicate_count);
+
+    return isfinite(interval->estimate) && isfinite(interval->half_width) ?
+	ANALYZE_OK : ANALYZE_ERROR;
+}
+
+int
+analyze_sampling_ratio_interval95(struct analyze_sampling_interval *interval,
+    const double *replicate_numerators,
+    const double *replicate_denominators, size_t replicate_count)
+{
+    if (!interval || !replicate_numerators || !replicate_denominators ||
+        replicate_count < ANALYZE_SAMPLING_MIN_REPLICATES)
+	return ANALYZE_ERROR;
+
+    double numerator_sum = 0.0;
+    double denominator_sum = 0.0;
+    for (size_t i = 0; i < replicate_count; i++) {
+	if (!isfinite(replicate_numerators[i]) ||
+	    !isfinite(replicate_denominators[i]))
+	    return ANALYZE_ERROR;
+	numerator_sum += replicate_numerators[i];
+	denominator_sum += replicate_denominators[i];
+    }
+    const double mean_numerator = numerator_sum / (double)replicate_count;
+    const double mean_denominator = denominator_sum / (double)replicate_count;
+    if (!isfinite(mean_numerator) || !isfinite(mean_denominator) ||
+	fabs(mean_denominator) <= DBL_MIN)
+	return ANALYZE_ERROR;
+
+    interval->estimate = mean_numerator / mean_denominator;
+    double squared_residual = 0.0;
+    for (size_t i = 0; i < replicate_count; i++) {
+	/* Linearization of Xbar/Ybar retains the paired X,Y covariance. */
+	const double residual = (replicate_numerators[i] -
+	    interval->estimate * replicate_denominators[i]) / mean_denominator;
+	squared_residual += residual * residual;
+    }
+    const double variance = squared_residual /
+	(double)(replicate_count - 1);
+    interval->half_width = ANALYZE_STUDENT_T_95_DF9 *
+	sqrt(variance / (double)replicate_count);
+
+    return isfinite(interval->estimate) && isfinite(interval->half_width) ?
+	ANALYZE_OK : ANALYZE_ERROR;
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -4,6 +4,7 @@ brlcad_add_test(NAME analyze_contour COMMAND analyze_contour)
 brlcad_addexec(analyze_raydiff raydiff.c "libanalyze;librt;libbu" TEST)
 brlcad_addexec(analyze_sp solid_partitions.c "libanalyze;librt;libbu" TEST)
 brlcad_addexec(analyze_nhit nhit.cpp "libanalyze;librt;libbu" TEST_USESDATA)
+brlcad_addexec(analyze_sampling sampling_test.c "libanalyze;libbn" TEST)

 #####################################
 #      analyze_densities testing    #
@@ -12,12 +13,14 @@ brlcad_addexec(analyze_densities density.cpp "libanalyze;libbu" TEST)

 brlcad_add_test(NAME analyze_densities_null        COMMAND analyze_densities)
 brlcad_add_test(NAME analyze_densities_std        COMMAND analyze_densities std)
+brlcad_add_test(NAME analyze_sampling_interval    COMMAND analyze_sampling)

 cmakefiles(
  CMakeLists.txt
  arbs.g
  contour.cpp
  raydiff.g
+  sampling_test.c
 )

 # Local Variables:
@@ -0,0 +1,98 @@
+/*                S A M P L I N G _ T E S T . C
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+
+#include "common.h"
+
+#include <float.h>
+#include <stdio.h>
+
+#include "analyze/gqa.h"
+#include "analyze/sampling.h"
+#include "bu/app.h"
+#include "vmath.h"
+
+
+int
+main(int UNUSED(argc), char *argv[])
+{
+    static const double estimates[ANALYZE_SAMPLING_MIN_REPLICATES] = {
+	1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
+    struct analyze_sampling_interval interval;
+    static const double ratio_numerators[ANALYZE_SAMPLING_MIN_REPLICATES] = {
+	1.0, 4.0, 9.0, 16.0, 25.0, 36.0, 49.0, 64.0, 81.0, 100.0};
+    static const double zero_denominators[ANALYZE_SAMPLING_MIN_REPLICATES] = {0.0};
+
+    bu_setprogname(argv[0]);
+
+    if (analyze_sampling_interval95(&interval, estimates,
+	    ANALYZE_SAMPLING_MIN_REPLICATES) != ANALYZE_OK ||
+	!NEAR_EQUAL(interval.estimate, 5.5, SMALL_FASTF) ||
+	!(interval.half_width > 0.0) || !(interval.half_width < DBL_MAX)) {
+	fprintf(stderr, "valid sampling interval was rejected\n");
+	return 1;
+    }
+    if (analyze_sampling_ratio_interval95(&interval, ratio_numerators,
+	    estimates, ANALYZE_SAMPLING_MIN_REPLICATES) != ANALYZE_OK ||
+	!NEAR_EQUAL(interval.estimate, 7.0, SMALL_FASTF) ||
+	!(interval.half_width > 0.0) || !(interval.half_width < DBL_MAX)) {
+	fprintf(stderr, "valid ratio sampling interval was rejected\n");
+	return 1;
+    }
+    if (analyze_sampling_ratio_interval95(&interval, ratio_numerators,
+	    zero_denominators, ANALYZE_SAMPLING_MIN_REPLICATES) != ANALYZE_ERROR) {
+	fprintf(stderr, "zero-denominator ratio interval was accepted\n");
+	return 1;
+    }
+    if (analyze_sampling_interval95(&interval, estimates,
+	    ANALYZE_SAMPLING_MIN_REPLICATES - 1) != ANALYZE_ERROR) {
+	fprintf(stderr, "undersized replicate set was accepted\n");
+	return 1;
+    }
+
+    {
+	const char *arguments[] = {
+	    "gqa", "--analyze", "--uncertainty", "--rays", "1600",
+	    "model.g", "object"};
+	if (analyze_gqa_model_argument(7, arguments) != 5) {
+	    fprintf(stderr, "analysis model argument was not located\n");
+	    return 1;
+	}
+    }
+    {
+	const char *arguments[] = {
+	    "gqa", "--analyze", "--rays", "model.g"};
+	if (analyze_gqa_model_argument(4, arguments) >= 0) {
+	    fprintf(stderr, "missing analysis option value was accepted\n");
+	    return 1;
+	}
+    }
+
+    return 0;
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -595,6 +595,11 @@ bv_vlblock_obj(struct bv_vlblock *vbp, struct bview *v, const char *name)
 	s = bv_obj_get(v, BV_VIEW_OBJS);
    }

+    /* Reset clears both fields, but replacement objects must retain their
+     * identity and view ownership before children are allocated. */
+    s->s_v = v;
+    bu_vls_sprintf(&s->s_name, "%s", name);
+
    for (size_t i = 0; i < vbp->nused; i++) {
 	if (!BU_LIST_IS_EMPTY(&(vbp->head[i]))) {
 	    struct bv_scene_obj *sc = bv_obj_get_child(s);
@@ -554,6 +554,8 @@ ged_init(struct ged *gedp)
    BU_GET(gedp->ged_cbs, struct ged_callback_state);
    gedp->ged_refresh_handler = NULL;
    gedp->ged_refresh_clientdata = NULL;
+    gedp->ged_run_on_main_thread = NULL;
+    gedp->ged_main_thread_clientdata = NULL;
    gedp->ged_output_handler = NULL;
    gedp->ged_create_vlist_scene_obj_callback = NULL;
    gedp->ged_create_vlist_display_list_callback = NULL;
@@ -1,7 +1,8 @@
-set(PLUGIN_LIBS libged libbu)
+set(PLUGIN_LIBS libged libanalyze libbu)

 ged_plugin_library(ged-gqa
  gqa.cpp
+  gqa_analyze.cpp
  LIBS ${PLUGIN_LIBS}
 )

@@ -19,14 +19,12 @@
 */
 /** @file libged/gqa.c
 *
- * performs a set of quantitative analyses on geometry.
+ * Performs quantitative analyses on geometry and optionally plots detected
+ * conditions.
 *
- * XXX need to look at gap computation
- *
- * plot the points where overlaps start/stop
- *
- * Designed to be a framework for 3d sampling of the geometry volume.
- * TODO: Need to move the sample pattern logic into LIBRT.
+ * This file retains the historical grid implementation for compatibility.
+ * The experimental samplers are shared by the standalone and GED commands
+ * through libanalyze, with Crofton ray generation supplied by librt.
 *
 */

@@ -294,7 +292,8 @@ static struct region_pair overlapList = {
 *
 * The standard units are millimeters, cubic millimeters, and grams.
 *
- * XXX this section should be extracted to libbu/units.c
+ * Keep this legacy command-output table local: its volume and mass aliases
+ * extend beyond libbu's general length-unit conversions.
 */
 struct cvt_tab {
    double val;
@@ -840,8 +839,6 @@ _gqa_overlap(struct application *ap,
 	bu_semaphore_release(state->sem_worker);
    }

-    /* XXX We should somehow flag the volume/weight calculations as invalid */
-
    /* since we have no basis to pick one over the other, just pick */
    return 1;	/* No further consideration to this partition */
 }
@@ -951,16 +948,12 @@ _gqa_hit(struct application *ap, struct partition *PartHeadp, struct seg *segs)
 	 */
 	if (analysis_flags & ANALYSIS_EXP_AIR) {

-	    /* FIXME: verify that the next partition is never
-	     * overlapping numerically with the current partition.
-	     * otherwise, we'll need to account for it. This debug
-	     * statement should be removed after confirming.
-	     * CSM@20220516
-	     */
+	    /* Boolean cleanup may leave the next partition beginning before
+	     * the current one exits.  Clip the exposed-air segment to that
+	     * entry rather than reporting a path through the overlap. */
 	    if (pp->pt_forw != PartHeadp) {
 		double next_dist = pp->pt_forw->pt_inhit->hit_dist - pp->pt_inhit->hit_dist;
 		if (next_dist < dist) {
-		    bu_log("DEBUG: next partition's entry is prior to current partition's exit\n");
 		    VJOIN1(opt, ap->a_ray.r_pt, pp->pt_forw->pt_inhit->hit_dist, ap->a_ray.r_dir);
 		}
 	    }
@@ -2551,7 +2544,14 @@ summary_reports(struct ged *gedp, struct cstate *state)
 	avg_mass /= num_views;
 	bu_vls_printf(gedp->ged_result_str, "  Average total volume: %g %s\n", avg_mass / units[VOL]->val, units[VOL]->name);
    }
-    if (analysis_flags & ANALYSIS_OVERLAPS) list_report(gedp, &overlapList);
+    if (analysis_flags & ANALYSIS_OVERLAPS) {
+	list_report(gedp, &overlapList);
+	if (BU_LIST_NON_EMPTY(&overlapList.l) &&
+	    (analysis_flags & (ANALYSIS_VOLUMES | ANALYSIS_WEIGHTS |
+		ANALYSIS_CENTROIDS | ANALYSIS_MOMENTS)))
+	    bu_vls_printf(gedp->ged_result_str,
+		"WARNING: overlaps invalidate reported volume and mass properties.\n");
+    }
    if (analysis_flags & ANALYSIS_ADJ_AIR) list_report(gedp, &adjAirList);
    if (analysis_flags & ANALYSIS_GAPS) list_report(gedp, &gapList);
    if (analysis_flags & ANALYSIS_EXP_AIR) list_report(gedp, &exposedAirList);
@@ -2596,9 +2596,13 @@ summary_reports(struct ged *gedp, struct cstate *state)
 }


+extern "C" int ged_gqa_analyze(struct ged *, int, const char **);
+
 extern "C" int
 ged_gqa_core(struct ged *gedp, int argc, const char *argv[])
 {
+    if (argc > 1 && !bu_strcmp(argv[1], "--analyze"))
+	return ged_gqa_analyze(gedp, argc, argv);
    int arg_count;
    struct rt_i *rtip;
    int i;
@@ -2632,16 +2636,12 @@ ged_gqa_core(struct ged *gedp, int argc, const char *argv[])
    azel_requested = 0;
    densityFileName = (char *)0;

-    /* FIXME: this is completely arbitrary, should probably be based
-     * on the model size.
-     */
+    /* Preserve the historical default for existing command lines.  The
+     * experimental interface derives its initial spacing from model scale. */
    gridSpacing = 50.0;

-    /* default grid spacing limit is based on the current distance
-     * tolerance, one order of magnitude greater.
-     *
-     * FIXME: should probably be based on the model size.
-     */
+    /* Preserve the historical tolerance-derived refinement limit.  The
+     * experimental interface instead uses dimensional stability. */
    gridSpacingLimit = 10.0 * wdbp->wdb_tol.dist;

    makeOverlapAssemblies = 0;
@@ -0,0 +1,130 @@
+/*                   G Q A _ A N A L Y Z E . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+
+#include "common.h"
+
+#include <utility>
+#include <vector>
+
+#include "analyze/gqa.h"
+#include "bv/vlist.h"
+#include "bu/log.h"
+#include "ged.h"
+#include "rt/vlist.h"
+#include "../ged_private.h"
+
+
+static int
+load_densities(struct analyze_densities **densities, char **source,
+    const char *filename, void *data)
+{
+    struct ged *gedp = static_cast<struct ged *>(data);
+    return _ged_read_densities(densities, source, gedp, filename, 0) ==
+        BRLCAD_OK ? ANALYZE_OK : ANALYZE_ERROR;
+}
+
+
+static void
+report_progress(const char *message, void *UNUSED(data))
+{
+    bu_log("%s\n", message);
+}
+
+
+struct PlotViewState {
+    struct ged *gedp = NULL;
+    std::vector<analyze_gqa_plot_line> lines;
+};
+
+static const char GQA_PLOT_VIEW_NAME[] = "gqa::analysis";
+
+
+struct PlotViewUpdate {
+    PlotViewState *state;
+    int action;
+    std::vector<analyze_gqa_plot_line> lines;
+};
+
+
+static void
+update_plot_view(void *data)
+{
+    PlotViewUpdate *update = static_cast<PlotViewUpdate *>(data);
+    PlotViewState &state = *update->state;
+    if (update->action == ANALYZE_GQA_PLOT_REPLACE)
+	state.lines = std::move(update->lines);
+    else
+	state.lines.insert(state.lines.end(), update->lines.begin(),
+	    update->lines.end());
+
+    struct bv_vlblock *vbp = rt_vlblock_init();
+    for (const auto &line : state.lines) {
+	struct bu_list *vhead = bv_vlblock_find(vbp,
+	    line.color[0], line.color[1], line.color[2]);
+	BV_ADD_VLIST(vbp->free_vlist_hd, vhead, line.start,
+	    BV_VLIST_LINE_MOVE);
+	BV_ADD_VLIST(vbp->free_vlist_hd, vhead, line.end,
+	    BV_VLIST_LINE_DRAW);
+    }
+    bv_vlblock_obj(vbp, state.gedp->ged_gvp, GQA_PLOT_VIEW_NAME);
+    bv_vlblock_free(vbp);
+    ged_refresh_cb(state.gedp);
+}
+
+
+static void
+report_plot(const struct analyze_gqa_plot_batch *batch, void *data)
+{
+    PlotViewState *state = static_cast<PlotViewState *>(data);
+    PlotViewUpdate update = {state, batch->action, {}};
+    if (batch->lines && batch->line_count)
+	update.lines.assign(batch->lines, batch->lines + batch->line_count);
+
+    struct ged *gedp = state->gedp;
+    gedp->ged_run_on_main_thread(gedp->ged_main_thread_clientdata,
+	update_plot_view, &update);
+}
+
+
+extern "C" int
+ged_gqa_analyze(struct ged *gedp, int argc, const char *argv[])
+{
+    GED_CHECK_DATABASE_OPEN(gedp, BRLCAD_ERROR);
+    PlotViewState plot_state;
+    plot_state.gedp = gedp;
+    struct analyze_gqa_context context = {
+        gedp->dbip, gedp->ged_result_str, load_densities, gedp,
+	report_progress, NULL,
+	gedp->ged_gvp && gedp->ged_run_on_main_thread ? report_plot : NULL,
+	&plot_state};
+    return analyze_gqa(&context, argc, argv) == ANALYZE_OK ?
+        BRLCAD_OK : BRLCAD_ERROR;
+}
+
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C++
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */
@@ -281,7 +281,7 @@ crofton_on_obj(struct db_i *dbip, const char *obj_name, size_t n_rays,
     * default).  A fixed ray count can be requested for reproducibility
     * or timing studies, but for reliable SA/volume estimates the
     * convergence path is strongly preferred.                           */
-    struct rt_crofton_params crp;
+    struct rt_crofton_params crp = {};
    if (n_rays > 0) {
 	crp.n_rays       = n_rays;
 	crp.stability_mm = 0.0;
@@ -6,6 +6,24 @@ brlcad_find_package(OPENNURBS REQUIRED SYSPATTERN openNURBS)
 brlcad_find_package(REGEX REQUIRED)
 find_package_gte(REQUIRED)

+set(OPENQMC_LIBS)
+if(BRLCAD_ENABLE_QMC)
+  brlcad_find_package(OpenQMC CONFIG)
+  if(OpenQMC_FOUND AND TARGET OpenQMC::OpenQMC)
+    add_definitions(-DBRLCAD_ENABLE_QMC)
+    set(OPENQMC_LIBS OpenQMC::OpenQMC)
+  else()
+    message(
+      WARNING
+      "QMC requested, but OpenQMC was not found - disabling"
+    )
+    set(
+      BRLCAD_ENABLE_QMC OFF CACHE BOOL
+      "Disabled because OpenQMC was not found" FORCE
+    )
+  endif()
+endif()
+
 if(BRLCAD_ENABLE_OPENCL)
  add_definitions(-DUSE_OPENCL)
 endif(BRLCAD_ENABLE_OPENCL)
@@ -435,6 +453,7 @@ set(
  ${RPCRT_LIB}
  ${STDCXX_LIBRARIES}
  manifold::manifold
+  ${OPENQMC_LIBS}
 )

 # Include directories needed by librt users
@@ -2450,7 +2450,11 @@ rt_ehy_surf_area(fastf_t *area, const struct rt_db_internal *ip)
     * an elliptic integral with no elementary closed form -- use Crofton.
     */
    if (!NEAR_EQUAL(eip->ehy_r1, eip->ehy_r2, RT_LEN_TOL)) {
-	do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(area, NULL, ip, &_p); } while (0);
+	do {
+	    static const struct rt_crofton_params params = {
+		50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	    rt_crofton_sample(area, NULL, ip, &params);
+	} while (0);
 	return;
    }

@@ -2098,7 +2098,11 @@ rt_ell_surf_area(fastf_t *area, const struct rt_db_internal *ip)
    default:
 	/* General triaxial ellipsoid: no closed-form solution exists.
 	 * Fall back to the Cauchy-Crofton ray-sampling estimator. */
-	do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(area, NULL, ip, &_p); } while (0);
+	do {
+	    static const struct rt_crofton_params params = {
+		50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	    rt_crofton_sample(area, NULL, ip, &params);
+	} while (0);
 	break;
    }
 }
@@ -2158,7 +2158,11 @@ rt_eto_volume(fastf_t *vol, const struct rt_db_internal *ip)
    if (eto_is_self_intersecting(tip)) {
 	struct rt_db_internal ip_meth = *ip;
 	ip_meth.idb_meth = &OBJ[ID_ETO];
-	do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(NULL, vol, &ip_meth, &_p); } while (0);
+	do {
+	    static const struct rt_crofton_params params = {
+		50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	    rt_crofton_sample(NULL, vol, &ip_meth, &params);
+	} while (0);
 	return;
    }

@@ -2189,7 +2193,11 @@ rt_eto_surf_area(fastf_t *area, const struct rt_db_internal *ip)
    if (eto_is_self_intersecting(tip)) {
 	struct rt_db_internal ip_meth = *ip;
 	ip_meth.idb_meth = &OBJ[ID_ETO];
-	do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(area, NULL, &ip_meth, &_p); } while (0);
+	do {
+	    static const struct rt_crofton_params params = {
+		50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	    rt_crofton_sample(area, NULL, &ip_meth, &params);
+	} while (0);
 	return;
    }

@@ -1702,7 +1702,11 @@ rt_hyp_surf_area(fastf_t *area, const struct rt_db_internal *ip)
 {
    if (!area || !ip)
 	return;
-    do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(area, NULL, ip, &_p); } while (0);
+    do {
+	static const struct rt_crofton_params params = {
+	    50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	rt_crofton_sample(area, NULL, ip, &params);
+    } while (0);
 }


@@ -69,6 +69,28 @@ struct submodel_specific {
 #define RT_SUBMODEL_SPECIFIC_MAGIC 0x73756253	/* subS */
 #define RT_CK_SUBMODEL_SPECIFIC(_p) BU_CKMAG(_p, RT_SUBMODEL_SPECIFIC_MAGIC, "submodel_specific")

+static int
+rt_submodel_is_unbounded(struct rt_i *rtip)
+{
+    struct region *regp;
+
+    for (BU_LIST_FOR(regp, region, &rtip->HeadRegion)) {
+	point_t region_min, region_max;
+
+	if (!regp->reg_treetop)
+	    continue;
+	if (rt_bound_tree(regp->reg_treetop, region_min, region_max) < 0)
+	    return 1;
+	for (int axis = 0; axis < 3; axis++) {
+	    if (INVALID(region_min[axis]) || INVALID(region_max[axis]))
+		return 1;
+	}
+    }
+
+    return 0;
+}
+
+
 static int
 rt_submodel_is_absolute_path(const char *path)
 {
@@ -306,15 +328,24 @@ rt_submodel_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rti
    bn_mat_inv(submodel->m2subm, sip->root2leaf);
    submodel->rtip = sub_rtip;

-    /* Propagate submodel bounding box back upwards, rotated&scaled. */
-    bg_rotate_bbox(stp->st_min, stp->st_max,
-		   submodel->subm2m,
-		   sub_rtip->mdl_min, sub_rtip->mdl_max);
+    /* mdl_min and mdl_max intentionally exclude infinite solids.  Preserve
+     * that distinction across the submodel boundary so a parent Boolean can
+     * either clip the submodel or recognize that its result is unbounded. */
+    if (rt_submodel_is_unbounded(sub_rtip)) {
+	VSETALL(stp->st_min, -INFINITY);
+	VSETALL(stp->st_max, INFINITY);
+	VSETALL(stp->st_center, 0.0);
+	stp->st_aradius = stp->st_bradius = INFINITY;
+    } else {
+	bg_rotate_bbox(stp->st_min, stp->st_max,
+		       submodel->subm2m,
+		       sub_rtip->mdl_min, sub_rtip->mdl_max);

-    VSUB2(diam, stp->st_max, stp->st_min);
-    VADD2SCALE(stp->st_center, stp->st_min, stp->st_max, 0.5);
-    VSCALE(radvec, diam, 0.5);
-    stp->st_aradius = stp->st_bradius = MAGNITUDE(radvec);
+	VSUB2(diam, stp->st_max, stp->st_min);
+	VADD2SCALE(stp->st_center, stp->st_min, stp->st_max, 0.5);
+	VSCALE(radvec, diam, 0.5);
+	stp->st_aradius = stp->st_bradius = MAGNITUDE(radvec);
+    }

    if (RT_G_DEBUG & (RT_DEBUG_DB|RT_DEBUG_SOLIDS)) {
 	bu_log("rt_submodel_prep(%s): finished loading database %s\n",
@@ -1245,7 +1245,11 @@ rt_superell_volume(fastf_t *volume, const struct rt_db_internal *ip)
    *volume = 2.0 * mag_a * mag_b * mag_c * sip->e * sip->n * (tgamma(sip->n/2.0 + 1.0) * tgamma(sip->n) / tgamma(3.0 * sip->n/2.0 + 1.0)) * (tgamma(sip->e / 2.0) * tgamma(sip->e / 2.0) / tgamma(sip->e));
 #else
    /* tgamma unavailable: fall back to Cauchy-Crofton estimation */
-    do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(NULL, volume, ip, &_p); } while (0);
+    do {
+	static const struct rt_crofton_params params = {
+	    50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	rt_crofton_sample(NULL, volume, ip, &params);
+    } while (0);
 #endif
 }

@@ -4010,7 +4010,11 @@ rt_tgc_surf_area(fastf_t *area, const struct rt_db_internal *ip)
 	    /* No closed-form formula exists for the general truncated
 	     * elliptic cone.  Fall back to the Cauchy-Crofton ray-sampling
 	     * estimator, which handles any TGC shape correctly. */
-	    do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(area, NULL, ip, &_p); } while (0);
+	    do {
+		static const struct rt_crofton_params params = {
+		    50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+		rt_crofton_sample(area, NULL, ip, &params);
+	    } while (0);
 	    break;
    }
 }
@@ -2021,7 +2021,11 @@ rt_tor_surf_area(fastf_t *area, const struct rt_db_internal *ip)
    if (tip->r_h > tip->r_a) {
 	struct rt_db_internal ip_meth = *ip;
 	ip_meth.idb_meth = &OBJ[ID_TOR];
-	do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(area, NULL, &ip_meth, &_p); } while (0);
+	do {
+	    static const struct rt_crofton_params params = {
+		50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	    rt_crofton_sample(area, NULL, &ip_meth, &params);
+	} while (0);
 	return;
    }

@@ -2042,7 +2046,11 @@ rt_tor_volume(fastf_t *vol, const struct rt_db_internal *ip)
    if (tip->r_h > tip->r_a) {
 	struct rt_db_internal ip_meth = *ip;
 	ip_meth.idb_meth = &OBJ[ID_TOR];
-	do { static const struct rt_crofton_params _p = {50000u, 0.0, 0.0}; rt_crofton_sample(NULL, vol, &ip_meth, &_p); } while (0);
+	do {
+	    static const struct rt_crofton_params params = {
+		50000u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	    rt_crofton_sample(NULL, vol, &ip_meth, &params);
+	} while (0);
 	return;
    }

@@ -34,6 +34,188 @@
 #include "wdb.h"


+struct region_visit_counts {
+    size_t segments;
+    double chord;
+};
+
+struct ray_visit_counts {
+    size_t rays;
+    size_t segments;
+    int valid;
+};
+
+struct progress_counts {
+    size_t calls;
+    size_t last_rays;
+    int valid;
+};
+
+static void
+visit_progress(size_t ray_count, size_t UNUSED(crossing_count),
+    double UNUSED(surface_area), double UNUSED(volume),
+    double UNUSED(stability_mm), int UNUSED(stability_evaluated),
+    double elapsed_ms, void *data)
+{
+    struct progress_counts *counts = (struct progress_counts *)data;
+    counts->calls++;
+    counts->last_rays = ray_count;
+    if (ray_count == 0 || elapsed_ms < 0.0)
+        counts->valid = 0;
+}
+
+static void
+visit_ray(const struct rt_crofton_ray *ray, void *data)
+{
+    struct ray_visit_counts *counts = (struct ray_visit_counts *)data;
+    counts->rays++;
+    for (size_t i = 0; i < ray->segment_count; i++) {
+	const struct rt_crofton_segment *segment = &ray->segments[i];
+	counts->segments++;
+	if (segment->ray_id != ray->ray_id ||
+	    !VNEAR_EQUAL(segment->ray_origin, ray->origin, SMALL_FASTF) ||
+	    !VNEAR_EQUAL(segment->ray_direction, ray->direction, SMALL_FASTF) ||
+	    !NEAR_EQUAL(segment->out_distance - segment->in_distance,
+		segment->thickness, SMALL_FASTF))
+	    counts->valid = 0;
+    }
+}
+
+static void
+visit_region_segment(const struct rt_crofton_segment *segment, void *data)
+{
+    struct region_visit_counts *counts = (struct region_visit_counts *)data;
+    const struct region *region = segment->region;
+    if (!region || !region->reg_name) return;
+    counts->segments++;
+    counts->chord += segment->thickness;
+}
+
+static int
+check_seeded_streams(struct rt_i *rtip, const char *label,
+    enum rt_crofton_sequence sequence)
+{
+    struct progress_counts progress = {0, 0, 1};
+    struct rt_crofton_params params = {
+        2048u, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, visit_progress, &progress};
+    struct rt_crofton_stats first_stats, repeated_stats, other_stats;
+    struct rt_crofton_stats session_stats = {0};
+    struct region_visit_counts first_counts = {0, 0.0};
+    struct region_visit_counts repeated_counts = {0, 0.0};
+    struct region_visit_counts other_counts = {0, 0.0};
+    struct region_visit_counts session_counts = {0, 0.0};
+    int first_ret = rt_crofton_visit_seeded(&first_stats, rtip,
+	&params, 0, NULL, NULL, sequence,
+	1234u, 7u, visit_region_segment, NULL, &first_counts);
+    int repeated_ret = rt_crofton_visit_seeded(&repeated_stats, rtip,
+	&params, 0, NULL, NULL, sequence,
+	1234u, 7u, visit_region_segment, NULL, &repeated_counts);
+    int other_ret = rt_crofton_visit_seeded(&other_stats, rtip,
+	&params, 0, NULL, NULL, sequence,
+	1234u, 8u, visit_region_segment, NULL, &other_counts);
+    struct rt_crofton_session *session = rt_crofton_session_create(rtip);
+    int session_ret = session ?
+        rt_crofton_session_visit_seeded(session, &session_stats, &params,
+            0, NULL, NULL, sequence, 1234u, 7u, visit_region_segment,
+            NULL, &session_counts) : -1;
+    rt_crofton_session_destroy(session);
+    const int session_match = session_ret == first_ret &&
+        session_stats.ray_count == first_stats.ray_count &&
+        session_stats.crossing_count == first_stats.crossing_count &&
+        EQUAL(session_stats.volume, first_stats.volume) &&
+        EQUAL(session_stats.surface_area, first_stats.surface_area) &&
+        session_counts.segments == first_counts.segments &&
+        EQUAL(session_counts.chord, first_counts.chord);
+    const int repeated_match = first_ret == repeated_ret &&
+	first_stats.ray_count == repeated_stats.ray_count &&
+	first_stats.crossing_count == repeated_stats.crossing_count &&
+	EQUAL(first_stats.volume, repeated_stats.volume) &&
+	EQUAL(first_stats.surface_area, repeated_stats.surface_area) &&
+	first_counts.segments == repeated_counts.segments &&
+	EQUAL(first_counts.chord, repeated_counts.chord);
+    const int independent_stream = other_ret >= 0 &&
+	(!EQUAL(first_stats.volume, other_stats.volume) ||
+	 !EQUAL(first_stats.surface_area, other_stats.surface_area));
+    const int stopped_at_limit =
+        first_stats.stop_reason == RT_CROFTON_STOP_RAYS &&
+        repeated_stats.stop_reason == RT_CROFTON_STOP_RAYS &&
+        other_stats.stop_reason == RT_CROFTON_STOP_RAYS &&
+        session_stats.stop_reason == RT_CROFTON_STOP_RAYS;
+    if (first_ret < 0 || !session_match || !repeated_match ||
+        !independent_stream || !stopped_at_limit || progress.calls != 4 ||
+        progress.last_rays != params.n_rays || !progress.valid) {
+	printf("  %-24s  invalid %s seeded Crofton streams\n", label,
+	    sequence == RT_CROFTON_SEQUENCE_QMC ? "QMC" : "random");
+	return 1;
+    }
+    return 0;
+}
+
+static int
+check_invalid_inputs(struct rt_i *rtip, const char *label)
+{
+    struct rt_crofton_stats stats = {0};
+    struct region_visit_counts counts = {0, 0.0};
+    struct rt_crofton_params params = {
+        1u, NAN, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+    point_t bbox_min, bbox_max;
+    VSETALL(bbox_min, 1.0);
+    VSETALL(bbox_max, 0.0);
+
+    const int nonfinite = rt_crofton_visit(&stats, rtip, &params,
+        0, NULL, NULL, RT_CROFTON_SEQUENCE_RANDOM,
+        visit_region_segment, NULL, &counts);
+    params.stability_mm = 0.0;
+    const int missing_bound = rt_crofton_visit(&stats, rtip, &params,
+        0, bbox_min, NULL, RT_CROFTON_SEQUENCE_RANDOM,
+        visit_region_segment, NULL, &counts);
+    const int reversed_bounds = rt_crofton_visit(&stats, rtip, &params,
+        0, bbox_min, bbox_max, RT_CROFTON_SEQUENCE_RANDOM,
+        visit_region_segment, NULL, &counts);
+    params.stability_metrics = 1u << 31;
+    const int unknown_metric = rt_crofton_visit(&stats, rtip, &params,
+        0, NULL, NULL, RT_CROFTON_SEQUENCE_RANDOM,
+        visit_region_segment, NULL, &counts);
+    if (nonfinite != -1 || missing_bound != -1 || reversed_bounds != -1 ||
+        unknown_metric != -1) {
+        printf("  %-24s  accepted invalid Crofton controls\n", label);
+        return 1;
+    }
+    return 0;
+}
+
+static int
+check_stopping_controls(struct rt_i *rtip, const char *label)
+{
+    const size_t expected_stability_rays = 80000u;
+    struct region_visit_counts counts = {0, 0.0};
+    struct rt_crofton_stats stats = {0};
+    struct rt_crofton_params params = {
+        100000u, 0.0, 10000.0, RT_CROFTON_STABILITY_VOLUME, NULL, NULL};
+    int ret = rt_crofton_visit(&stats, rtip, &params, 0, NULL, NULL,
+        RT_CROFTON_SEQUENCE_RANDOM, visit_region_segment, NULL, &counts);
+    if (ret <= 0 || stats.stop_reason != RT_CROFTON_STOP_RAYS ||
+        stats.ray_count != params.n_rays) {
+        printf("  %-24s  time-only run did not honor its ray limit\n", label);
+        return 1;
+    }
+
+    params.n_rays = 0;
+    params.stability_mm = MAX_FASTF;
+    params.time_ms = 1000.0;
+    counts.segments = 0;
+    counts.chord = 0.0;
+    ret = rt_crofton_visit(&stats, rtip, &params, 0, NULL, NULL,
+        RT_CROFTON_SEQUENCE_RANDOM, visit_region_segment, NULL, &counts);
+    if (ret <= 0 || stats.stop_reason != RT_CROFTON_STOP_STABILITY ||
+        stats.ray_count != expected_stability_rays) {
+        printf("  %-24s  invalid doubling-checkpoint stop (%zu rays)\n",
+            label, stats.ray_count);
+        return 1;
+    }
+    return check_invalid_inputs(rtip, label);
+}
+
 static double
 rel_err(double estimated, double exact)
 {
@@ -43,12 +225,88 @@ rel_err(double estimated, double exact)
 }


+static int
+check_halfspace_clipping(void)
+{
+    const fastf_t cube_extent = 1000.0;
+    const fastf_t split_offset = 0.5 * cube_extent;
+    const size_t ray_count = 20000u;
+    const double exact_volume = 500000000.0;
+    const double maximum_relative_error = 0.05;
+    struct db_i *dbip;
+    struct rt_wdb *wdbp;
+    struct rt_i *rtip = NULL;
+    point_t minimum = VINIT_ZERO;
+    point_t maximum = {cube_extent, cube_extent, cube_extent};
+    vect_t normal = {1.0, 0.0, 0.0};
+    struct wmember members;
+    double volume = 0.0;
+    int failures = 0;
+
+    dbip = db_create_inmem();
+    if (!dbip) return 1;
+    wdbp = wdb_dbopen(dbip, RT_WDB_TYPE_DB_INMEM);
+    if (!wdbp) {
+	db_close(dbip);
+	return 1;
+    }
+
+    if (mk_rpp(wdbp, "half_cube.s", minimum, maximum) ||
+	mk_half(wdbp, "clip.s", normal, split_offset)) {
+	failures++;
+	goto done;
+    }
+    BU_LIST_INIT(&members.l);
+    if (!mk_addmember("half_cube.s", &members.l, NULL, WMOP_UNION) ||
+	!mk_addmember("clip.s", &members.l, NULL, WMOP_INTERSECT)) {
+	mk_freemembers(&members.l);
+	failures++;
+	goto done;
+    }
+    if (mk_lcomb(wdbp, "half_cube.r", &members, 1, NULL, NULL, NULL, 0)) {
+	failures++;
+	goto done;
+    }
+
+    db_update_nref(dbip);
+    rtip = rt_i_create(dbip);
+    if (!rtip || rt_gettree(rtip, "half_cube.r")) {
+	failures++;
+	goto done;
+    }
+    rt_prep_parallel(rtip, 1);
+    {
+	struct rt_crofton_params params = {
+	    ray_count, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
+	const int result = rt_crofton_shoot(NULL, &volume, NULL, NULL, NULL,
+	    NULL, NULL, rtip, &params, NULL, NULL);
+	if (result < 1 ||
+	    rel_err(volume, exact_volume) > maximum_relative_error)
+	    failures++;
+    }
+
+done:
+    if (rtip) rt_i_destroy(rtip);
+    db_close(dbip);
+    if (failures)
+	printf("  half-space clipping       Crofton estimate failed\n");
+    return failures;
+}
+
+
 static int
 crofton_segments_equal(const struct rt_crofton_segment *left,
 		       const struct rt_crofton_segment *right)
 {
    return left->ray_id == right->ray_id &&
+	left->region == right->region &&
+	left->in_solid == right->in_solid &&
+	left->out_solid == right->out_solid &&
+	NEAR_EQUAL(left->in_distance, right->in_distance, SMALL_FASTF) &&
+	NEAR_EQUAL(left->out_distance, right->out_distance, SMALL_FASTF) &&
 	NEAR_EQUAL(left->thickness, right->thickness, SMALL_FASTF) &&
+	VNEAR_EQUAL(left->ray_origin, right->ray_origin, SMALL_FASTF) &&
+	VNEAR_EQUAL(left->ray_direction, right->ray_direction, SMALL_FASTF) &&
 	VNEAR_EQUAL(left->in_point, right->in_point, SMALL_FASTF) &&
 	VNEAR_EQUAL(left->in_normal, right->in_normal, SMALL_FASTF) &&
 	VNEAR_EQUAL(left->out_point, right->out_point, SMALL_FASTF) &&
@@ -63,7 +321,7 @@ verify_crofton_estimates(void)

    printf("\n--- Crofton estimator verification ---\n");

-    struct rt_crofton_params cparams = {0, 0.0, 0.0};
+    struct rt_crofton_params cparams = {0, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};

 #define CROFTON_CHECK(label, ip_ptr, analytic_sa, analytic_vol) \
    do { \
@@ -378,7 +636,7 @@ run_convergence_case(struct db_i *dbip,
 	rt_prep_parallel(rtip, 1);

 	/* Convergence-based stopping with a per-target wall-clock budget. */
-	struct rt_crofton_params p = { 0u, 0.05, 1000.0 };
+	struct rt_crofton_params p = { 0u, 0.05, 1000.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL };
 	if (target_pct <= 7.0)
 	    p.time_ms = 1500.0;
 	if (target_pct <= 5.0)
@@ -397,7 +655,7 @@ run_convergence_case(struct db_i *dbip,
 	    const size_t ray_offset = 4096u;
 	    const size_t continuation_rays = 256u;
 	    struct rt_crofton_params sample_params =
-		{ray_offset + continuation_rays, 0.0, 0.0};
+		{ray_offset + continuation_rays, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
 	    int sample_ret = rt_crofton_collect(&samples, rtip,
 		&sample_params, 0, NULL, NULL);
 	    int full_samples_valid = sample_ret > 0 && samples.segments &&
@@ -426,10 +684,43 @@ run_convergence_case(struct db_i *dbip,
 		}
 	    }

+	    {
+		struct rt_crofton_stats stats;
+		struct region_visit_counts counts = {0, 0.0};
+		int visit_ret = rt_crofton_visit(&stats, rtip, &sample_params,
+		    0, NULL, NULL, RT_CROFTON_SEQUENCE_RANDOM,
+		    visit_region_segment, NULL, &counts);
+		if (visit_ret <= 0 || counts.segments * 2 != stats.crossing_count ||
+		    stats.ray_count != sample_params.n_rays ||
+		    counts.chord <= 0.0 || stats.volume <= 0.0) {
+		    printf("  %-24s  invalid streaming Crofton output\n", label);
+		    failures++;
+		}
+
+		struct ray_visit_counts ray_counts = {0, 0, 1};
+		visit_ret = rt_crofton_visit_rays(&stats, rtip, &sample_params,
+		    0, NULL, NULL, RT_CROFTON_SEQUENCE_RANDOM,
+		    visit_ray, NULL, &ray_counts);
+		if (visit_ret <= 0 || ray_counts.rays != sample_params.n_rays ||
+		    ray_counts.segments * 2 != stats.crossing_count ||
+		    !ray_counts.valid) {
+		    printf("  %-24s  invalid ray streaming Crofton output\n",
+			label);
+		    failures++;
+		}
+	    }
+
+	    failures += check_seeded_streams(rtip, label,
+		RT_CROFTON_SEQUENCE_RANDOM);
+	    if (rt_crofton_qmc_available())
+		failures += check_seeded_streams(rtip, label,
+		    RT_CROFTON_SEQUENCE_QMC);
+	    failures += check_stopping_controls(rtip, label);
+
 	    struct rt_crofton_result offset_samples =
 		RT_CROFTON_RESULT_INIT;
 	    struct rt_crofton_params offset_params =
-		{continuation_rays, 0.0, 0.0};
+		{continuation_rays, 0.0, 0.0, RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
 	    sample_ret = rt_crofton_collect(&offset_samples, rtip,
 		&offset_params, ray_offset, NULL, NULL);
 	    int offset_samples_valid = sample_ret > 0 &&
@@ -571,6 +862,7 @@ main(int argc, char *argv[])

    int failures = 0;
    failures += verify_crofton_estimates();
+    failures += check_halfspace_clipping();
    failures += test_crofton_convergence_timing();

    printf("\n=== Summary: %d failure(s) ===\n", failures);
@@ -321,7 +321,8 @@ compare_paths(const char  *label,

    /* fixed ray count for deterministic regression behavior. */
    static const struct rt_crofton_params acc_p =
-        { CROFTON_ACCURACY_RAYS, 0.0, 0.0 };
+        {CROFTON_ACCURACY_RAYS, 0.0, 0.0,
+	 RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
    struct crofton_result dda_acc = run_crofton(rtip, &acc_p);

    printf("\n    %-12s  %14s  %14s  %10s  %10s\n",
@@ -4273,7 +4273,9 @@ scan_input_g(const char *g_path, const char *g_root)
 	    if (cr_rtip) {
 		if (rt_gettree(cr_rtip, dp->d_namep) == 0) {
 		    rt_prep_parallel(cr_rtip, 1);
-		    struct rt_crofton_params crp = { 2000u, 0.0, 0.0 };
+		    struct rt_crofton_params crp = {
+			2000u, 0.0, 0.0,
+			RT_CROFTON_STABILITY_DEFAULT, NULL, NULL};
 		    cr = rt_crofton_shoot(&csa, &cv, NULL, NULL, NULL, NULL,
 			NULL, cr_rtip, &crp, NULL, NULL);
 		}
@@ -426,6 +426,8 @@ mged_attach(struct mged_state *s, const char *wp_name, int argc, const char *arg
    int opt_argc;
    char **opt_argv;
    struct mged_dm *o_dm;
+    const char *dm_name;
+    const char *dm_lname;

    if (!wp_name) {
 	return TCL_ERROR;
@@ -495,8 +497,8 @@ mged_attach(struct mged_state *s, const char *wp_name, int argc, const char *arg
    mged_link_vars(s->mged_curr_dm);

    Tcl_ResetResult(s->interp);
-    const char *dm_name = dm_get_dm_name(DMP);
-    const char *dm_lname = dm_get_dm_lname(DMP);
+    dm_name = dm_get_dm_name(DMP);
+    dm_lname = dm_get_dm_lname(DMP);
    if (dm_name && dm_lname) {
 	Tcl_AppendResult(s->interp, "ATTACHING ", dm_name, " (", dm_lname,	")\n", (char *)NULL);
    }
@@ -1099,7 +1099,7 @@ cmd_ged_gqa(ClientData clientData, Tcl_Interp *interpreter, int argc, const char
    int i;
    int option;
    int ret;
-    bool plot_in_view = false;
+    bool legacy_plot_in_view = false;
    struct cmdtab *ctp = (struct cmdtab *)clientData;
    MGED_CK_CMD(ctp);
    struct mged_state *s = ctp->s;
@@ -1125,7 +1125,7 @@ cmd_ged_gqa(ClientData clientData, Tcl_Interp *interpreter, int argc, const char
    bu_optind = 1;
    while ((option = bu_getopt(argc, (char * const *)argv, gqa_options)) != -1) {
 	if (option == 'A' && bu_optarg && strchr(bu_optarg, 'p'))
-	    plot_in_view = true;
+	    legacy_plot_in_view = true;
    }
    i = bu_optind;
    bu_optarg = saved_optarg;
@@ -1157,10 +1157,10 @@ cmd_ged_gqa(ClientData clientData, Tcl_Interp *interpreter, int argc, const char
 	gd_rt_cmd_len += ged_who_argv(s->gedp, vp, (const char **)&gd_rt_cmd[args]);
    }

-    /* In-view plotting creates display objects and invokes MGED's display
-     * manager callbacks.  Keep that work on the GUI thread while preserving
-     * asynchronous execution for analyses that only produce reports/files. */
-    if (plot_in_view)
+    /* Legacy -A p mutates view objects directly and must remain on the GUI
+     * thread.  The experimental --plot-view path uses GED's main-thread
+     * callback and may run asynchronously with the other analyses. */
+    if (legacy_plot_in_view)
 	ret = (*ctp->ged_func)(s->gedp, gd_rt_cmd_len, (const char **)gd_rt_cmd);
    else
 	ret = run_ged_async(s, [&]() -> int { return (*ctp->ged_func)(s->gedp, gd_rt_cmd_len, (const char **)gd_rt_cmd); });
@@ -87,6 +87,7 @@ f_facedef(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[
    struct rt_arb_internal *arbo;
    plane_t planes[6];
    int status;
+    int arb_type;
    struct bu_vls error_msg;

    if (argc < 2) {
@@ -125,7 +126,7 @@ f_facedef(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[
    RT_ARB_CK_MAGIC(arb);

    /* find new planes to account for any editing */
-    int arb_type = rt_arb_std_type(&MEDIT(s)->es_int, MEDIT(s)->tol);
+    arb_type = rt_arb_std_type(&MEDIT(s)->es_int, MEDIT(s)->tol);
    if (rt_arb_calc_planes(&error_msg, arb, arb_type, planes, &s->tol.tol)) {
 	Tcl_AppendResult(interp, bu_vls_addr(&error_msg),
 			 "Unable to determine plane equations\n", (char *)NULL);
@@ -113,9 +113,6 @@ extern struct _mged_variables default_mged_variables;
 /* defined in color_scheme.c */
 extern struct _color_scheme default_color_scheme;

-/* defined in grid.c */
-extern struct bv_grid_state default_grid_state;
-
 /* defined in axes.c */
 extern struct _axes_state default_axes_state;

@@ -77,6 +77,10 @@
 /* Needed to define struct bv_scene_obj */
 #include "bv/defines.h"

+__BEGIN_DECLS
+extern struct bv_grid_state default_grid_state;
+__END_DECLS
+
 __BEGIN_DECLS
 #include "tcl.h"
 #ifdef HAVE_TK
@@ -602,6 +602,15 @@ mged_refresh_handler(void *clientdata)
    mged_run_on_gui_thread(s, mged_refresh_handler_impl, clientdata);
 }

+static void
+mged_main_thread_handler(void *clientdata, ged_main_thread_callback_t callback,
+	void *data)
+{
+    struct mged_state *s = (struct mged_state *)clientdata;
+    MGED_CK_STATE(s);
+    mged_run_on_gui_thread(s, callback, data);
+}
+
 /*
 * Initialize mged, configure the path, set up the tcl interpreter.
 */
@@ -660,6 +669,8 @@ mged_setup(struct mged_state *s)
    s->gedp->ged_output_handler = mged_output_handler;
    s->gedp->ged_refresh_clientdata = (void *)s;
    s->gedp->ged_refresh_handler = mged_refresh_handler;
+    s->gedp->ged_main_thread_clientdata = (void *)s;
+    s->gedp->ged_run_on_main_thread = mged_main_thread_handler;
    s->gedp->vlist_ctx = (void *)s;
    s->gedp->ged_create_vlist_scene_obj_callback = createDListSolid;
    s->gedp->ged_create_vlist_display_list_callback = createDListAll;
@@ -25,13 +25,16 @@

 #include "common.h"

-#include <vector>
-#include <set>
-#include <string>
+#include <cerrno>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 #include <iostream>
 #include <limits>
-#include <cstdio>
-#include <cstring>
+#include <set>
+#include <string>
+#include <vector>

 /* needed on mac in c90 mode */
 #ifndef HAVE_DECL_FSEEKO
@@ -50,6 +53,7 @@ extern "C" b_off_t ftello(FILE *);
 #include "bu/process.h"
 #include "bu/units.h"
 #include "analyze.h"
+#include "../libbu/json.hpp"

 #include "linenoise.hpp"

@@ -422,6 +426,18 @@ nirt_app_exec(struct nirt_state *ns, struct bu_vls *iline, struct bu_vls *state_
    return nret;
 }

+static bool
+parse_replay_id(const char *value, unsigned long long &id)
+{
+    if (!value || !value[0] || value[0] == '-')
+	return false;
+
+    char *end = NULL;
+    errno = 0;
+    id = std::strtoull(value, &end, 10);
+    return errno == 0 && end != value && *end == '\0';
+}
+

 int
 main(int argc, const char **argv)
@@ -453,6 +469,12 @@ main(int argc, const char **argv)
    struct db_i *dbip;
    struct nirt_io_data io_data = IO_DATA_NULL;
    struct nirt_state *ns = NULL;
+    std::string replay_file;
+    unsigned long long replay_id = 0;
+    bool replay_requested = false;
+    bool replay_id_set = false;
+    std::vector<const char *> filtered_argv;
+    fastf_t replay_base2local = 1.0;

    /* These bu_opt_desc_opts settings approximate the old struct nirt_state help formatting */
    struct bu_opt_desc_opts dopts = { BU_OPT_ASCII, 1, 15, 65, NULL, NULL, NULL, 1, NULL, NULL };
@@ -478,9 +500,51 @@ main(int argc, const char **argv)
 	bu_vls_printf(&launch_cmd, "%s", argv[argc-1]);
    }

+    /* Extract replay controls before libbu parses the normal nirt options.  They
+     * are application-level controls and should not alter the established
+     * nirt option structure or its existing command line behavior. */
+    filtered_argv.push_back(argv[0]);
+    for (int ic = 1; ic < argc; ic++) {
+	const char *arg = argv[ic];
+	if (BU_STR_EQUAL(arg, "--replay")) {
+	    if (ic + 1 >= argc)
+		bu_exit(EXIT_FAILURE,
+		    "--replay requires a JSON file\n");
+	    replay_file = argv[++ic];
+	    replay_requested = true;
+	} else if (BU_STR_EQUAL(arg, "--ray-id")) {
+	    if (ic + 1 >= argc)
+		bu_exit(EXIT_FAILURE,
+		    "--ray-id requires an integer\n");
+	    if (!parse_replay_id(argv[++ic], replay_id))
+		bu_exit(EXIT_FAILURE,
+		    "--ray-id requires a nonnegative integer\n");
+	    replay_requested = true;
+	    replay_id_set = true;
+	} else if (bu_strncmp(arg, "--replay=", 9) == 0) {
+	    replay_file = arg + 9;
+	    replay_requested = true;
+	} else if (bu_strncmp(arg, "--ray-id=", 9) == 0) {
+	    if (!parse_replay_id(arg + 9, replay_id))
+		bu_exit(EXIT_FAILURE,
+		    "--ray-id requires a nonnegative integer\n");
+	    replay_requested = true;
+	    replay_id_set = true;
+	} else {
+	    filtered_argv.push_back(arg);
+	}
+    }
+
+    if (replay_requested && replay_file.empty())
+	bu_exit(EXIT_FAILURE, "--ray-id requires --replay FILE\n");
+    if (replay_requested && !replay_id_set)
+	bu_exit(EXIT_FAILURE, "--replay requires --ray-id N\n");
+
    /* Let libbu know where we are */
    bu_setprogname(argv[0]);

+    argv = filtered_argv.data();
+    argc = (int)filtered_argv.size();
    argv++; argc--;

    ac = bu_opt_parse(&optparse_msg, argc, (const char **)argv, d);
@@ -720,6 +784,7 @@ main(int argc, const char **argv)
 	ret = EXIT_FAILURE;
 	goto done;
    }
+    replay_base2local = dbip->dbi_base2local;
    db_close(dbip); /* nirt will now manage its own copies of the dbip */

    /* Report Database info */
@@ -745,6 +810,106 @@ main(int argc, const char **argv)
 	}
    }

+    if (replay_requested) {
+	try {
+	    std::ifstream replay_stream(replay_file.c_str());
+	    if (!replay_stream) {
+		bu_vls_sprintf(&msg, "Unable to open replay file %s\n",
+		    replay_file.c_str());
+		nirt_err(&io_data, bu_vls_cstr(&msg));
+		ret = EXIT_FAILURE;
+		goto done;
+	    }
+	    nlohmann::json replay_json;
+	    replay_stream >> replay_json;
+	    if (!replay_json.is_object() ||
+		!replay_json.contains("schema_version") ||
+		!replay_json["schema_version"].is_number_unsigned() ||
+		(replay_json["schema_version"].get<unsigned int>() != 1 &&
+		 replay_json["schema_version"].get<unsigned int>() != 2) ||
+		!replay_json.contains("rays") || !replay_json["rays"].is_array()) {
+		nirt_err(&io_data,
+		    "Replay file must use schema version 1 or 2 and contain a rays array\n");
+		ret = EXIT_FAILURE;
+		goto done;
+	    }
+	    const nlohmann::json *ray = NULL;
+	    for (const auto &candidate : replay_json["rays"]) {
+		if (candidate.contains("ray_id") &&
+		    candidate["ray_id"].is_number_unsigned() &&
+		    candidate["ray_id"].get<unsigned long long>() == replay_id) {
+		    if (ray) {
+			bu_vls_sprintf(&msg,
+			    "Ray ID %llu occurs more than once in %s\n",
+			    replay_id, replay_file.c_str());
+			nirt_err(&io_data, bu_vls_cstr(&msg));
+			ret = EXIT_FAILURE;
+			goto done;
+		    }
+		    ray = &candidate;
+		}
+	    }
+	    if (!ray || !ray->contains("origin_mm") || !ray->contains("direction") ||
+		!(*ray)["origin_mm"].is_array() || !(*ray)["direction"].is_array() ||
+		(*ray)["origin_mm"].size() != 3 || (*ray)["direction"].size() != 3) {
+		bu_vls_sprintf(&msg,
+		    "Ray %llu was not found or is malformed in %s\n",
+		    replay_id, replay_file.c_str());
+		nirt_err(&io_data, bu_vls_cstr(&msg));
+		ret = EXIT_FAILURE;
+		goto done;
+	    }
+	    point_t origin;
+	    vect_t direction;
+	    for (int coord = 0; coord < 3; coord++) {
+		if (!(*ray)["origin_mm"][coord].is_number() ||
+		    !(*ray)["direction"][coord].is_number()) {
+		    nirt_err(&io_data, "Replay ray coordinates must be numeric\n");
+		    ret = EXIT_FAILURE;
+		    goto done;
+		}
+		origin[coord] = (*ray)["origin_mm"][coord].get<double>() *
+		    replay_base2local;
+		direction[coord] = (*ray)["direction"][coord].get<double>();
+		if (!std::isfinite(origin[coord]) ||
+		    !std::isfinite(direction[coord])) {
+		    nirt_err(&io_data,
+			"Replay ray contains non-finite coordinates\n");
+		    ret = EXIT_FAILURE;
+		    goto done;
+		}
+	    }
+	    const fastf_t direction_magnitude = MAGNITUDE(direction);
+	    if (!std::isfinite(direction_magnitude) ||
+		direction_magnitude <= VUNITIZE_TOL) {
+		nirt_err(&io_data, "Replay ray has a zero direction\n");
+		ret = EXIT_FAILURE;
+		goto done;
+	    }
+	    bu_vls_sprintf(&msg, "Replaying ray %llu from %s\n",
+		replay_id, replay_file.c_str());
+	    nirt_msg(&io_data, bu_vls_cstr(&msg));
+	    bu_vls_sprintf(&ncmd,
+		"backout 0; xyz %.17g %.17g %.17g; "
+		"dir %.17g %.17g %.17g; s",
+		origin[0], origin[1], origin[2],
+		direction[0], direction[1], direction[2]);
+	    if (nirt_exec(ns, bu_vls_cstr(&ncmd)) < 0) {
+		nirt_err(&io_data, "Replay shot failed\n");
+		ret = EXIT_FAILURE;
+		goto done;
+	    }
+	    ret = EXIT_SUCCESS;
+	    goto done;
+	} catch (const std::exception &e) {
+	    bu_vls_sprintf(&msg, "Unable to parse replay file %s: %s\n",
+		replay_file.c_str(), e.what());
+	    nirt_err(&io_data, bu_vls_cstr(&msg));
+	    ret = EXIT_FAILURE;
+	    goto done;
+	}
+    }
+
    /* If we're supposed to read matrix input from stdin instead of interacting, do that */
    if (optv.read_matrix) {
 	while ((buf = rt_read_cmd(stdin)) != (char *) 0) {
@@ -889,4 +1054,3 @@ done:
 // c-file-style: "stroustrup"
 // End:
 // ex: shiftwidth=4 tabstop=8
-
Author	SHA1	Message	Date
Clifford Yapp	dac8f504bf	Complete experimental GQA bounds and visualization Derive finite sampling bounds from prepared region trees so bounded half-space intersections and subtractions remain analyzable. Reject unbounded selections explicitly, propagate unbounded submodels, and keep infinite solids out of the Crofton sampling sphere without removing them from Boolean evaluation. Add an unsampled bounds measurement with text and JSON reporting. Use the same conservative bounds for sampling, adaptive dimensional targets, and automatic visualization scale. Restore Plot3 output for volume, gaps, overlaps, adjacent air, and exposed air. Provide exact all-rays and bounded representative policies. Representatives retain deterministic spatial samples, draw finite grid or Crofton support footprints, and preserve diagnostic severity and region identity. Publish immutable plot batches only between parallel ray batches. Add an MGED main-thread bridge that incrementally replaces or appends the named gqa::analysis overlay while keeping all view mutation on the display-owning thread. Repair named vlblock replacement so refreshed objects preserve their identity and view ownership. Document the behavior, limitations, and supporting sampling references, including Owen's replication guidance. Add the analytic and complex model comparison report. Expand regression coverage for bounds, half-spaces, submodels, plot policies, diagnostic files, and live MGED output.	2026-09-21 16:06:58 +00:00
Clifford Yapp	37f583285d	fixes	2026-09-21 06:33:46 -04:00
Clifford Yapp	ce1307152b	Add experimental GQA analysis modes Add a separate --analyze interface while preserving existing GQA behavior. Move the new measurement engine into libanalyze so the standalone and GED commands share bu_opt parsing, sampling, density loading, reporting, and schema-versioned JSON output. Provide axis-aligned and rotated three-view grids. The rotated mode builds an orthonormal view triad, clips boundary cells to projected model bounds, includes transverse covariance in mass moments, and shifts successive grids to reduce alignment aliasing. Run grid rays in parallel with worker-local accumulators and publish only complete three-view passes. Add time-first defaults, scale-aware starting spacing, dimensional stability checks, directional spread reporting, and progress updates. Retain explicit spacing, refinement, ray, seed, and deadline controls for repeatable studies. Extend Crofton sampling with selectable random and QMC sequences, streaming region and ray visitors, reusable parallel sessions, diagnostic overlap records, progress callbacks, and explicit stopping statistics. Estimate sampling intervals from independent randomized replications and keep pilot planning separate from production results. Add invalid-ray replay support to nirt, update the user and developer manuals, and expand regression coverage for parsing, clipped and rotated grids, uncertainty, diagnostics, output safety, and reproducibility. Reject malformed Crofton controls before sampling and make legacy overlap invalidation explicit in its report.	2026-09-20 12:02:59 +00:00
Clifford Yapp	7880489c2b	Fix MGED Unity build declarations	2026-09-20 12:02:59 +00:00
Clifford Yapp	55e788215a	Support using OpenQMC for Cauchy/Crofton sampling instead of PRNG	2026-09-20 12:02:59 +00:00