Thomas K Ales [MSE]
263 lines
11 KiB
Python
Executable File
263 lines
11 KiB
Python
Executable File
"""Middle-row SAW quality check: acquire one row per angle, then read the
|
||
alignment off the frequencies it produces.
|
||
|
||
Two halves of one test mode, kept together because neither is much use
|
||
without the other:
|
||
|
||
*Acquisition* — ``middle_row_plan`` reduces a full ScanPlan to a single row
|
||
per angle, the row-wise middle of the ROI. ScanEngine runs the result
|
||
exactly like any other scan and writes it as a v11 .sras file
|
||
(``sras_format.VERSION_SAW_CHECK``), so a check costs one row-time per angle
|
||
instead of the hours a full multi-angle scan takes.
|
||
|
||
*Analysis* — ``frequency_traces`` turns such a file back into one peak-SAW-
|
||
frequency trace per angle, and ``alignment_summary`` reduces those to the
|
||
numbers the operator is actually asking about. Both are Qt-free; the plotting
|
||
lives in saw_check_viewer.py.
|
||
|
||
Why the middle row answers an alignment question: ``scan_geometry.build_plan``
|
||
centres every angle's rotated bounding box on the same nominal ROI centre, so
|
||
each angle's middle row crosses that one point on the sample. Every angle
|
||
therefore measures the same material, and a spread in the per-angle
|
||
frequencies is a property of the rig (or of a genuinely anisotropic sample),
|
||
not of where each row happened to land.
|
||
"""
|
||
from __future__ import annotations
|
||
|
||
from dataclasses import dataclass, field, replace
|
||
|
||
import numpy as np
|
||
|
||
from core.scan_geometry import ScanGeometryError, ScanPlan
|
||
from core.sras_analysis import ChannelCalibration, compute_rf_image
|
||
from core.sras_format import SrasFile
|
||
|
||
# Rules of thumb for the read-out, not physics. A well-aligned rig on an
|
||
# isotropic sample reads the same frequency at every angle, so the spread of
|
||
# the per-angle medians is the alignment signal — but an anisotropic sample
|
||
# genuinely varies with angle, so a wide spread is a prompt to look at the
|
||
# curves, never a verdict on its own.
|
||
SPREAD_GOOD_PCT = 1.0
|
||
SPREAD_MARGINAL_PCT = 3.0
|
||
# Below this fraction of unmasked pixels a trace is too sparse to read.
|
||
VALID_FRACTION_FLOOR = 0.5
|
||
|
||
|
||
# ── Acquisition side ─────────────────────────────────────────────────────────
|
||
|
||
def middle_row_plan(plan: ScanPlan) -> ScanPlan:
|
||
"""Reduce a scan plan to its row-wise middle row at every angle.
|
||
|
||
Each angle keeps the geometry the full scan would have used — same
|
||
x_start, x_delta and n_frames from its own rotated bounding box — and
|
||
scans only the middle entry of its row list, so the check samples exactly
|
||
what the scan would along that row.
|
||
|
||
An even row count has no exact middle; the upper of the two central rows
|
||
is taken (``n_rows // 2``), which is also the row the viewer picks when it
|
||
reads the middle row out of a full scan.
|
||
"""
|
||
if plan.n_angles == 0:
|
||
raise ScanGeometryError("Cannot build a SAW check from a plan with no angles")
|
||
|
||
per_angle = []
|
||
for pa in plan.per_angle:
|
||
if not pa.y_positions:
|
||
raise ScanGeometryError(
|
||
f"Angle {pa.angle_deg:.1f}° has no rows, so it has no middle row to check"
|
||
)
|
||
per_angle.append(replace(pa, n_rows=1,
|
||
y_positions=[pa.y_positions[middle_row_index(pa.n_rows)]]))
|
||
return replace(plan, per_angle=per_angle)
|
||
|
||
|
||
def middle_row_index(n_rows: int) -> int:
|
||
"""The row this check calls the middle one. One rule, two callers."""
|
||
return max(0, n_rows // 2)
|
||
|
||
|
||
# ── Analysis side ────────────────────────────────────────────────────────────
|
||
|
||
@dataclass
|
||
class AngleTrace:
|
||
"""One angle's peak SAW frequency along its middle row.
|
||
|
||
``freq_mhz`` is NaN wherever the pixel was masked out (CH4 DC below the
|
||
threshold), so the gaps stay gaps instead of reading as 0 MHz.
|
||
"""
|
||
angle_idx: int
|
||
angle_deg: float
|
||
row_idx: int
|
||
y_mm: float
|
||
x_mm: np.ndarray # absolute stage X of each frame
|
||
freq_mhz: np.ndarray # NaN where masked
|
||
_valid: np.ndarray = field(init=False, repr=False)
|
||
|
||
def __post_init__(self):
|
||
self._valid = np.isfinite(self.freq_mhz)
|
||
|
||
@property
|
||
def offset_mm(self) -> np.ndarray:
|
||
"""X relative to the centre of this row.
|
||
|
||
Every angle's row is centred on the same ROI centre, so plotting
|
||
against this puts all the angles' curves over the same piece of
|
||
sample — which is the whole point of the comparison.
|
||
"""
|
||
if len(self.x_mm) == 0:
|
||
return self.x_mm
|
||
return self.x_mm - 0.5 * (self.x_mm[0] + self.x_mm[-1])
|
||
|
||
@property
|
||
def n_valid(self) -> int:
|
||
return int(self._valid.sum())
|
||
|
||
@property
|
||
def valid_fraction(self) -> float:
|
||
return self.n_valid / len(self.freq_mhz) if len(self.freq_mhz) else 0.0
|
||
|
||
@property
|
||
def median_mhz(self) -> float:
|
||
return float(np.median(self.freq_mhz[self._valid])) if self.n_valid else float("nan")
|
||
|
||
@property
|
||
def std_mhz(self) -> float:
|
||
return float(np.std(self.freq_mhz[self._valid])) if self.n_valid > 1 else float("nan")
|
||
|
||
@property
|
||
def drift_mhz_per_mm(self) -> float:
|
||
"""Least-squares slope of frequency along the row.
|
||
|
||
A flat trace means the response did not change across the ROI; a
|
||
sloped one is the signature of a tilt or a defocus the angle spread
|
||
alone would not show.
|
||
"""
|
||
if self.n_valid < 2:
|
||
return float("nan")
|
||
x = self.offset_mm[self._valid]
|
||
if np.ptp(x) == 0:
|
||
return float("nan")
|
||
return float(np.polyfit(x, self.freq_mhz[self._valid], 1)[0])
|
||
|
||
|
||
@dataclass
|
||
class AlignmentSummary:
|
||
"""What the per-angle traces say about the alignment, in scalars."""
|
||
n_angles: int
|
||
median_mhz: float
|
||
spread_mhz: float # max − min of the per-angle medians
|
||
spread_pct: float # that spread as a % of the overall median
|
||
best_angle_deg: float # angle reading the highest median
|
||
worst_angle_deg: float # angle reading the lowest median
|
||
worst_drift_mhz_per_mm: float
|
||
worst_drift_angle_deg: float
|
||
min_valid_fraction: float
|
||
|
||
@property
|
||
def level(self) -> str:
|
||
""""good" / "marginal" / "poor" — see the module's threshold note."""
|
||
if self.n_angles == 0 or not np.isfinite(self.spread_pct):
|
||
return "poor"
|
||
if self.min_valid_fraction < VALID_FRACTION_FLOOR:
|
||
return "poor"
|
||
if self.spread_pct <= SPREAD_GOOD_PCT:
|
||
return "good"
|
||
if self.spread_pct <= SPREAD_MARGINAL_PCT:
|
||
return "marginal"
|
||
return "poor"
|
||
|
||
def describe(self) -> str:
|
||
if self.n_angles == 0:
|
||
return "No angle produced a usable frequency trace."
|
||
if self.min_valid_fraction < VALID_FRACTION_FLOOR:
|
||
return (f"Only {self.min_valid_fraction * 100:.0f} % of the worst angle's row "
|
||
f"is above the DC threshold — check the detection beam and the "
|
||
f"threshold before reading the spread.")
|
||
return (f"Per-angle medians span {self.spread_mhz:.3f} MHz "
|
||
f"({self.spread_pct:.2f} % of {self.median_mhz:.3f} MHz), "
|
||
f"lowest at {self.worst_angle_deg:.1f}°, highest at {self.best_angle_deg:.1f}°. "
|
||
f"Largest drift along a row: {self.worst_drift_mhz_per_mm:+.3f} MHz/mm "
|
||
f"at {self.worst_drift_angle_deg:.1f}°.")
|
||
|
||
|
||
def frequency_traces(sras: SrasFile, *, dc_threshold_mv: float = 0.0,
|
||
subtract_background: bool = False,
|
||
gate_start_ns: float | None = None,
|
||
gate_end_ns: float | None = None,
|
||
calib: ChannelCalibration | None = None,
|
||
on_progress=lambda done, total: None) -> list[AngleTrace]:
|
||
"""Peak SAW frequency along the middle row of every angle in ``sras``.
|
||
|
||
Works on a v11 check (one row per angle, so the middle row is the only
|
||
row) and on a full scan alike — the same middle row the check would
|
||
have acquired is pulled out of the scan, which is what lets a finished
|
||
scan be re-examined with the check's own read-out.
|
||
|
||
``subtract_background`` takes each angle's own background out of its
|
||
frames (v6/v10 files have only the one, which every angle then shares).
|
||
Doing it per angle is the point of the per-angle capture: comparing
|
||
angles is exactly what this read-out is for, so they must not be
|
||
referenced against one background taken at whichever angle came first.
|
||
|
||
Angles with nothing on disk (an aborted file) are skipped rather than
|
||
reported as flat zero.
|
||
"""
|
||
calib = calib if calib is not None else ChannelCalibration.from_preambles(sras.preambles)
|
||
freq_axis = sras.freq_axis_mhz(sras.header.samples_per_frame)
|
||
time_axis = sras.time_axis_ns()
|
||
statuses = sras.angle_status()
|
||
|
||
traces: list[AngleTrace] = []
|
||
for st in statuses:
|
||
on_progress(st.index, len(statuses))
|
||
if st.n_rows_available < 1:
|
||
continue
|
||
pa = sras.per_angle[st.index]
|
||
row = middle_row_index(st.n_rows_available)
|
||
view = sras.load_angle(st.index, n_rows=st.n_rows_available)[row:row + 1]
|
||
|
||
background = sras.background_array(st.index) if subtract_background else None
|
||
img = compute_rf_image(view, calib, freq_axis, dc_threshold_mv,
|
||
background=background,
|
||
gate_start_ns=gate_start_ns, gate_end_ns=gate_end_ns,
|
||
time_axis_ns=time_axis)
|
||
# compute_rf_image zeroes masked pixels and its FFT never peaks in the
|
||
# suppressed DC bin, so 0 MHz means "no reading" and nothing else.
|
||
freq = img[0].astype(np.float64)
|
||
freq[freq <= 0.0] = np.nan
|
||
|
||
traces.append(AngleTrace(
|
||
angle_idx=st.index, angle_deg=pa.angle_deg, row_idx=row,
|
||
y_mm=pa.y_positions[row] if row < len(pa.y_positions) else float("nan"),
|
||
x_mm=sras.x_axis_mm(st.index), freq_mhz=freq,
|
||
))
|
||
on_progress(len(statuses), len(statuses))
|
||
return traces
|
||
|
||
|
||
def alignment_summary(traces: list[AngleTrace]) -> AlignmentSummary:
|
||
"""Reduce per-angle traces to the alignment read-out."""
|
||
usable = [t for t in traces if t.n_valid > 0]
|
||
if not usable:
|
||
nan = float("nan")
|
||
return AlignmentSummary(0, nan, nan, nan, nan, nan, nan, nan, 0.0)
|
||
|
||
medians = np.array([t.median_mhz for t in usable])
|
||
overall = float(np.median(medians))
|
||
spread = float(medians.max() - medians.min())
|
||
drifts = [(abs(t.drift_mhz_per_mm), t) for t in usable
|
||
if np.isfinite(t.drift_mhz_per_mm)]
|
||
worst_drift = max(drifts, key=lambda d: d[0])[1] if drifts else None
|
||
|
||
return AlignmentSummary(
|
||
n_angles=len(usable),
|
||
median_mhz=overall,
|
||
spread_mhz=spread,
|
||
spread_pct=spread / overall * 100.0 if overall else float("nan"),
|
||
best_angle_deg=usable[int(np.argmax(medians))].angle_deg,
|
||
worst_angle_deg=usable[int(np.argmin(medians))].angle_deg,
|
||
worst_drift_mhz_per_mm=worst_drift.drift_mhz_per_mm if worst_drift else float("nan"),
|
||
worst_drift_angle_deg=worst_drift.angle_deg if worst_drift else float("nan"),
|
||
min_valid_fraction=min(t.valid_fraction for t in usable),
|
||
)
|