Thomas Ales
00a7afade0
- SrasFile: _parse_v6 now retains per-angle x_delta and the verbatim preamble/background byte spans, and gains public data_offset, y_pos_per_angle, and iter_angle_blocks() (which now owns the ragged block-offset walk used three separate places before). - sras_edit_scans: the 70-line re-parse of the v6 header sections (_read_v6_sections/_reread_span/_v6_angle_offsets) collapses into a _write_v6 that consumes SrasFile directly — verified byte-identical round-trip on v6 int8/int16 and legacy files. Eight print-and-exit pairs become _die(). - tools/make_test_sras imports the struct layouts from sras_format and the rotation matrix from sras_compute instead of re-declaring them (byte assembly stays independent of the reader). - sras_compute: block_mean_2d, pixel_pitch_mm, nominal_delta_deg made public (they were GUI-facing); registration_workers() and default_max_workers() wrap the remaining private reach-throughs from sras_workers. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
327 lines
14 KiB
Python
327 lines
14 KiB
Python
#!/usr/bin/env python3
|
|
"""Generate small synthetic .sras files for testing.
|
|
|
|
Writes v6 files (per-angle geometry, ragged waveform blocks) matching
|
|
scan_format.md, with deterministic pseudo-random waveform content so a test
|
|
can compute expected DC/FFT images independently of the reader under test.
|
|
|
|
Usage:
|
|
python tools/make_test_sras.py out.sras [--angles 3] [--seed 0]
|
|
"""
|
|
|
|
import argparse
|
|
import struct
|
|
import sys
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
|
|
sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
|
|
|
|
# Single source of truth for the byte layout: the reader's own constants.
|
|
# The byte *assembly* below stays independent, so a writer bug can't be
|
|
# masked by a matching reader bug.
|
|
from sras_format import HDR_FMT as HDR_FMT_LEGACY # noqa: E402
|
|
from sras_format import GEO_FMT_V6, HDR_FMT_V6 # noqa: E402
|
|
from sras_compute import _rotation_matrix as _rot # noqa: E402
|
|
|
|
# Per-angle (n_rows, n_frames) — deliberately different per angle so ragged
|
|
# geometry handling is actually exercised.
|
|
_GEOMETRY = [(5, 7), (4, 11), (6, 9), (3, 13), (7, 6)]
|
|
|
|
_SAMPLE_RATE_HZ = 6.25e9
|
|
_VELOCITY_MM_S = 20.0
|
|
_LASER_FREQ_HZ = 1000.0
|
|
_ROW_SPACING_MM = 0.05
|
|
|
|
|
|
def _preamble(ymult_v: float, yoff_adc: float, yzero_v: float) -> bytes:
|
|
"""A Tektronix WFMOutpre string in verbose (keyword) form — the reader
|
|
pulls YMULT/YOFF/YZERO out of it by name, so the keywords must be
|
|
present literally. YMULT/YZERO are in volts, as the scope reports them."""
|
|
return (
|
|
":WFMOUTPRE:BYT_NR 1;BIT_NR 8;ENCDG BIN;BN_FMT RI;BYT_OR MSB;"
|
|
'WFID "Ch1, DC coupling";NR_PT 2500;PT_FMT Y;'
|
|
"XINCR 1.6000E-10;XZERO 0.0E0;XUNIT \"s\";"
|
|
f"YMULT {ymult_v:.6E};YOFF {yoff_adc:.6E};YZERO {yzero_v:.6E};"
|
|
'YUNIT "V"'
|
|
).encode("utf-8")
|
|
|
|
|
|
def build(n_angles: int, seed: int, samples_per_frame: int,
|
|
geometry: list[tuple[int, int]] | None = None,
|
|
bps: int = 1) -> tuple[bytes, dict]:
|
|
rng = np.random.default_rng(seed)
|
|
src_geom = geometry or _GEOMETRY
|
|
geom = [src_geom[a % len(src_geom)] for a in range(n_angles)]
|
|
n_ch = 3
|
|
|
|
angles_deg = np.linspace(0.0, 60.0, n_angles, dtype=np.float32)
|
|
# Distinct calibration per channel so a swapped-channel bug is visible.
|
|
cal = [
|
|
(1.5625e-3, -87.04, 0.0),
|
|
(2.0000e-3, -60.00, 1.0e-3),
|
|
(2.5000e-3, -40.00, -2.0e-3),
|
|
]
|
|
|
|
out = bytearray()
|
|
out += struct.pack(
|
|
HDR_FMT_V6, b"SRAS", 6, n_angles,
|
|
0.0, 0.0, 1.0, 1.0, _ROW_SPACING_MM,
|
|
_VELOCITY_MM_S, _LASER_FREQ_HZ,
|
|
samples_per_frame, _SAMPLE_RATE_HZ, bps, n_ch,
|
|
)
|
|
out += angles_deg.astype(">f4").tobytes()
|
|
|
|
x_starts = []
|
|
for a, (n_rows, n_frames) in enumerate(geom):
|
|
x_start = -0.5 + 0.1 * a
|
|
x_starts.append(x_start)
|
|
out += struct.pack(GEO_FMT_V6, x_start, 1.0, n_frames, n_rows)
|
|
|
|
y_positions = []
|
|
for a, (n_rows, _) in enumerate(geom):
|
|
y = (0.2 * a + np.arange(n_rows) * _ROW_SPACING_MM).astype(np.float32)
|
|
y_positions.append(y)
|
|
out += y.astype(">f4").tobytes()
|
|
|
|
for ymult_v, yoff, yzero_v in cal:
|
|
p = _preamble(ymult_v, yoff, yzero_v)
|
|
out += struct.pack(">H", len(p)) + p
|
|
|
|
background = rng.integers(-8, 9, size=samples_per_frame, dtype=np.int8)
|
|
out += struct.pack(">I", samples_per_frame) + background.tobytes()
|
|
|
|
# Waveform data. CH1 gets a sinusoid at a per-pixel frequency so the FFT
|
|
# peak is predictable; CH3/CH4 get per-pixel DC levels so the mean is too.
|
|
t = np.arange(samples_per_frame)
|
|
waveforms = []
|
|
for a, (n_rows, n_frames) in enumerate(geom):
|
|
block = np.empty((n_rows, n_ch, n_frames, samples_per_frame), dtype=np.int8)
|
|
for r in range(n_rows):
|
|
for f in range(n_frames):
|
|
bin_idx = 3 + ((a + r + f) % 17)
|
|
phase = 2 * np.pi * bin_idx * t / samples_per_frame
|
|
block[r, 0, f] = np.clip(
|
|
np.round(60 * np.sin(phase)), -128, 127).astype(np.int8)
|
|
block[r, 1, f] = np.int8((a * 7 + r * 3 + f) % 100 - 50)
|
|
block[r, 2, f] = np.int8((a * 5 + r * 11 + f * 2) % 120 - 60)
|
|
waveforms.append(block)
|
|
# bps=2 stores the same values big-endian int16, exercising the
|
|
# reader's >i2 memmap path.
|
|
out += (block.astype(">i2") if bps == 2 else block).tobytes()
|
|
|
|
meta = {
|
|
"n_angles": n_angles,
|
|
"geometry": geom,
|
|
"angles_deg": angles_deg,
|
|
"x_starts": x_starts,
|
|
"y_positions": y_positions,
|
|
"cal": cal,
|
|
"background": background,
|
|
"waveforms": waveforms,
|
|
"samples_per_frame": samples_per_frame,
|
|
"sample_rate_hz": _SAMPLE_RATE_HZ,
|
|
}
|
|
return bytes(out), meta
|
|
|
|
|
|
def write(path: Path, n_angles: int = 3, seed: int = 0,
|
|
samples_per_frame: int = 64,
|
|
geometry: list[tuple[int, int]] | None = None,
|
|
bps: int = 1) -> dict:
|
|
payload, meta = build(n_angles, seed, samples_per_frame, geometry, bps=bps)
|
|
path.write_bytes(payload)
|
|
return meta
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Rotating-sample scan: one shape, imaged at several known rotations
|
|
# ---------------------------------------------------------------------------
|
|
#
|
|
# The scan the angle-alignment path actually has to solve: every angle images
|
|
# the *same* sample at a different known rotation and offset, and a correct
|
|
# alignment stacks them all back into one shape. Two properties are
|
|
# deliberately hostile:
|
|
#
|
|
# * every angle gets a different window size and a different, meaningless
|
|
# stage x_start / y0 — alignment must ignore per-angle stage coordinates
|
|
# entirely, so any code that reads them will visibly fail here;
|
|
# * the pixel grid is strongly anisotropic (5 µm along x, 50 µm along y),
|
|
# like the real instrument, so any registration that rotates raw indices
|
|
# instead of millimetres shears the image and cannot converge.
|
|
|
|
_ROT_DX_MM = 0.005 # x pitch, from velocity/laser_freq below
|
|
_ROT_DY_MM = 0.05 # row spacing
|
|
_ROT_BG_MV = 4.0
|
|
_ROT_FG_MV = 160.0
|
|
|
|
|
|
# How far the sample sits from the rotation axis. Non-zero on purpose: on the
|
|
# real instrument every angle's scan window is centred on the rotation axis
|
|
# while the sample is not, so each scan sees the sample somewhere else along a
|
|
# circle. That offset is exactly what a wrong rotation pivot turns into a ring
|
|
# of scans instead of a stack, so a centred test sample would hide the bug.
|
|
_ROT_SAMPLE_OFFSET_MM = (0.55, 0.40)
|
|
|
|
|
|
def _sample_shape_mv(u: np.ndarray, v: np.ndarray) -> np.ndarray:
|
|
"""An asymmetric test sample in its own mm frame, chirally distinct at
|
|
every rotation (no 180° ambiguity) and with structure at several radii so
|
|
rotation is well determined."""
|
|
u = u - _ROT_SAMPLE_OFFSET_MM[0]
|
|
v = v - _ROT_SAMPLE_OFFSET_MM[1]
|
|
img = np.full(u.shape, _ROT_BG_MV, dtype=np.float32)
|
|
img[((u / 0.85) ** 2 + (v / 0.40) ** 2) <= 1.0] = _ROT_FG_MV # bar
|
|
img[(np.abs(u - 0.55) <= 0.22) & (np.abs(v - 0.62) <= 0.22)] = _ROT_FG_MV # nub
|
|
img[((u + 0.75) ** 2 + (v + 0.30) ** 2) <= 0.20 ** 2] = _ROT_FG_MV # dot
|
|
return img
|
|
|
|
|
|
def write_rotating(path: Path, n_angles: int = 5, samples_per_frame: int = 4,
|
|
seed: int = 0) -> dict:
|
|
"""Write a v6 file whose CH4 DC image is one sample seen at n_angles known
|
|
rotations, and return the ground truth each angle should register to.
|
|
|
|
``truth[a] = (rotation_deg, (shift_x_mm, shift_y_mm))`` is the rigid map
|
|
from angle *a*'s local mm (origin at its own array center) to angle 0's —
|
|
exactly what ``register_angle_to_reference`` is supposed to recover.
|
|
"""
|
|
rng = np.random.default_rng(seed)
|
|
n_ch, bps = 3, 1
|
|
cal = [(1.5625e-3, -87.04, 0.0), (2.0e-3, -60.0, 1.0e-3), (2.5e-3, -40.0, -2.0e-3)]
|
|
ymult_mv, yoff, yzero_mv = cal[2][0] * 1000, cal[2][1], cal[2][2] * 1000
|
|
|
|
stage_angles, geom, x_starts, y_starts, thetas, offsets = [], [], [], [], [], []
|
|
for a in range(n_angles):
|
|
stage = -37.0 * a # what the rotation stage reports
|
|
stage_angles.append(stage)
|
|
# The true image rotation is the negative of the stage's reported
|
|
# angle: the stage's positive sense is the opposite of math-positive
|
|
# (x toward y) in scan mm. Nothing may depend on knowing that — the
|
|
# registration search tries both signs.
|
|
thetas.append(-stage)
|
|
offsets.append((0.0, 0.0) if a == 0
|
|
else (float(rng.uniform(-0.3, 0.3)), float(rng.uniform(-0.3, 0.3))))
|
|
# A different window per angle, all centred on the same array center —
|
|
# the real instrument grows each angle's axis-aligned bounding box to
|
|
# cover the rotated ROI. Sized so the off-axis sample stays inside every
|
|
# window at every angle, keeping the expected result unambiguous.
|
|
geom.append((88 + 8 * a, 780 + 60 * a))
|
|
# Meaningless per-angle stage positions: correct alignment never reads
|
|
# them, so scattering them proves it.
|
|
x_starts.append(float(20.0 + rng.uniform(-6.0, 6.0)))
|
|
y_starts.append(float(30.0 + rng.uniform(-6.0, 6.0)))
|
|
|
|
out = bytearray()
|
|
out += struct.pack(
|
|
HDR_FMT_V6, b"SRAS", 6, n_angles,
|
|
x_starts[0], y_starts[0], 1.0, 1.0, _ROT_DY_MM,
|
|
_VELOCITY_MM_S, _VELOCITY_MM_S / _ROT_DX_MM, # velocity/freq -> 5 µm pitch
|
|
samples_per_frame, _SAMPLE_RATE_HZ, bps, n_ch,
|
|
)
|
|
out += np.array(stage_angles, dtype=">f4").tobytes()
|
|
for a, (n_rows, n_frames) in enumerate(geom):
|
|
out += struct.pack(GEO_FMT_V6, x_starts[a], 1.0, n_frames, n_rows)
|
|
for a, (n_rows, _) in enumerate(geom):
|
|
out += (y_starts[a] + np.arange(n_rows) * _ROT_DY_MM).astype(">f4").tobytes()
|
|
for ymult_v, yoff_a, yzero_v in cal:
|
|
p = _preamble(ymult_v, yoff_a, yzero_v)
|
|
out += struct.pack(">H", len(p)) + p
|
|
background = rng.integers(-8, 9, size=samples_per_frame, dtype=np.int8)
|
|
out += struct.pack(">I", samples_per_frame) + background.tobytes()
|
|
|
|
truth, dc4_images = {}, []
|
|
for a, (n_rows, n_frames) in enumerate(geom):
|
|
# Local mm of every pixel, measured from this angle's own array center.
|
|
lx = (np.arange(n_frames) - (n_frames - 1) / 2.0) * _ROT_DX_MM
|
|
ly = (np.arange(n_rows) - (n_rows - 1) / 2.0) * _ROT_DY_MM
|
|
gx, gy = np.meshgrid(lx, ly)
|
|
# local = R(theta) @ sample + offset, so sample = R(theta)^T @ (local - offset)
|
|
rel = np.stack([gx - offsets[a][0], gy - offsets[a][1]], axis=-1)
|
|
s = rel @ _rot(thetas[a]) # == rel @ R^T.T == R^T @ rel
|
|
dc4 = _sample_shape_mv(s[..., 0], s[..., 1])
|
|
dc4_images.append(dc4)
|
|
|
|
inv = _rot(-thetas[a])
|
|
truth[a] = (-thetas[a],
|
|
tuple(float(v) for v in -(inv @ np.array(offsets[a]))))
|
|
|
|
adc4 = np.clip(np.round((dc4 - yzero_mv) / ymult_mv + yoff), -128, 127).astype(np.int8)
|
|
block = np.zeros((n_rows, n_ch, n_frames, samples_per_frame), dtype=np.int8)
|
|
block[:, 2] = adc4[:, :, None] # CH4 carries the sample
|
|
block[:, 1] = 10 # CH3 flat
|
|
block[:, 0] = rng.integers(-40, 41, size=(n_rows, n_frames, samples_per_frame),
|
|
dtype=np.int8) # CH1 noise
|
|
out += block.tobytes()
|
|
|
|
path.write_bytes(bytes(out))
|
|
return {"n_angles": n_angles, "geometry": geom, "stage_angles_deg": stage_angles,
|
|
"truth": truth, "dc4_mv": dc4_images, "x_starts": x_starts,
|
|
"y_starts": y_starts, "dx_mm": _ROT_DX_MM, "dy_mm": _ROT_DY_MM}
|
|
|
|
|
|
|
|
|
|
def write_legacy(path: Path, version: int = 4, n_angles: int = 2,
|
|
n_rows: int = 4, n_frames: int = 10,
|
|
samples_per_frame: int = 32, seed: int = 0) -> dict:
|
|
"""Write a v2/v3/v4 file: uniform geometry, one flat waveform block.
|
|
|
|
Used to exercise sras_average.py, which only handles the legacy formats.
|
|
"""
|
|
rng = np.random.default_rng(seed)
|
|
n_ch, bps = 3, 1
|
|
|
|
out = bytearray()
|
|
out += struct.pack(
|
|
HDR_FMT_LEGACY, b"SRAS", version, n_angles, n_rows,
|
|
-0.5, 1.0, _VELOCITY_MM_S, _LASER_FREQ_HZ,
|
|
n_frames, samples_per_frame, _SAMPLE_RATE_HZ, bps, n_ch,
|
|
)
|
|
angles = np.linspace(0.0, 45.0, n_angles, dtype=np.float32)
|
|
out += angles.astype(">f4").tobytes()
|
|
y = (np.arange(n_rows) * _ROW_SPACING_MM).astype(np.float32)
|
|
out += y.astype(">f4").tobytes()
|
|
|
|
if version >= 3:
|
|
for ymult_v, yoff, yzero_v in ((1.5625e-3, -87.04, 0.0),
|
|
(2.0e-3, -60.0, 1.0e-3),
|
|
(2.5e-3, -40.0, -2.0e-3))[:n_ch]:
|
|
p = _preamble(ymult_v, yoff, yzero_v)
|
|
out += struct.pack(">H", len(p)) + p
|
|
|
|
background = rng.integers(-8, 9, size=samples_per_frame, dtype=np.int8)
|
|
if version >= 4:
|
|
out += struct.pack(">I", samples_per_frame) + background.tobytes()
|
|
|
|
data = rng.integers(-100, 101,
|
|
size=(n_angles, n_rows, n_ch, n_frames, samples_per_frame),
|
|
dtype=np.int8)
|
|
out += data.tobytes()
|
|
path.write_bytes(bytes(out))
|
|
return {"version": version, "n_angles": n_angles, "n_rows": n_rows,
|
|
"n_frames": n_frames, "samples_per_frame": samples_per_frame,
|
|
"n_channels": n_ch, "data": data, "angles_deg": angles,
|
|
"y_positions": y, "background": background}
|
|
|
|
|
|
def main():
|
|
p = argparse.ArgumentParser(description=__doc__)
|
|
p.add_argument("output")
|
|
p.add_argument("--angles", type=int, default=3)
|
|
p.add_argument("--seed", type=int, default=0)
|
|
p.add_argument("--spf", type=int, default=64, help="samples per frame")
|
|
args = p.parse_args()
|
|
|
|
out = Path(args.output)
|
|
meta = write(out, args.angles, args.seed, args.spf)
|
|
print(f"Wrote {out} ({out.stat().st_size:,} bytes)")
|
|
print(f" angles : {meta['n_angles']}")
|
|
print(f" geometry : {meta['geometry']}")
|
|
print(f" spf : {meta['samples_per_frame']}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|