Add scan-editing CLI and alignment tests; extend Manual Alignment correlation
Continues the Manual Alignment work: refines the FFT cross-correlation and mask handling, adds sras_edit_scans.py (drop/renumber bad angle scans), tools/test_alignment.py (registration ground-truth suite), and a rotating test fixture in make_test_sras.py. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This commit is contained in:
Thomas Ales
@@ -125,6 +125,138 @@ def write(path: Path, n_angles: int = 3, seed: int = 0,
|
||||
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 _rot(theta_deg: float) -> np.ndarray:
|
||||
t = np.radians(theta_deg)
|
||||
c, s = np.cos(t), np.sin(t)
|
||||
return np.array([[c, -s], [s, c]])
|
||||
|
||||
|
||||
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}
|
||||
|
||||
|
||||
HDR_FMT_LEGACY = ">4sBHHffffIIdBB"
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user