Files
Thomas Ales [M S E] a8b962e317 Add min peak frequency floor (CACH v4) and Batch Export View as Images
Two strands of in-progress work, committed together because they overlap
in sras_workers.py and main_window.py.

Min peak frequency floor:
- CACH tail bumped to version 4, adding u32 min_freq_khz provenance in
  fixed-point kHz (a float32 20.1 reads back as 20.10000038 and would
  report a spurious mismatch forever). v1-v3 tails read as no floor.
- Stored FFT caches are accepted when the reader's floor is at or above
  the stored one, since a higher floor is re-applicable by masking.
- Floor plumbed through compute_rf_image, BatchCacheWorker and the viewer.

Batch Export View as Images:
- New sras_render.py holds draw_view_image, shared by the Qt canvas and
  the headless exporter so a PNG cannot drift from what the GUI shows.
  Deliberately Qt-free so it is importable in a pool subprocess.
- BatchExportImagesWorker renders the current view settings across many
  files, process-pooled with an inline fallback, reporting per-file
  output names so the caller can flag same-stem collisions.
- _axes_extent extracted into sras_format for both render paths.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
2026-08-12 10:17:43 -05:00

428 lines
25 KiB
Markdown
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
# sras-viewer design notes
Rationale that outgrew code comments. Each section is referenced by a short
pointer comment at the relevant definition, so the code stays scannable and
the reasoning stays findable.
## Memory budget and row chunking (`sras_compute.py`)
DC images are computed over row chunks so the float32 working buffers for one
chunk stay under a memory budget. A fixed row count (the original design)
works fine for small legacy scans but is catastrophic for a v6 scan with a
large per-angle frame/sample count — e.g. a 7500-frame × 2500-sample angle
needs ~2.4 GB for a single 32-row chunk.
With chunks running concurrently the budget has to cover *all* live chunks at
once. On a large scan `chunk_rows` is already clamped to its floor of one row
(one row alone is ~75 MB of float32 at 7507×2500), so shrinking the per-chunk
size cannot buy more concurrency — the worker count must be derived from the
budget instead: `_plan_chunks` picks the worker count *first* and sizes the
chunk to it. Sizing the chunk first is the trap: a single chunk would always
consume the whole budget and leave room for exactly one worker, precisely on
the large scans that need concurrency most.
The 1024 MB default (`SRAS_MEM_BUDGET_MB`) is the measured knee on a 16-core
machine against a 7507-frame × 2500-sample angle: 512 MB left ~20% of the
speedup on the table, and 1536+ MB cost ~0.4 GB more resident memory for no
further gain.
A caller that itself runs several computations concurrently (angle-level
parallelism, `plan_angle_level`) must pass *both* `max_workers=1` and its
share of the budget. Capping the workers alone is not enough: the chunk would
still be sized against the whole budget, and N concurrent callers would each
allocate all of it.
## FFT peak search: block-parallel direct transform (`sras_compute.py`)
The displayed RF value per pixel is the argmax of the zero-padded power
spectrum of that pixel's CH1 waveform. This used to run through `_peak_bins_zoom`,
a coarse-rfft-plus-local-fine-DFT refinement that avoided ever materialising
a padded spectrum — at the pad factor of 40 needed for mapping resolution, a
full padded spectrum is ~9 GB per scan row, which used to collapse the old
row-chunk planner to one worker and make synthesis single-threaded. That
refinement was removed once pyFFTW became the sole, mandatory FFT backend
(SciPy dropped as a compute backend entirely): `_peak_bins` now always runs
the real full transform, and the memory problem zoom dodged is instead
solved by bounding *per-block* spectrum memory rather than avoiding full
spectra altogether.
`_peak_bins` runs the full transform (`_block_rfft`, a cached pyFFTW
`builders.rfft` plan, FFTW_MEASURE, wisdom persisted under
`~/.cache/sras-viewer/`) and argmaxes the power spectrum, in blocks fanned
out task-parallel over a persistent thread pool (`_fft_pool()`) — each pool
thread runs one single-threaded transform at a time, so aggregate
parallelism equals the pool's worker count. This block-streaming structure
is not just about parallelism: it is also what keeps memory bounded, by
never materialising more than one block's worth of full padded spectrum at
a time, regardless of how many waveforms a chunk holds.
The block size is the part that has to adapt to pad factor.
`_fft_block_for(spf, n_len)` derives waveforms-per-task from a fixed
per-thread byte budget (`_FFT_PLAN_BYTES_BUDGET`, `SRAS_FFT_PLAN_BUDGET_MB`,
default 16 MB) rather than a fixed constant, because a cached pyFFTW plan's
input+output buffers are *permanent* per-thread memory (the plan cache is
never evicted) — with a fixed 512-waveform block, pad 40 at a 2500-sample
frame costs ~210 MB per pool thread (~3.3 GB total across 16 threads);
`_fft_block_for` bounds that to ~16 MB per thread (~260 MB across 16
threads) at the same pad factor, while still reproducing the old tuned 512
exactly at natural resolution (pad 1), where it cost nothing to begin with.
`_FFT_BLOCK_MAX` (512) and `_FFT_BLOCK_MIN` (32) cap and floor the result:
the ceiling is the measured knee on a 16-core machine at natural resolution
(smaller blocks serialise on GIL-held numpy dispatch, larger ones lose cache
residency and task granularity); the floor keeps task granularity from
collapsing at extreme pad factors, at the cost of exceeding the byte budget
there.
The outer row-chunk sizing (`_plan_fft_rows`) needed no companion change.
It only ever budgets the raw float32 waveform *read* buffer, which this
change doesn't touch — spectrum memory is bounded independently by
`_fft_block_for`, and since the pool only ever runs as many blocks
concurrently as it has workers, peak transient spectrum memory during a
chunk's FFT phase is `_MAX_WORKERS * block * bytes_per_wf`, the same bound
whether the chunk holds 10 rows or 10,000. Queuing more rows into one chunk
to keep the read-row pool busy therefore can't blow up spectrum memory.
`tools/check_equivalence.py`'s golden-hash harness remains the end-to-end
regression baseline for this path, unaffected by this change.
## Row-averaged FFT: same-row, distance-weighted SNR cleanup (`sras_compute.py`)
`compute_rf_image`'s `row_avg_n` parameter averages each pixel's CH1
waveform with its up-to-n same-row neighbors before the FFT peak search, to
improve SNR on noisy scans. Never crosses rows: pixel pitch is strongly
anisotropic and varies by scan (5 µm × 50 µm on a typical scan, but as
stretched as 5 µm × 1 mm on others), so a physically meaningful "neighbor"
set can't be a fixed-shape 2-D window — but the X pitch *within one row* is
a single file-wide constant (`SrasFile.pixel_x_mm`), so restricting to the
row axis sidesteps the anisotropy question entirely rather than solving it
with an elliptical or physically-scaled 2-D kernel.
`_row_average_weights` is a Gaussian in pixel-index distance, not physical
mm distance — deliberately: within one row those are the same function up
to a fixed scale factor (`pixel_x_mm` is constant along a row), so the
kernel itself needs no pitch at all. `pixel_x_mm` is used for real exactly
once, in the GUI's options dialog, to show the window's physical width —
not in the kernel math, where it would only ever cancel out.
`_row_average_waveforms` is a masked/renormalized convolution (two
`correlate1d` calls, numerator and denominator, divided) rather than a
single fixed-normalized convolution, because a masked neighbor must
contribute *zero weight*, not a zero-amplitude sample at full weight — the
latter would bias every average near a masked run or a row's own edge
toward zero. The same two-correlation trick handles row-edge truncation for
free: `mode="constant", cval=0.0` zero-pads both the numerator and the
denominator beyond a row's own ends, so the output renormalizes by whatever
weight sum actually landed inside the row, no separate edge case.
Background subtraction stays exactly where it already was (subtracted once
from the fully-assembled `waves` buffer) rather than being threaded into the
per-neighbor gather. This is exact, not an approximation: because
`_row_average_waveforms`'s denominator is always the *actual* sum of
included, valid weights (never a fixed total), `Σwᵢ·(rawᵢ−bg) / Σwᵢ`
distributes to `avg − bg·(Σwᵢ/Σwᵢ) = avg − bg` regardless of which or how
many neighbors were included — subtracting background from the averaged
waveform is identical to subtracting it from every neighbor first, for any
window, at any row edge, with any number of masked-out neighbors.
No cross-row halo is needed: `compute_rf_image`'s chunk loop already splits
on rows only, and `read_row` already reads one row's complete
`(n_frames, spf)` slice at a time — averaging happens entirely inside that
one row's own frame axis, so a chunk boundary (which falls between rows)
can never truncate a window. Only a row's own start/end can, and that's the
same edge case the masked convolution already handles.
The averaging step doubles the live per-row scratch memory (a full-width
`(n_frames, spf)` buffer on top of the existing compacted `waves` buffer),
so `compute_rf_image` halves its byte budget when `row_avg_n > 0` before
`_plan_fft_rows` runs — see "Memory budget and row chunking" above. On the
largest real scans `_plan_chunks` is already
clamped to its floor of one row regardless, so this costs no concurrency
where it matters most; it mainly protects moderate-sized scans from an
unexpected regression.
Persistence: `cached_rf_image` (the extracted fast-path check) requires
`sras.precomputed_row_avg_n == row_avg_n` exactly, so a raw request can
never be silently served a row-averaged cache or vice versa, and a request
at one window size can never be served a cache at another — see
`scan_format.md`'s Cache Tail / CACH tail version history sections for the
on-disk `row_avg_n` field this depends on.
## Serving a stored cache: provenance, not just presence
A stored `peak_freq_mhz` image is only interchangeable with a live compute
for the *exact* settings it was computed under. Three of them are baked
irreversibly into the numbers — background subtraction, row-averaging window,
and zero-padding — so all three are recorded in the `SFFT` block and checked
by `cached_rf_image` before it hands the image back. Getting this wrong is
not a slow display, it is a *wrong* display, which is why the check is a
single predicate in one place rather than spread across callers.
The min peak frequency floor is the fourth recorded setting, and it sits
between the DC threshold (fully re-applicable at display time) and the
baked-in three: it is *tighten-only* re-applicable. A floor at or above the
stored one is served by masking stored pixels below it to the 0.0 "no valid
peak" sentinel — invalidated, not re-resolved; only a real recompute can
find their true above-floor peak. A floor *below* the stored one is a
genuine mismatch: the stored search never looked at those bins. Because a
served-then-masked image is a lossy stand-in for a real floored compute,
`cache_file` passes `use_stored=False` so a batch recompute always runs the
real FFT — otherwise re-batching a cached file would silently bake the
masked copy of its own old image back into the store.
Padding is the subtlest of the three, because a padded FFT looks like it
should be a refinement of the unpadded one. It isn't: zero-padding
interpolates between the natural bins, so it resolves a different peak
frequency for the same waveform. `precomputed_pad_factor` exists so a padded
view can be served from a cache computed at *its* pad while still refusing
one computed at any other, including pad 1. Before it existed the store was
pad 1 by definition and any `n_fft` was rejected outright — correct, but it
meant a user working at a pad factor got nothing at all from batch-computing
a file, which is most of the point of the feature. `pad_factor_for` maps an
`n_fft` request onto the integer factor a store could have recorded, and
returns 0 for a request that is not a whole multiple of `samples_per_frame`
— unmatchable by construction, since only integer factors are representable.
The batch actions therefore have to cache at the *viewer's* current pad
factor, not a fixed one: a cache stored at a pad nobody is viewing at is
dead weight. When the two do diverge (the user changes the pad after
batching), `_cache_mismatch_notes` says so in the scan info panel, because
the symptom otherwise is just "the file I pre-computed got slow again" with
no visible cause.
That divergence note is informational only, not a warning of an impending
recompute. The *display* path (`_stored_fft_image`) never asks
`cached_rf_image` whether a stored image matches the window's live
bg-sub/pad/row-averaging controls — it asks whether the image matches its
*own* recorded settings (`sras.precomputed_bg_sub`/`precomputed_pad_factor`/
`precomputed_row_avg_n`), which is always true whenever a stored image
exists. Two settings are taken live instead: the DC threshold and the min
peak frequency floor, both cheaply re-appliable as masks on serve. So for
bg-sub/pad/row-avg, presence alone decides whether a stored image is shown;
the floor is the one case where a live control can gate it — a live floor
*below* the stored one can't be answered by masking, so `cached_rf_image`
declines and the caller falls through to a real compute. They still matter for two things: a genuinely never-computed
angle's first live compute, and an explicit batch recompute — both of which
read the live controls and produce new stored data, at which point it's the
new data's *own* settings that get self-matched from then on. This is what
keeps a view switch (angle, channel, or flipping bg-sub/pad) from ever
discarding precomputed data — only an explicit batch recompute does, and it
already reloads the file afterward so the new data displays immediately.
Row-averaging has no live control to diverge from in the first place (it's
only ever set inside the batch dialog), so it never appears in the
divergence note — the "Cached images" line's own `row-averaged n=…` phrase
already covers it.
### Two caches, in cost order
`_refresh_display` consults this window's in-session `_fft_cache`/`_dc_cache`
dicts first, then the open file's stored v5/v7 blocks, and only then
dispatches a `ComputeWorker`. The second tier is what makes a batch-computed
file worth having; without it every angle change queued a worker and a
progress popup for an image already sitting on disk — precisely the cost the
batch was run to avoid. (The stored image *was* reachable before, but only
from inside `ComputeWorker`, i.e. after paying for the thread and the popup.)
The file tier asks with `allow_dc_recompute=False`. If applying the DC4 mask
would mean reading a whole CH4 channel, it declines rather than blocking the
GUI thread, and the fall-through worker reaches the same stored image via
`compute_rf_image` and pays for the mask off-thread. So the GUI thread never
does I/O, and the slow path is still a fast path.
## Angle alignment coordinate frames (`sras_compute.py`)
Alignment puts every angle's images onto one shared, zero-padded pixel grid
using a rigid transform only — rotation + translation, never scale.
Angle 0 (the reference) is the sole coordinate authority: it is the only
angle whose stage XY (`x_start_mm` / `y_positions_mm`) is ever read, and the
shared canvas is literally an extension of angle 0's own pixel grid, so the
aligned view carries angle 0's real X/Y axes. Every *other* angle is placed
purely by content — its rotation and translation come from cross-correlating
its CH4 image against angle 0's (`register_angle_to_reference`) — and its own
stage XY is deliberately never consulted. That is not an oversight: the
rotation stage moves the sample relative to the scan window, so where a
window sat in stage coordinates says nothing about where the sample is, and
an earlier design that pivoted each angle on a signal-weighted centroid of
its own window put every angle on a ~20 mm circle around the optical center
instead of stacking them into one shape.
Only two coordinate frames exist:
* **local mm** — one angle's own physical frame: origin at the *center of its
own pixel array*, x along +column, y along +row, scaled by that angle's own
pitches. Carries no stage position whatsoever.
* **ref mm** — the reference angle's local mm. A registration result
`(rotation_deg, shift_mm)` is exactly the rigid map from an angle's local
mm to ref mm: `q = R(rotation_deg) @ l + shift_mm`. Stage coordinates
re-enter once, at the very end, when the canvas origin is converted to
angle 0's stage mm (`AlignmentResult.canvas_origin_mm`).
Rotation is done in mm, never on raw pixel indices: the x pitch
(`SrasFile.pixel_x_mm`, 5 µm on a real scan) and the y/row pitch (50 µm)
differ by 10×, so rotating the raw index grid would shear the image — an
unwanted anisotropic scale. Registration runs on a resampled *isotropic* grid
for the same reason, and every affine maps shared-grid index → mm → undo
rotation/shift → that angle's own local mm → that angle's own raw index,
matching the output→input convention `scipy.ndimage.affine_transform` wants.
### Cropping the canvas is index translation, not a second transform
`crop_alignment_result` restricts an `AlignmentResult` to a rectangular window
of its canvas by folding the crop into each angle's existing affine rather than
composing a new one. From `_affine_out_to_src`, `matrix = D @ Rinv @ A_out`
depends only on the pitches and the rotation, and `A_out @ [row0, col0]` is
exactly the mm displacement of the new origin, so
```
matrix @ [r', c'] + (offset + matrix @ [row0, col0])
== matrix @ [r' + row0, c' + col0] + offset
```
identically. `matrix` is untouched and `offset` — which already absorbs the
origin — absorbs the crop too.
Two things follow, and both are relied on. `apply_alignment`, `reproject_mask`
and the aligned exporter all work on a cropped result with no special-casing:
resampling a cropped result is *exactly* a slice of resampling the full one
(`tests/test_align_export.py::test_crop_is_a_window_of_the_full_canvas` asserts
bit equality). And because the crop offset is a whole number of canvas pixels,
`canvas_for_params`' snap invariant — the reference angle lands on integer
canvas pixels — survives the crop, which is what keeps the reference exportable
as a verbatim block.
## Aligned export (`sras_align_export.py`)
`write_aligned_sras` bakes an alignment into a new v6 file: every angle
resampled onto the cropped shared canvas, so all of them end up with identical
geometry and the file opens already aligned. It is the only place in the
codebase that *resamples* waveform data — `sras_edit_scans` and `sras_average`
copy waveform bytes verbatim — which is why it is its own top-level module
rather than part of `sras_format` (scoped to the versioned binary spec, per the
sidecar section's own rule) or `sras_compute` (imported by every
multiprocessing child).
**Nearest neighbour, never interpolation.** Each output pixel gets exactly one
source pixel's three waveforms, verbatim. Averaging two neighbouring CH1
packets would synthesise a waveform the instrument never measured, whose FFT
peak is the peak of neither — meaningless for a technique whose entire output is
that peak frequency. The cost is that some source pixels are duplicated and
others dropped, which is the same trade `apply_alignment`'s `order=0` already
makes for the display.
**The rounding rule is `floor(x + 0.5)`, not `np.rint`.** `scipy.ndimage`'s
`order=0` rounds halves away from zero while `np.rint` rounds them to even. The
canvas is snapped to the reference's own pixel grid, so an angle whose row pitch
differs from the reference's lands on exact half-integers across whole rows —
this is the common case, not a corner case. Getting it wrong shifts those rows
by one source pixel relative to what the Aligned View drew.
**Out-of-bounds is tested on the fractional coordinate, not the rounded index.**
`scipy`'s `mode="constant"` writes `cval` wherever the coordinate leaves the
range of sample *centres*, `[0, n-1]` — a coordinate of −0.4 rounds to a
perfectly valid index 0 and is still padding. Testing the rounded index instead
puts a one-pixel rim of real data everywhere the preview shows padding.
**...but with a tolerance (`_EDGE_TOL`).** The affine is built from a chain of
mm-space multiplications, so an exactly-integer transform comes out a few times
1e-13 off: the reference angle's offset is `-20 - 7e-15`, not `-20`. A bare
`>= 0.0` therefore rejects that angle's entire first row, and `<= n-1` its last
column — for the *reference* angle, whose whole job is to pass through as an
exact integer crop. The tolerance is ~7 orders of magnitude above that noise and
~7 below the half-pixel scale at which a rounding decision means anything, so it
can only ever change pixels whose scipy answer was itself decided by noise.
**Padding is the per-channel ADC code nearest 0 mV, not 0.** Zero ADC decodes to
`(0 - yoff) * ymult + yzero`, which on real calibration is around +100 mV —
above any sensible CH4 mask threshold, so a zero fill would paint a solid
rectangle of "valid" pixels around the sample and corrupt every DC image and ROI
statistic downstream.
**Source rows are served from sliding in-RAM bands** (`_SourceReader`). A
rotated angle maps one output row to a *diagonal* across the source array, so
the pixels of a single output row come from hundreds of different source rows —
~1.4 MB each on a full-size scan. Indexing a memmap pixel-by-pixel in output
order re-faults nearly the whole angle per output row: terabytes of paging for a
gigabyte of data. Reading a contiguous band per output chunk, with the band
advancing monotonically, costs roughly 2× the source size in total reads.
**Writes go to `.part` and are `os.replace`d into position.** Not politeness: a
truncated .sras is not detectably broken, because `_parse_v6` drops incomplete
trailing angle blocks and opens what is left as an aborted scan. A half-written
export left in place would silently look like a real file with fewer angles.
**The Angle Table is carried over unchanged.** Alignment removes the *spatial*
rotation of the sample; it does not change which acoustic propagation direction
each angle measured, and that direction is the scientific content of a
multi-angle scan. Zeroing the table would make the export self-consistent for
re-registration and useless for anisotropy work. The consequence is that
re-registering an export needs `seed_deg=0.0` to put 0° inside the coarse sweep,
since `nominal_delta_deg` is still non-zero — which is exactly what the seed
parameter exists for.
## Alignment wizard (`sras_viewer/align_wizard.py`)
A `QWizard` rather than another dialog because the three steps are genuinely
sequential and the last one is destructive: correlate, choose a crop, write a
file. It replaces both former Fusion actions, so it also absorbs the old
`ManualAlignmentDialog`'s by-eye nudge editor — otherwise a scan the search
cannot fit would have no fallback at all.
Shared state lives on the wizard object, not in `registerField`: the pages pass
numpy arrays, `ManualAngleParams` and an `AlignmentResult` between them, none of
which are scalar widget properties.
`IndependentPages` is deliberately left **off**. With it set Qt never calls
`cleanupPage`, and `cleanupPage` is how the ROI page discards a crop when the
user goes back to re-correlate — a crop is indexed in canvas pixels, and a new
rotation means a different canvas, so stale indices would silently be
reinterpreted against the wrong grid. `geometry_generation` is the belt-and-
braces check for the same hazard.
The mask-stack preview shares the **final** canvas's origin and uses a pitch
that is an integer multiple of it, unlike the old manual dialog's padded,
unsnapped preview canvas. That is what lets the crop page convert a rectangle
drawn in millimetres into an exact integer window of the real canvas, with no
second coordinate frame to reconcile.
"Fit to full overlap" uses `largest_rect_at_least`, a largest-rectangle sweep,
not a bounding box of the fully-covered pixels. The full-overlap region of
several rotated scans is roughly a disc, and its bounding box has corners no
angle covers — offering that as the crop would hand the user the padding they
were trying to avoid.
Every background launch follows the two rules `_run_worker`'s docstring
establishes: disable the trigger *before* the call (so a re-entrant click cannot
start a second thread over the first), and never ignore the returned bool.
Progress is an inline `QProgressBar` on the page rather than a `QProgressDialog`
— a window-modal popup over a wizard both looks wrong and reintroduces the
event-loop pumping hazard that ordering exists to avoid. `reject()` refuses to
close while a job is in flight, since the running worker's signals are connected
to bound methods of the pages Qt would be deleting.
## Manual-alignment sidecar (`sras_compute.py`)
`<name>.sras.align.json` lives next to the scan file. The code lives in
`sras_compute`, not `sras_format`: `sras_format` is scoped to the versioned
binary .sras spec itself (see `scan_format.md`), while a manual alignment is
a viewer-computed *derived* artifact, analogous in kind to `AlignmentResult`
— so it belongs with the alignment math it serialises. json + pathlib are
stdlib, so this adds no dependency to a module whose load-bearing constraint
is staying free of Qt/matplotlib for cheap multiprocessing-child imports.
### Schema history
The stored `rotation_deg`/`shift_mm` are meaningless without the frame they
were measured in, so `_SIDECAR_SCHEMA_VERSION` is bumped whenever that frame
changes. Each bump makes older files describe a different (and, for the bugs
each bump fixed, actively wrong) transform than the same numbers would today;
loading one unchanged would silently reproduce the very "scans show up
everywhere" symptom the bump fixed — so older sidecars are treated as absent
rather than migrated.
* **1 → 2** — pivot moved from the scan-window bbox center to a
content-derived centroid, and the rotation sign convention was corrected.
* **2 → 3** — the content centroid was abandoned entirely: rotation is now
about each angle's own array center, mapped onto the reference's array
center, with `shift_mm` in the reference's local mm frame. No angle but the
reference contributes stage coordinates any more.