From ca0c736c28013b249f537ccf705b8e3007fdf3b7 Mon Sep 17 00:00:00 2001 From: Thomas Ales Date: Fri, 31 Jul 2026 09:36:12 -0500 Subject: [PATCH] Fix alignment rotation pivot and sign convention Rotation was pivoting on each angle's raw scan-window bbox center, which ties alignment to wherever the window happened to sit in microscope/global XY space rather than to the sample itself, leaving a residual orbital drift between angles. Pivot on each angle's own CH4-threshold mask centroid instead (its own sample footprint), for both the automatic phase- correlation path and the manual dialog's Auto De-rotate/live preview. Also negate the rotation used everywhere (_theta_deg): the GR stage's reported angle increases in the opposite rotational sense from this module's CCW math convention, so de-rotating by the raw angles_deg delta was making misalignment worse rather than better. Co-Authored-By: Claude Sonnet 5 --- sras_compute.py | 158 +++++++++++++++++++++++++++++++++++----------- sras_viewer.py | 25 ++++++-- tools/test_gui.py | 29 +++++++++ 3 files changed, 172 insertions(+), 40 deletions(-) diff --git a/sras_compute.py b/sras_compute.py index 37f7b16..803d2c1 100644 --- a/sras_compute.py +++ b/sras_compute.py @@ -442,14 +442,75 @@ def _rotation_matrix(theta_deg: float) -> np.ndarray: def _theta_deg(sras: SrasFile, angle_idx: int, ref_idx: int) -> float: - return float(sras.angles_deg[angle_idx] - sras.angles_deg[ref_idx]) + """CCW rotation, in degrees and in _rotation_matrix's convention, that + maps angle_idx's own local mm frame onto ref_idx's. + + This is the *negative* of the raw angles_deg delta: the GR rotation + stage's reported angle increases in the opposite rotational sense from + this module's math-positive (CCW, x toward y) convention in scan mm + space. Rotating by +(angles_deg[a] - angles_deg[ref]) therefore turns + misalignment the wrong way — confirmed empirically (Auto De-rotate made + real scans worse, not better, before this negation). + """ + 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(): + 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()) + + +def compute_pivot_points_mm(sras: SrasFile, dc_threshold_mv: float, + masks: 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. + + Pivoting on each angle's own content — instead of on wherever its + scanned window happened to sit in microscope/global XY space — is what + makes alignment purely relative *between scans* rather than to a global + coordinate system: the rotation stage's true mechanical axis need not + coincide with the scan window's geometric center, and the sample need + not be perfectly centered on that axis either, so a bbox-center pivot + 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. + """ + masks = masks 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) + return pivots def _corners_in_ref_frame(sras: SrasFile, angle_idx: int, ref_idx: int, + pivot_mm: dict[int, tuple[float, float]], shift_mm=(0.0, 0.0), theta_deg: float | None = None) -> np.ndarray: - """Angle *angle_idx*'s bbox corners, rotated about its own centroid into - the reference frame and translated by *shift_mm*. Shape (4, 2). + """Angle *angle_idx*'s bbox corners, rotated about its own alignment + pivot (pivot_mm[angle_idx] — see compute_pivot_points_mm) into the + reference frame and translated by *shift_mm*. Shape (4, 2). theta_deg overrides the analytic angles_deg-derived rotation used by default — the manual-alignment path (union_canvas_mm) passes a @@ -458,8 +519,8 @@ def _corners_in_ref_frame(sras: SrasFile, angle_idx: int, ref_idx: int, """ theta = _theta_deg(sras, angle_idx, ref_idx) if theta_deg is None else theta_deg R = _rotation_matrix(theta) - c_a = np.array(_bbox_center_mm(sras, angle_idx)) - c_ref = np.array(_bbox_center_mm(sras, ref_idx)) + c_a = np.array(pivot_mm[angle_idx]) + c_ref = np.array(pivot_mm[ref_idx]) shift = np.asarray(shift_mm, dtype=np.float64) return np.array([R @ (corner - c_a) + c_ref + shift for corner in _bbox_corners_mm(sras, angle_idx)]) @@ -469,6 +530,7 @@ def _build_affine_canvas_to_raw(sras: SrasFile, angle_idx: int, ref_idx: int, shift_mm: tuple[float, float], canvas_dx: float, canvas_dy: float, canvas_origin_mm: tuple[float, float], + pivot_mm: dict[int, tuple[float, float]], theta_deg: float | None = None ) -> tuple[np.ndarray, np.ndarray]: """matrix, offset s.t. raw_index = matrix @ [row_out, col_out] + offset, @@ -479,18 +541,21 @@ def _build_affine_canvas_to_raw(sras: SrasFile, angle_idx: int, ref_idx: int, [lx;ly] = R(theta)^T @ ([X;Y]-c_ref-shift) + c_a # undo rotation+shift -> angle a's local mm [row;col] = D @ ([lx;ly] - [x_start_a; y0_a]) # local mm -> angle a's raw idx - where theta = angles_deg[angle_idx] - angles_deg[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 raw-bbox mm - centroid (the rotation pivot — this keeps rotated content centered, - minimizing required canvas padding), and A_out/D are the index<->mm - scaling matrices for the canvas pitch and this angle's own native pitch - respectively. + 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. """ theta = _theta_deg(sras, angle_idx, ref_idx) if theta_deg is None else theta_deg Rinv = _rotation_matrix(theta).T - cx_a, cy_a = _bbox_center_mm(sras, angle_idx) - cx_ref, cy_ref = _bbox_center_mm(sras, ref_idx) + cx_a, cy_a = pivot_mm[angle_idx] + cx_ref, cy_ref = pivot_mm[ref_idx] dx_a, dy_a = _pixel_pitch_mm(sras, angle_idx) x0_a = float(sras.x_start_mm[angle_idx]) y0_a = float(sras.y_positions_mm(angle_idx)[0]) @@ -554,15 +619,17 @@ def _phase_correlate_shift(ref_img: np.ndarray, mov_img: np.ndarray) -> tuple[in def _working_canvas_for_pair(sras: SrasFile, ref_idx: int, a_idx: int, - dx: float, dy: float, margin_frac: float = 0.3 + dx: float, dy: float, + pivot_mm: dict[int, tuple[float, float]], + margin_frac: float = 0.3 ) -> tuple[tuple[float, float], tuple[int, int]]: """Union of the reference's own raw bbox and angle a's raw bbox rotated - (about its own center) into the ref frame with zero shift, padded by - margin_frac on each side — sized generously so the true phase- - correlation shift lands well inside the canvas (see + (about its own alignment pivot) into the ref frame with zero shift, + padded by margin_frac on each side — sized generously so the true + phase-correlation shift lands well inside the canvas (see _phase_correlate_shift's wraparound note).""" pts = np.vstack([_bbox_corners_mm(sras, ref_idx), - _corners_in_ref_frame(sras, a_idx, ref_idx)]) + _corners_in_ref_frame(sras, a_idx, ref_idx, pivot_mm)]) x_min, y_min = pts.min(axis=0) x_max, y_max = pts.max(axis=0) pad_x, pad_y = (x_max - x_min) * margin_frac, (y_max - y_min) * margin_frac @@ -615,6 +682,11 @@ def compute_angle_alignment(sras: SrasFile, ref_angle_idx: int, masks = dict(enumerate(_parallel_map(mask_for, range(n), n_workers))) + # 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) + # ---- Step 2: coarse correlation stage (downsample first, then rotate) # Downsampling before affine_transform (not after) is what keeps this # tractable for a v6 scan with thousands of rows/frames per angle. @@ -630,15 +702,16 @@ def compute_angle_alignment(sras: SrasFile, ref_angle_idx: int, tick(25, 50) return (0.0, 0.0) work_origin, work_shape = _working_canvas_for_pair( - sras, ref_angle_idx, a, dx_c, dy_c, margin_frac=0.3) + sras, ref_angle_idx, a, dx_c, dy_c, pivot_mm, margin_frac=0.3) # Coarse canvas->raw affine at native pitch, then rescale by /factor # so it maps coarse-canvas idx -> coarse (downsampled) raw idx — # exact for block-mean downsampling (up to the trimmed remainder). m_a, o_a = _build_affine_canvas_to_raw( - sras, a, ref_angle_idx, (0.0, 0.0), dx_c, dy_c, work_origin) + sras, a, ref_angle_idx, (0.0, 0.0), dx_c, dy_c, work_origin, pivot_mm) m_ref, o_ref = _build_affine_canvas_to_raw( - sras, ref_angle_idx, ref_angle_idx, (0.0, 0.0), dx_c, dy_c, work_origin) + sras, ref_angle_idx, ref_angle_idx, (0.0, 0.0), dx_c, dy_c, work_origin, + pivot_mm) rotated_a = scipy_ndimage.affine_transform( masks_small[a], m_a / factor, offset=o_a / factor, output_shape=work_shape, order=0, mode="constant", cval=0.0) @@ -653,8 +726,9 @@ def compute_angle_alignment(sras: SrasFile, ref_angle_idx: int, shifts_mm = dict(enumerate(_parallel_map(shift_for, range(n), n_workers))) # ---- Step 3: union bounding box over all angles (rotation+shift applied) - corners = np.vstack([_corners_in_ref_frame(sras, a, ref_angle_idx, shifts_mm[a]) - for a in range(n)]) + corners = np.vstack([ + _corners_in_ref_frame(sras, a, ref_angle_idx, pivot_mm, shifts_mm[a]) + for a in range(n)]) x_min, y_min = corners.min(axis=0) x_max, y_max = corners.max(axis=0) n_cols = int(np.ceil((x_max - x_min) / dx_ref)) + 1 @@ -665,7 +739,8 @@ def compute_angle_alignment(sras: SrasFile, ref_angle_idx: int, per_angle: dict[int, AngleTransform] = {} for a in range(n): matrix, offset = _build_affine_canvas_to_raw( - sras, a, ref_angle_idx, shifts_mm[a], dx_ref, dy_ref, canvas_origin_mm) + sras, a, ref_angle_idx, shifts_mm[a], dx_ref, dy_ref, canvas_origin_mm, + pivot_mm) per_angle[a] = AngleTransform( _theta_deg(sras, a, ref_angle_idx), shifts_mm[a], matrix, offset) if progress_cb: @@ -698,6 +773,7 @@ _compute_angle_alignment = compute_angle_alignment def union_canvas_mm(sras: SrasFile, ref_angle_idx: int, dx: float, dy: float, per_angle_params: dict[int, ManualAngleParams], + pivot_mm: dict[int, tuple[float, float]], margin_frac: float = 0.0 ) -> tuple[tuple[float, float], tuple[int, int]]: """Shared-canvas origin (mm) and (n_rows, n_cols) at pitch (dx, dy) that @@ -720,7 +796,7 @@ def union_canvas_mm(sras: SrasFile, ref_angle_idx: int, dx: float, dy: float, n = sras.n_angles corners = np.vstack([ _corners_in_ref_frame( - sras, a, ref_angle_idx, + sras, a, ref_angle_idx, pivot_mm, per_angle_params.get(a, ManualAngleParams()).shift_mm, theta_deg=per_angle_params.get(a, ManualAngleParams()).rotation_deg) for a in range(n)]) @@ -741,7 +817,8 @@ def reproject_mask(sras: SrasFile, angle_idx: int, ref_angle_idx: int, shift_mm: tuple[float, float], canvas_dx: float, canvas_dy: float, canvas_origin_mm: tuple[float, float], - canvas_shape: tuple[int, int]) -> np.ndarray: + canvas_shape: tuple[int, int], + pivot_mm: dict[int, tuple[float, float]]) -> np.ndarray: """Resample one angle's binary/float mask onto an arbitrary canvas via an explicit rotation+shift — the single building block ManualAlignmentDialog's live preview repeatedly calls (once per @@ -753,7 +830,7 @@ def reproject_mask(sras: SrasFile, angle_idx: int, ref_angle_idx: int, """ matrix, offset = _build_affine_canvas_to_raw( sras, angle_idx, ref_angle_idx, shift_mm, canvas_dx, canvas_dy, - canvas_origin_mm, theta_deg=rotation_deg) + canvas_origin_mm, pivot_mm, theta_deg=rotation_deg) return scipy_ndimage.affine_transform( mask.astype(np.float32, copy=False), matrix, offset=offset, output_shape=canvas_shape, order=0, mode="constant", cval=0.0) @@ -761,32 +838,41 @@ def reproject_mask(sras: SrasFile, angle_idx: int, ref_angle_idx: int, def build_manual_alignment(sras: SrasFile, ref_angle_idx: int, dc_threshold_mv: float, - per_angle_params: dict[int, ManualAngleParams] + per_angle_params: dict[int, ManualAngleParams], + pivot_mm: dict[int, tuple[float, float]] | None = None ) -> AlignmentResult: """Build a full, full-resolution AlignmentResult from user-supplied per-angle rotation+shift — the Manual Alignment counterpart to compute_angle_alignment, skipping its mask/phase-correlation search entirely (every angle's transform here is exactly what the caller - supplied). Pure matrix/bbox math, no image data touched anywhere in this - function, so it is cheap enough to call synchronously on the GUI thread. - The reference angle's params are always forced to identity, regardless - of what per_angle_params holds for it — it defines the shared origin and - must never be transformed. + supplied). Pure matrix/bbox math once pivot_mm is known, so it is cheap + enough to call synchronously on the GUI thread. The reference angle's + params are always forced to identity, regardless of what + per_angle_params holds for it — it defines the shared origin and must + 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 + 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). """ 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) canvas_origin_mm, canvas_shape = union_canvas_mm( - sras, ref_angle_idx, dx_ref, dy_ref, params, margin_frac=0.0) + sras, ref_angle_idx, dx_ref, dy_ref, params, pivot_mm, margin_frac=0.0) per_angle: dict[int, AngleTransform] = {} for a in range(n): p = params[a] matrix, offset = _build_affine_canvas_to_raw( sras, a, ref_angle_idx, p.shift_mm, dx_ref, dy_ref, - canvas_origin_mm, theta_deg=p.rotation_deg) + canvas_origin_mm, pivot_mm, theta_deg=p.rotation_deg) per_angle[a] = AngleTransform(p.rotation_deg, p.shift_mm, matrix, offset) return AlignmentResult(ref_angle_idx, dc_threshold_mv, canvas_shape, diff --git a/sras_viewer.py b/sras_viewer.py index e6f2730..2c26b71 100644 --- a/sras_viewer.py +++ b/sras_viewer.py @@ -860,6 +860,7 @@ class ManualAlignmentDialog(QDialog): self._downsample_factor = 1 self._dc4_mv: dict[int, np.ndarray] = {} self._masks_small: dict[int, np.ndarray] = {} + self._pivot_mm: dict[int, tuple[float, float]] = {} self._preview_layers: dict[int, np.ndarray] = {} self._preview_origin_mm = (0.0, 0.0) self._preview_shape = (1, 1) @@ -1063,6 +1064,7 @@ 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() self._rebuild_preview_canvas() self._set_controls_enabled(True) self.lbl_status.setText("Ready.") @@ -1079,6 +1081,19 @@ 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 # ------------------------------------------------------------------ @@ -1096,14 +1111,14 @@ class ManualAlignmentDialog(QDialog): dx_c, dy_c = dx_ref * factor, dy_ref * factor origin, shape = compute.union_canvas_mm( self._sras, self._ref_angle_idx, dx_c, dy_c, self._angle_params, - margin_frac=self._PREVIEW_MARGIN_FRAC) + self._pivot_mm, margin_frac=self._PREVIEW_MARGIN_FRAC) self._preview_origin_mm, self._preview_shape = origin, shape self._preview_dx_mm, self._preview_dy_mm = dx_c, dy_c self._preview_layers = { a: compute.reproject_mask( self._sras, a, self._ref_angle_idx, self._masks_small[a], self._angle_params[a].rotation_deg, self._angle_params[a].shift_mm, - dx_c, dy_c, origin, shape) + dx_c, dy_c, origin, shape, self._pivot_mm) for a in range(self._sras.n_angles) } self._redraw_overlay() @@ -1117,7 +1132,7 @@ class ManualAlignmentDialog(QDialog): self._sras, a, self._ref_angle_idx, self._masks_small[a], self._angle_params[a].rotation_deg, self._angle_params[a].shift_mm, self._preview_dx_mm, self._preview_dy_mm, - self._preview_origin_mm, self._preview_shape) + self._preview_origin_mm, self._preview_shape, self._pivot_mm) self._redraw_overlay() def _redraw_overlay(self): @@ -1212,6 +1227,7 @@ class ManualAlignmentDialog(QDialog): if not self._masks_ready: return self._recompute_masks_small() + self._recompute_pivot_mm() self._rebuild_preview_canvas() # ------------------------------------------------------------------ @@ -1237,7 +1253,8 @@ class ManualAlignmentDialog(QDialog): resolved = dict(self._angle_params) # already concrete floats try: path = save_manual_alignment(self._sras, self._ref_angle_idx, threshold, resolved) - result = build_manual_alignment(self._sras, self._ref_angle_idx, threshold, resolved) + result = build_manual_alignment(self._sras, self._ref_angle_idx, threshold, + resolved, self._pivot_mm) except OSError as exc: QMessageBox.warning(self, "Save Alignment Failed", str(exc)) return diff --git a/tools/test_gui.py b/tools/test_gui.py index 915dbd9..8db05f6 100644 --- a/tools/test_gui.py +++ b/tools/test_gui.py @@ -240,6 +240,35 @@ 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 + 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", + 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", + 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", + np.allclose(reused_pivot, expected_corner)) + + # --- 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 " + "math convention)", + all(np.isclose(compute._theta_deg(s, a, 0), + -(float(s.angles_deg[a]) - float(s.angles_deg[0]))) + for a in range(s.n_angles))) + # --- Open, seeding from the still-live automatic AlignmentResult -------- win._on_manual_alignment() check("dialog opened", win._manual_align_dialog is not None)