diff --git a/sras_compute.py b/sras_compute.py index 803d2c1..6b49a6f 100644 --- a/sras_compute.py +++ b/sras_compute.py @@ -455,27 +455,41 @@ def _theta_deg(sras: SrasFile, angle_idx: int, ref_idx: int) -> float: return -float(sras.angles_deg[angle_idx] - sras.angles_deg[ref_idx]) -def _mask_centroid_mm(sras: SrasFile, angle_idx: int, - mask: np.ndarray) -> tuple[float, float]: - """Centroid (mean x, mean y) of *mask*'s True pixels, in angle_idx's own - local mm frame — the alignment pivot used in place of the raw scan- - window bbox center (see compute_pivot_points_mm). Falls back to the bbox - center if the mask has no pixels above threshold, since a centroid of - nothing is undefined.""" - if not mask.any(): +def _signal_centroid_mm(sras: SrasFile, angle_idx: int, + dc4_mv: np.ndarray) -> tuple[float, float]: + """Intensity-weighted centroid (mean x, mean y, weighted by CH4 signal + after subtracting this angle's own minimum) in angle_idx's own local mm + frame — the alignment pivot used in place of the raw scan-window bbox + center (see compute_pivot_points_mm). + + Weighting by the continuous DC signal, rather than a binary >= + dc_threshold_mv mask, means the pivot never depends on how well one + shared threshold happens to suit this particular angle: real signal + levels vary scan to scan, so a threshold tuned for one angle can leave + another angle's binary mask empty — and a centroid of an empty mask has + nothing to fall back to *except* the raw bbox center, silently + reproducing the exact "aligned to the scan window, not the sample" + problem this pivot exists to avoid. Falls back to the bbox center only + in the fully-degenerate case of a perfectly flat signal (nothing to + weight by at all). + """ + weights = dc4_mv - dc4_mv.min() + total = float(weights.sum()) + if total <= 0.0: return _bbox_center_mm(sras, angle_idx) x = sras.x_axis_mm(angle_idx) y = sras.y_positions_mm(angle_idx) - rows, cols = np.nonzero(mask) - return float(x[cols].mean()), float(y[rows].mean()) + cx = float((x * weights.sum(axis=0)).sum() / total) + cy = float((y * weights.sum(axis=1)).sum() / total) + return cx, cy -def compute_pivot_points_mm(sras: SrasFile, dc_threshold_mv: float, - masks: dict[int, np.ndarray] | None = None +def compute_pivot_points_mm(sras: SrasFile, + dc4_mv: dict[int, np.ndarray] | None = None ) -> dict[int, tuple[float, float]]: """Per-angle alignment pivot, in each angle's own local mm frame: the - centroid of its own CH4-threshold mask (its own sample footprint), - rather than the raw scan-window bbox center. + CH4-signal-weighted centroid of its own footprint (see + _signal_centroid_mm), rather than the raw scan-window bbox center. Pivoting on each angle's own content — instead of on wherever its scanned window happened to sit in microscope/global XY space — is what @@ -486,21 +500,23 @@ def compute_pivot_points_mm(sras: SrasFile, dc_threshold_mv: float, leaves a residual orbital motion between angles that a content-centroid pivot does not. - *masks* lets a caller that has already computed CH4 masks at this - threshold (compute_angle_alignment's Step 1) reuse them instead of - recomputing the DC image; angles missing from it are computed fresh via - dc_image_mv, which prefers a stored v5/v7 cache over recomputing from - raw waveforms. + Deliberately independent of dc_threshold_mv (the RF-mask / overlay + threshold): that value is a display/masking choice and must never + silently change where alignment pivots. + + *dc4_mv* lets a caller that has already computed each angle's CH4 mV + image (compute_angle_alignment's Step 1, or ManualAlignmentDialog's own + cache) reuse it instead of recomputing; angles missing from it are + computed fresh via dc_image_mv, which prefers a stored v5/v7 cache over + recomputing from raw waveforms. """ - masks = masks or {} + dc4_mv = dc4_mv or {} pivots: dict[int, tuple[float, float]] = {} for a in range(sras.n_angles): - mask = masks.get(a) - if mask is None: - mask = dc_image_mv(sras, a, CH4_IDX) >= dc_threshold_mv - else: - mask = mask > 0 # compute_angle_alignment's masks are float32 0.0/1.0 - pivots[a] = _mask_centroid_mm(sras, a, mask) + img = dc4_mv.get(a) + if img is None: + img = dc_image_mv(sras, a, CH4_IDX) + pivots[a] = _signal_centroid_mm(sras, a, img) return pivots @@ -544,13 +560,14 @@ def _build_affine_canvas_to_raw(sras: SrasFile, angle_idx: int, ref_idx: int, where theta is _theta_deg(angle_idx, ref_idx) unless *theta_deg* overrides it (see _corners_in_ref_frame's note, used by the manual- alignment path), c_ref/c_a are each angle's own alignment pivot - (pivot_mm — see compute_pivot_points_mm; the centroid of its own - CH4-threshold mask, not the raw scan-window bbox center — this keeps - the sample itself centered post-rotation, minimizing required canvas - padding and, more importantly, keeping alignment relative to the sample - rather than to wherever the scan window sat in microscope/global XY - space), and A_out/D are the index<->mm scaling matrices for the canvas - pitch and this angle's own native pitch respectively. + (pivot_mm — see compute_pivot_points_mm; the CH4-signal-weighted + centroid of its own footprint, not the raw scan-window bbox center — + this keeps the sample itself centered post-rotation, minimizing + required canvas padding and, more importantly, keeping alignment + relative to the sample rather than to wherever the scan window sat in + microscope/global XY space), and A_out/D are the index<->mm scaling + matrices for the canvas pitch and this angle's own native pitch + respectively. """ theta = _theta_deg(sras, angle_idx, ref_idx) if theta_deg is None else theta_deg Rinv = _rotation_matrix(theta).T @@ -671,21 +688,25 @@ def compute_angle_alignment(sras: SrasFile, ref_angle_idx: int, _ticks.append(1) progress_cb(base + int(min(len(_ticks), n) / n * span)) - # ---- Step 1: binarized CH4 mask per angle, native per-angle grid ----- - def mask_for(a: int) -> np.ndarray: + # ---- Step 1: CH4 DC image + binarized mask per angle, native grid ----- + def dc4_for(a: int) -> np.ndarray: dc4 = adc_to_mv( compute_dc_image(sras, a, CH4_IDX, max_workers=1, budget=angle_budget), *sras.cal(CH4_IDX)) - m = (dc4 >= dc_threshold_mv).astype(np.float32) tick(0, 25) - return m + return dc4 - masks = dict(enumerate(_parallel_map(mask_for, range(n), n_workers))) + dc4_mv = dict(enumerate(_parallel_map(dc4_for, range(n), n_workers))) + masks = {a: (img >= dc_threshold_mv).astype(np.float32) + for a, img in dc4_mv.items()} - # Each angle's own alignment pivot — the centroid of its own mask, not - # the raw scan-window bbox center (see compute_pivot_points_mm) — reuses - # the masks just computed above, so this is free (no extra DC compute). - pivot_mm = compute_pivot_points_mm(sras, dc_threshold_mv, masks=masks) + # Each angle's own alignment pivot — the CH4-signal-weighted centroid, + # not the raw scan-window bbox center (see compute_pivot_points_mm). + # Reuses the dc4_mv images just computed above, so this is free, and is + # deliberately independent of dc_threshold_mv (see that function's + # docstring) so a threshold that happens to leave some angle's binary + # mask empty can't silently degrade the pivot back to the bbox center. + pivot_mm = compute_pivot_points_mm(sras, dc4_mv=dc4_mv) # ---- Step 2: coarse correlation stage (downsample first, then rotate) # Downsampling before affine_transform (not after) is what keeps this @@ -852,17 +873,20 @@ def build_manual_alignment(sras: SrasFile, ref_angle_idx: int, never be transformed. Pass an already-computed *pivot_mm* (e.g. ManualAlignmentDialog's own - cache, built once from its live CH4 masks) to skip recomputing every + cache, built once from its live CH4 images) to skip recomputing every angle's DC image here; otherwise it's computed fresh via compute_pivot_points_mm, which is still fine for a one-off Save or sidecar restore (just not free on a very large, not-yet-cached scan). + dc_threshold_mv itself plays no part in the pivot (see that function's + docstring) — it's stored on the returned AlignmentResult purely as a + record of the RF-mask threshold in effect at the time. """ n = sras.n_angles dx_ref, dy_ref = _pixel_pitch_mm(sras, ref_angle_idx) params = {a: per_angle_params.get(a, ManualAngleParams()) for a in range(n)} params[ref_angle_idx] = ManualAngleParams() if pivot_mm is None: - pivot_mm = compute_pivot_points_mm(sras, dc_threshold_mv) + pivot_mm = compute_pivot_points_mm(sras) canvas_origin_mm, canvas_shape = union_canvas_mm( sras, ref_angle_idx, dx_ref, dy_ref, params, pivot_mm, margin_frac=0.0) @@ -904,15 +928,26 @@ def sidecar_path(sras_path) -> Path: return p.with_name(p.name + ".align.json") +# Bumped from 1 -> 2 when the rotation pivot changed from the raw scan- +# window bbox center to a content-derived centroid, *and* _theta_deg's sign +# convention was corrected — either change alone makes a version-1 file's +# stored rotation_deg/shift_mm numbers describe a different (and, for the +# pivot bug, actively wrong) transform than they would today. Loading one +# unchanged would silently reproduce exactly the "scans show up everywhere" +# symptom these fixes address, so version-1 sidecars are treated as absent +# rather than migrated. +_SIDECAR_SCHEMA_VERSION = 2 + + def save_manual_alignment(sras: SrasFile, ref_angle_idx: int, dc_threshold_mv: float, per_angle: dict[int, ManualAngleParams]) -> Path: """Write the sidecar JSON for sras.path (overwriting any existing one) and return the path written. - Schema (schema_version 1): + Schema (schema_version 2): { - "schema_version": 1, + "schema_version": 2, "ref_angle_idx": , "dc_threshold_mv": , "per_angle": { @@ -925,7 +960,7 @@ def save_manual_alignment(sras: SrasFile, ref_angle_idx: int, """ path = sidecar_path(sras.path) payload = { - "schema_version": 1, + "schema_version": _SIDECAR_SCHEMA_VERSION, "ref_angle_idx": ref_angle_idx, "dc_threshold_mv": dc_threshold_mv, "per_angle": { @@ -949,7 +984,7 @@ def load_manual_alignment(sras: SrasFile) -> ManualAlignmentSidecar | None: return None try: raw = json.loads(path.read_text()) - if raw.get("schema_version") != 1: + if raw.get("schema_version") != _SIDECAR_SCHEMA_VERSION: return None per_angle = { int(a): ManualAngleParams( diff --git a/sras_viewer.py b/sras_viewer.py index 2c26b71..ea9a663 100644 --- a/sras_viewer.py +++ b/sras_viewer.py @@ -1064,7 +1064,12 @@ class ManualAlignmentDialog(QDialog): max_dim = max(max(img.shape) for img in self._dc4_mv.values()) self._downsample_factor = max(1, int(np.ceil(max_dim / self._MAX_PREVIEW_DIM))) self._recompute_masks_small() - self._recompute_pivot_mm() + # Alignment pivot: the CH4-signal-weighted centroid of each angle's + # own footprint (see compute.compute_pivot_points_mm) — computed once + # from the full-res CH4 images and deliberately independent of the + # mask threshold, so it never needs recomputing when that changes + # (unlike _masks_small, which is purely for the overlay's visuals). + self._pivot_mm = compute.compute_pivot_points_mm(self._sras, self._dc4_mv) self._rebuild_preview_canvas() self._set_controls_enabled(True) self.lbl_status.setText("Ready.") @@ -1072,7 +1077,9 @@ class ManualAlignmentDialog(QDialog): def _recompute_masks_small(self): """Threshold + downsample every angle's already-in-memory full-res CH4 mV image. Cheap (a compare + block-mean), so this re-runs in - full whenever the mask-threshold spin box changes — no re-fetch.""" + full whenever the mask-threshold spin box changes — no re-fetch. + Purely for the overlay's visuals — the alignment pivot does not + depend on this threshold (see _pivot_mm / compute_pivot_points_mm).""" threshold = self.spin_mask_threshold_mv.value() factor = self._downsample_factor self._masks_small = { @@ -1081,19 +1088,6 @@ class ManualAlignmentDialog(QDialog): for a, img in self._dc4_mv.items() } - def _recompute_pivot_mm(self): - """Each angle's alignment pivot: the centroid of its own full- - resolution CH4-threshold mask (its own sample footprint), not the - raw scan-window bbox center — see compute.compute_pivot_points_mm. - Recomputed alongside _recompute_masks_small whenever the mask - threshold changes, from the full-res masks (not the downsampled - preview ones) since this feeds the canonical geometry, not just the - preview.""" - threshold = self.spin_mask_threshold_mv.value() - masks = {a: img >= threshold for a, img in self._dc4_mv.items()} - self._pivot_mm = compute.compute_pivot_points_mm( - self._sras, threshold, masks=masks) - # ------------------------------------------------------------------ # Preview canvas: full rebuild vs. incremental single-layer refresh # ------------------------------------------------------------------ @@ -1227,7 +1221,6 @@ class ManualAlignmentDialog(QDialog): if not self._masks_ready: return self._recompute_masks_small() - self._recompute_pivot_mm() self._rebuild_preview_canvas() # ------------------------------------------------------------------ diff --git a/tools/test_gui.py b/tools/test_gui.py index 8db05f6..83a0366 100644 --- a/tools/test_gui.py +++ b/tools/test_gui.py @@ -240,27 +240,36 @@ def main(): print("\nmanual alignment (Fusion)") check("manual alignment action enabled", win._manual_align_act.isEnabled()) - # --- Alignment pivot is the mask centroid, not the raw bbox center ----- - corner_mask = np.zeros(s.image_shape(0), dtype=bool) - corner_mask[0, 0] = True + # --- Alignment pivot is a signal-weighted centroid, not the raw bbox -- + # center, and is independent of any DC threshold (so a threshold that + # happens to leave a real angle's binary mask empty can't silently + # degrade the pivot back to the bbox center). + corner_signal = np.zeros(s.image_shape(0), dtype=np.float32) + corner_signal[0, 0] = 1.0 # single spike -> weighted centroid is exact expected_corner = (float(s.x_axis_mm(0)[0]), float(s.y_positions_mm(0)[0])) - centroid = compute._mask_centroid_mm(s, 0, corner_mask) - check("mask centroid of a single corner pixel is that corner exactly", + centroid = compute._signal_centroid_mm(s, 0, corner_signal) + check("signal-weighted centroid of a single spike pixel is that pixel exactly", np.allclose(centroid, expected_corner), f"{centroid} vs {expected_corner}") bbox_center = compute._bbox_center_mm(s, 0) - check("mask centroid differs from the raw scan-window bbox center", + check("signal centroid differs from the raw scan-window bbox center", not np.allclose(centroid, bbox_center), f"centroid {centroid} vs bbox center {bbox_center}") - # A threshold that would yield an *empty* mask if recomputed from scratch - # (real DC values never reach 999 mV) — so this only matches expected_corner - # if compute_pivot_points_mm actually reused the pre-computed masks dict - # instead of silently recomputing (and falling back to the bbox center). - reused_pivot = compute.compute_pivot_points_mm( - s, 999.0, masks={0: corner_mask.astype(np.float32)})[0] - check("compute_pivot_points_mm reuses a pre-computed masks dict", + # compute_pivot_points_mm should reuse a pre-computed dc4_mv dict rather + # than recomputing from the real DC4 image (which has no such spike and + # would give a different answer if silently recomputed). + reused_pivot = compute.compute_pivot_points_mm(s, dc4_mv={0: corner_signal})[0] + check("compute_pivot_points_mm reuses a pre-computed dc4_mv dict", np.allclose(reused_pivot, expected_corner)) + # A perfectly flat signal carries no information to weight by, so it + # falls back to the bbox center rather than producing a NaN/degenerate + # centroid. + flat_signal = np.full(s.image_shape(0), 5.0, dtype=np.float32) + flat_centroid = compute._signal_centroid_mm(s, 0, flat_signal) + check("a perfectly flat signal falls back to the bbox center", + np.allclose(flat_centroid, bbox_center)) + # --- Rotation sign convention: negative of the raw angles_deg delta ---- check("_theta_deg negates the raw angles_deg delta (GR stage's positive " "angle is the opposite rotational sense from this module's CCW " @@ -339,7 +348,8 @@ def main(): sidecar = compute.sidecar_path(s.path) check("sidecar file written", sidecar.exists()) raw = json.loads(sidecar.read_text()) if sidecar.exists() else {} - check("sidecar schema_version is 1", raw.get("schema_version") == 1) + check("sidecar schema_version is current", + raw.get("schema_version") == compute._SIDECAR_SCHEMA_VERSION) check("sidecar per_angle round-trips the dialog's resolved params", all(raw.get("per_angle", {}).get(str(a), {}).get("rotation_deg") == dlg._angle_params[a].rotation_deg for a in range(s.n_angles))) @@ -350,6 +360,14 @@ def main(): check("Aligned View auto-enabled after Save", win.chk_aligned_view.isEnabled() and win.chk_aligned_view.isChecked()) + # --- An old-schema sidecar (pre-pivot/sign fix) is treated as absent ------ + stale = dict(raw) + stale["schema_version"] = compute._SIDECAR_SCHEMA_VERSION - 1 + sidecar.write_text(json.dumps(stale)) + check("a sidecar with an old schema_version is not loaded", + compute.load_manual_alignment(s) is None) + sidecar.write_text(json.dumps(raw)) # restore for the rest of this section + # --- Clear (with confirmation) -------------------------------------------- with patch("sras_viewer.QMessageBox.question", return_value=QMessageBox.StandardButton.Yes):