sras-viewer/tools/make_test_sras.py

#!/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()