diff --git a/changelog.md b/changelog.md index bb5c30b0..73cca0f5 100644 --- a/changelog.md +++ b/changelog.md @@ -1,6 +1,14 @@ # Changelog -Last change: 28-MAY-2026 CEST +Last change: 11-JUN-2026 CEST + +## 0.10.2 +- Rotations use quaterions for unambiguous workflow, fixing bugs [#673](https://github.com/jungmannlab/picasso/issues/673), [#674](https://github.com/jungmannlab/picasso/issues/674) and [#675](https://github.com/jungmannlab/picasso/issues/675) +- All 3 angles in "Rotate by angle" in the 3D rotation window are accumulated into a single widget +- Render by property fixed for large files ([#677](https://github.com/jungmannlab/picasso/issues/677)), possibly related to [#672](https://github.com/jungmannlab/picasso/issues/672) +- Fixed "Best fitting combination" button in SPINNA ([#676](https://github.com/jungmannlab/picasso/issues/676)) +- Fixed 3D calibration when frame range is user-specified +- Fixed redoing 3D calibration when identification parameters change ## 0.10.1 diff --git a/docs/render.rst b/docs/render.rst index 85edc7a2..edf10572 100644 --- a/docs/render.rst +++ b/docs/render.rst @@ -67,8 +67,6 @@ Note that to build animations, the user must have ``ffmpeg`` installed on their Rotation around z-axis is available by pressing Ctrl/Command. Rotation axis can be frozen by pressing x/y/z to freeze around the corresponding axes (to freeze around the z-axis, Ctrl/Command must be pressed as well). -There are several things to keep in mind when using the rotation window. Firstly, using individual localization precision is very slow and is not recommended as a default blur method. Also, the size of the rotation window can be altered, however, if it becomes too large, rendering may start to lag. - RESI ---- .. image:: ../docs/render_resi.png diff --git a/picasso/gui/localize.py b/picasso/gui/localize.py index 33717b15..85d6848c 100644 --- a/picasso/gui/localize.py +++ b/picasso/gui/localize.py @@ -2103,6 +2103,8 @@ def _clean_up_external_ids(self) -> None: self.last_identification_info = { "Box Size": self.parameters_dialog.box_spinbox.value(), "Min. Net Gradient": self.parameters_dialog.mng_slider.value(), + "ROI": self.view.roi, + "Frame bounds": self.frame_range, } self.ready_for_fit = True self.draw_frame() @@ -2551,6 +2553,7 @@ def on_fit_finished( playsound(sound_path, block=False) base, ext = os.path.splitext(self.movie_path) if calibrate_z: + self.parameters_dialog.gpufit_checkbox.setDisabled(False) step, ok = QtWidgets.QInputDialog.getDouble( self, "3D Calibration", @@ -2572,6 +2575,7 @@ def on_fit_finished( step, self.parameters_dialog.magnification_factor.value(), path=path, + frame_bounds=self.frame_range, ) dt = time.time() - t0 if dt > lib.SOUND_NOTIFICATION_DURATION: @@ -2902,11 +2906,34 @@ def localize(self, calibrate_z: bool = False) -> None: Default is False. """ self.parameters_dialog.gpufit_checkbox.setDisabled(True) - if self.identifications is not None and calibrate_z is True: + if ( + calibrate_z + and self.identifications is not None + and self.ready_for_fit + and not self.identifications_outdated() + ): + # reuse existing identifications (e.g., loaded picks/locs) self.fit(calibrate_z=calibrate_z) else: self.identify(fit_afterwards=True, calibrate_z=calibrate_z) + def identifications_outdated(self) -> bool: + """Check whether the identification settings (box size, min. net + gradient, ROI, frame range) have changed since + ``self.identifications`` was created.""" + last = self.last_identification_info + if last is None: + return True + if any( + last.get(key) != value for key, value in self.parameters.items() + ): + return True + if last.get("ROI") != self.view.roi: + return True + if last.get("Frame bounds") != self.frame_range: + return True + return False + def select_locs_columns(self) -> None: """Select only the columns that are checked in the corresponding dialog.""" diff --git a/picasso/gui/render.py b/picasso/gui/render.py index 640bc868..b985aded 100644 --- a/picasso/gui/render.py +++ b/picasso/gui/render.py @@ -55,8 +55,8 @@ ) from .rotation import RotationWindow -# PyImarisWrite works on windows only -from ..ext.bitplane import IMSWRITER +# Optional modules with external/hardware dependencies live in ext +from ..ext.bitplane import IMSWRITER # PyImarisWrite works on Windows only if IMSWRITER: from ..ext.bitplane import numpy_to_imaris @@ -10995,14 +10995,6 @@ def show_drift(self) -> None: self.plot_window.plot(drift) self.plot_window.show() - def sync_groups(self) -> None: - """Remove localizations whose group field is not found in all - channels.""" - if len(self.locs_paths) < 2: - return - self.locs = lib.sync_groups(self.locs) - self.update_scene(resample_locs=True) - def undrift_aim(self) -> None: """Undrift with Adaptive Intersection Maximization (AIM). @@ -11252,6 +11244,14 @@ def _apply_drift(self, channel: int, drift: pd.DataFrame) -> None: self.render_index[channel] = None self.update_scene(resample_locs=True) + def sync_groups(self) -> None: + """Remove localizations whose group field is not found in all + channels.""" + if len(self.locs_paths) < 2: + return + self.locs = lib.sync_groups(self.locs) + self.update_scene(resample_locs=True) + def unfold_groups_square(self) -> None: """Shifts grouped localizations onto a square grid with a chosen number of columns and spacing. Localizations can be grouped, for @@ -11995,6 +11995,7 @@ def initUI(self, plugins_loaded: bool) -> None: # menu bar - Postprocess postprocess_menu = self.menu_bar.addMenu("Postprocess") + undrift_aim_action = postprocess_menu.addAction("Undrift by AIM") undrift_aim_action.setShortcut("Ctrl+U") undrift_aim_action.triggered.connect(self.view.undrift_aim) @@ -12011,6 +12012,7 @@ def initUI(self, plugins_loaded: bool) -> None: undrift_from_picked2d_action.triggered.connect( self.view.undrift_from_picked2d ) + undrift_action = postprocess_menu.addAction("Undrift by RCC") undrift_action.triggered.connect(self.view.undrift_rcc) drift_action = postprocess_menu.addAction("Undo drift") @@ -12829,9 +12831,9 @@ def open_rotated_locs(self) -> None: self.tools_settings_dialog.pick_side_length.setValue( self.view.infos[0][-1]["Pick size (nm)"] ) - self.window_rot.view_rot.angx = self.view.infos[0][-1]["angx"] - self.window_rot.view_rot.angy = self.view.infos[0][-1]["angy"] - self.window_rot.view_rot.angz = self.view.infos[0][-1]["angz"] + self.window_rot.view_rot.load_saved_rotation( + self.view.infos[0][-1] + ) self.rot_win() def resize_view_to_fov(self, w: float, h: float) -> None: diff --git a/picasso/gui/rotation.py b/picasso/gui/rotation.py index f7a59ce7..fc1bb310 100644 --- a/picasso/gui/rotation.py +++ b/picasso/gui/rotation.py @@ -311,9 +311,14 @@ class AnimationDialog(lib.Dialog): Click to delete the last saved position. fps : QSpinBox Contains frames per second used in the animation. - positions : list - Contains all positions in the animation sequence; each includes - 3 rotations angles and viewport. + positions : list of dict + Contains all positions in the animation sequence. Each holds + the rotation (scipy Rotation, key ``"R"``), the accumulated + rotation vector (key ``"rotvec"``, radians, scipy convention, + keeps full turns), the rotation accumulated since the previous + position (key ``"segment_rotvec"``, radians, scipy convention; + may encode turns beyond 180 degrees) and the viewport (key + ``"viewport"``). rows : list of dict One entry per saved position, holding the row's widgets: ``p_label``, ``angle_label``, ``show_btn``, ``d_label``, @@ -352,7 +357,10 @@ def __init__(self, window: QtWidgets.QMainWindow) -> None: # Header: current position header = QtWidgets.QHBoxLayout() cp_label = QtWidgets.QLabel("Current position:") - cp_label.setToolTip("Current rotation angles in x, y, z (deg).") + cp_label.setToolTip( + "Current rotation in x, y, z (deg). The angles keep track " + "of full turns, e.g. 720 encodes two full rotations." + ) header.addWidget(cp_label) angx = np.round(self.window.view_rot.angx * 180 / np.pi, 1) angy = np.round(self.window.view_rot.angy * 180 / np.pi, 1) @@ -453,25 +461,35 @@ def add_position(self, freeze: bool = False) -> None: True when the new position is the same as the last one, i.e., when self.stay is clicked. """ - # check that the viewport or angle(s) have changed - cond1 = cond2 = cond3 = False + view = self.window.view_rot + # rotation accumulated since the last position (keeps full + # turns, e.g. 720 deg); zero when freezing in place + segment_rotvec = ( + np.zeros(3) if freeze else view.rotation_since_anchor() + ) + + # check that the viewport or the rotation have changed if not freeze and self.positions: - cond1 = self.window.view_rot.angx == self.positions[-1][0] - cond2 = self.window.view_rot.angy == self.positions[-1][1] - cond3 = self.window.view_rot.angz == self.positions[-1][2] - cond4 = self.window.view_rot.viewport == self.positions[-1][3] - if all([cond1, cond2, cond3, cond4]): + last = self.positions[-1] + same_rotation = ( + view.rotation * last["R"].inv() + ).magnitude() < 1e-9 + no_turns = np.linalg.norm(segment_rotvec) < 1e-9 + same_viewport = view.viewport == last["viewport"] + if same_rotation and no_turns and same_viewport: return # add a new position to the attribute self.positions.append( - [ - self.window.view_rot.angx, - self.window.view_rot.angy, - self.window.view_rot.angz, - self.window.view_rot.viewport, - ] + { + "R": view.rotation, + "rotvec": view._rotvec.copy(), + "segment_rotvec": segment_rotvec, + "viewport": view.viewport, + } ) + # track the next segment's rotation from this position onwards + view.reset_rotation_anchor() index = len(self.positions) - 1 grid_row = index @@ -483,9 +501,9 @@ def add_position(self, freeze: bool = False) -> None: ) self.rows_layout.addWidget(p_label, grid_row, 0) - angx = np.round(self.window.view_rot.angx * 180 / np.pi, 1) - angy = np.round(self.window.view_rot.angy * 180 / np.pi, 1) - angz = np.round(self.window.view_rot.angz * 180 / np.pi, 1) + angx = np.round(view.angx * 180 / np.pi, 1) + angy = np.round(view.angy * 180 / np.pi, 1) + angz = np.round(view.angz * 180 / np.pi, 1) angle_label = QtWidgets.QLabel("{}, {}, {}".format(angx, angy, angz)) self.rows_layout.addWidget(angle_label, grid_row, 1) @@ -507,14 +525,12 @@ def add_position(self, freeze: bool = False) -> None: duration.setDecimals(2) self.rows_layout.addWidget(duration, grid_row, 4) - # calculate recommended duration - if not freeze and not all([cond1, cond2, cond3]): - dx = self.positions[-1][0] - self.positions[-2][0] - dy = self.positions[-1][1] - self.positions[-2][1] - dz = self.positions[-1][2] - self.positions[-2][2] - dmax = np.max(np.abs([dx, dy, dz])) + # calculate recommended duration from the total rotation + # (including full turns) at the requested rotation speed + total_angle = float(np.linalg.norm(segment_rotvec)) + if not freeze and total_angle > 1e-9: rot_speed = self.rot_speed.value() * np.pi / 180 - duration.setValue(dmax / rot_speed) + duration.setValue(total_angle / rot_speed) self.rows.append( { @@ -553,10 +569,11 @@ def retrieve_position(self, i: int) -> None: Index of the position to be displayed. """ if i <= len(self.positions) - 1: - self.window.view_rot.angx = self.positions[i][0] - self.window.view_rot.angy = self.positions[i][1] - self.window.view_rot.angz = self.positions[i][2] - self.window.view_rot.update_scene(viewport=self.positions[i][3]) + position = self.positions[i] + self.window.view_rot.set_rotation( + position["R"], rotvec=position["rotvec"] + ) + self.window.view_rot.update_scene(viewport=position["viewport"]) def build_animation(self) -> None: """Create an animation as an .mp4 file using the positions from @@ -592,12 +609,17 @@ def build_animation(self) -> None: ) adjust_display_pixel = disp_dlg.dynamic_disp_px.isChecked() intensities = self.window.window.view.read_relative_intensities() + positions = [(p["R"], p["viewport"]) for p in self.positions] + segment_rotations = [ + p["segment_rotvec"] for p in self.positions[1:] + ] render.build_animation( path, locs, infos, - positions=self.positions, + positions=positions, durations=durations, + segment_rotations=segment_rotations, disp_px_size=disp_dlg.disp_px_size.value(), image_size=( self.window.view_rot.width(), @@ -619,6 +641,62 @@ def build_animation(self) -> None: progress.close() +class RotateByAngleDialog(lib.Dialog): + """Choose rotations angles. + + ... + + Attributes + ---------- + angx, angy, angz: QtWidgets.QDoubleSpinBoxes + Store the rotation angles input by the user. + """ + + def __init__(self, window: QtWidgets.QMainWindow) -> None: + super().__init__(window) + self.window = window + self.setWindowTitle(f"Enter rotation angles") + layout = QtWidgets.QFormLayout(self) + self.angx = QtWidgets.QDoubleSpinBox() + self.angx.setValue(0) + self.angx.setRange(-999999, 999999) + self.angx.setSingleStep(1) + layout.addRow("Angle x (deg)", self.angx) + self.angy = QtWidgets.QDoubleSpinBox() + self.angy.setValue(0) + self.angy.setRange(-999999, 999999) + self.angy.setSingleStep(1) + layout.addRow("Angle y (deg)", self.angy) + self.angz = QtWidgets.QDoubleSpinBox() + self.angz.setValue(0) + self.angz.setRange(-999999, 999999) + self.angz.setSingleStep(1) + layout.addRow("Angle z (deg)", self.angz) + + # OK and Cancel buttons + self.buttons = QtWidgets.QDialogButtonBox( + QtWidgets.QDialogButtonBox.StandardButton.Ok + | QtWidgets.QDialogButtonBox.StandardButton.Cancel, + QtCore.Qt.Orientation.Horizontal, + self, + ) + layout.addRow(self.buttons) + self.buttons.accepted.connect(self.accept) + self.buttons.rejected.connect(self.reject) + + @staticmethod + def getParams(parent: QtWidgets.QMainWindow | None = None) -> tuple: + """Create the dialog and return the requested rot. angles.""" + dialog = RotateByAngleDialog(parent) + result = dialog.exec() + return ( + dialog.angx.value(), + dialog.angy.value(), + dialog.angz.value(), + result == QtWidgets.QDialog.DialogCode.Accepted, + ) + + class ViewRotation(QtWidgets.QLabel): """Display rotated super-resolution datasets. @@ -629,8 +707,10 @@ class ViewRotation(QtWidgets.QLabel): Attributes ---------- angx, angy, angz : float - Current rotation angle around x, y, and z axes (read/write - properties; backed by ``self._R``). + Accumulated rotation around the x, y and z axes (radians) in + the codebase's angle convention (read-only properties; backed + by ``self._rotvec``). They keep track of full turns, e.g. read + 720 degrees after two full rotations. block_x, block_y, block_z : bool True if rotate only around x, y, or z axis respectively. group_color : lib.IntArray1D @@ -660,8 +740,21 @@ class ViewRotation(QtWidgets.QLabel): qimage : QImage Current image of rendered locs, picks and other drawings. _R : scipy.spatial.transform.Rotation - Source of truth for the current rotation. ``angx/angy/angz`` - are derived from this via ``_codebase_euler``. + Source of truth for the current rotation (quaternion-backed). + _rotvec : np.ndarray + Accumulated rotation vector (radians, scipy convention) - the + sum of all applied rotation vectors (a path integral), so full + turns are preserved, e.g. two full rotations around z read + (0, 0, 4 pi). For rotations around a single axis it represents + ``_R`` exactly; in general the exact orientation is ``_R``. + ``angx/angy/angz`` are derived from this. + _anchor_R : scipy.spatial.transform.Rotation + Reference orientation for per-segment rotation tracking (set + when an animation position is added or shown). + _anchor_rotvec : np.ndarray + Unwrapped rotation vector (radians, scipy convention) of the + rotation accumulated since ``_anchor_R``; used to encode + rotations beyond 180 degrees in animation segments. _last_mouse_x, _last_mouse_y : int Previous mouse position (Qt coords) during a trackball drag. viewport : tuple @@ -674,6 +767,9 @@ def __init__(self, window: QtWidgets.QMainWindow) -> None: super().__init__(window) self.window = window self._R = Rotation.identity() + self._rotvec = np.zeros(3) + self._anchor_R = Rotation.identity() + self._anchor_rotvec = np.zeros(3) self._pan_z = 0.0 self._last_mouse_x = 0 self._last_mouse_y = 0 @@ -692,37 +788,139 @@ def __init__(self, window: QtWidgets.QMainWindow) -> None: self.block_z = False self.setFocusPolicy(QtCore.Qt.FocusPolicy.ClickFocus) - # --- rotation angles --- # - def _codebase_euler(self) -> tuple[float, float, float]: - e = self._R.as_euler("XYZ") - return -float(e[0]), float(e[1]), float(e[2]) - + # --- rotation state --- # + # The codebase's angle convention (see ``render.rotation_matrix``) + # flips the sign of the x rotation relative to scipy's right-hand + # rule, i.e., rotation_matrix(a, 0, 0).as_rotvec() == (-a, 0, 0). + # The angx/angy/angz properties report the accumulated rotation + # vector in that convention. + # + # ``_rotvec`` accumulates the applied rotation vectors (a path + # integral) rather than re-deriving angles from ``_R``. An absolute + # orientation cannot count full turns - spinning a tilted structure + # by 360 degrees returns to the same orientation - so this is the + # only way the displayed angles can keep track of the number of + # turns (e.g. read 720 after two full rotations). For rotations + # around a single axis the accumulated values match the orientation + # exactly; for composed rotations around different axes they are a + # record of the applied rotations and the exact orientation is + # given by ``_R``. @property def angx(self) -> float: - return self._codebase_euler()[0] - - @angx.setter - def angx(self, value: float) -> None: - _, b, c = self._codebase_euler() - self._R = render.rotation_matrix(float(value), b, c) + return -float(self._rotvec[0]) @property def angy(self) -> float: - return self._codebase_euler()[1] - - @angy.setter - def angy(self, value: float) -> None: - a, _, c = self._codebase_euler() - self._R = render.rotation_matrix(a, float(value), c) + return float(self._rotvec[1]) @property def angz(self) -> float: - return self._codebase_euler()[2] + return float(self._rotvec[2]) - @angz.setter - def angz(self, value: float) -> None: - a, b, _ = self._codebase_euler() - self._R = render.rotation_matrix(a, b, float(value)) + @property + def rotation(self) -> Rotation: + """Current rotation of the localizations.""" + return self._R + + def set_rotation( + self, + R: Rotation, + rotvec: np.ndarray | None = None, + ) -> None: + """Set the absolute rotation and reset the per-segment anchor. + + Parameters + ---------- + R : scipy.spatial.transform.Rotation + New rotation. + rotvec : np.ndarray, optional + Accumulated rotation vector to restore as the displayed + angles (radians, scipy convention), e.g. saved earlier + from ``self._rotvec``; keeps full-turn information. If + None, the shortest rotation vector of ``R`` is used. + """ + self._R = R + if rotvec is None: + self._rotvec = R.as_rotvec() + else: + self._rotvec = np.asarray(rotvec, dtype=float).copy() + self.reset_rotation_anchor() + + def reset_rotation_anchor(self) -> None: + """Start tracking the rotation accumulated from the current + orientation (used for animation segments).""" + self._anchor_R = self._R + self._anchor_rotvec = np.zeros(3) + + def rotation_since_anchor(self) -> np.ndarray: + """Unwrapped rotation vector (radians, scipy convention) + accumulated since the last anchor; magnitude may exceed pi if + rotated beyond 180 degrees.""" + return self._anchor_rotvec.copy() + + def apply_rotation( + self, + rotvec: np.ndarray, + frame: str = "world", + ) -> None: + """Compose a rotation onto the current one and track turns. + + Parameters + ---------- + rotvec : np.ndarray + Rotation vector (radians, scipy convention). Its magnitude + may exceed pi; full turns are preserved in the displayed + angles and in the per-segment tracking. + frame : {"world", "object"}, optional + Frame in which ``rotvec`` is given. "world": the fixed + world/screen frame (rotation composed as ``delta * R``). + "object": the data's own (rotated) frame, i.e., rotation + around the data axes shown by the axes icon (composed as + ``R * delta``). Default is "world". + """ + rotvec = np.asarray(rotvec, dtype=float) + magnitude = float(np.linalg.norm(rotvec)) + if magnitude < 1e-12: + return + # accumulate the displayed angles (path integral; an absolute + # orientation cannot count full turns, so the applied rotation + # vectors are summed instead) + self._rotvec = self._rotvec + rotvec + # apply in sub-steps small enough for unambiguous unwrapping of + # the per-segment rotation + n_steps = max(1, int(np.ceil(magnitude / (np.pi / 2)))) + step = Rotation.from_rotvec(rotvec / n_steps) + for _ in range(n_steps): + if frame == "object": + self._R = self._R * step + else: + self._R = step * self._R + relative = self._R * self._anchor_R.inv() + self._anchor_rotvec = render.closest_rotvec( + relative, self._anchor_rotvec + ) + + def load_saved_rotation(self, info: dict) -> None: + """Restore the rotation saved by + ``RotationWindow.save_locs_rotated``. + + New files store a quaternion together with the accumulated + rotation angles (keeping full turns); legacy files store Euler + angles (angx, angy, angz, radians, codebase convention). + """ + angles = ( + info.get("angx", 0.0), + info.get("angy", 0.0), + info.get("angz", 0.0), + ) + if "Quaternion (x, y, z, w)" in info: + R = Rotation.from_quat(info["Quaternion (x, y, z, w)"]) + # angx is stored in the codebase convention (x negated + # relative to scipy's rotation vector) + rotvec = np.array([-angles[0], angles[1], angles[2]]) + self.set_rotation(R, rotvec=rotvec) + else: # legacy: Euler angles + self.set_rotation(render.rotation_matrix(*angles)) def sizeHint(self) -> QtCore.QSize: return QtCore.QSize(*self._size_hint) @@ -851,9 +1049,6 @@ def render_scene( viewport : tuple, optional Viewport to be rendered ``((y_min, x_min), (y_max, x_max))``. If None, takes current viewport. - ang : tuple, optional - Rotation angles to be rendered. If None, takes the current - angles. autoscale : bool, optional If True, optimally adjust contrast. use_cache : bool, optional @@ -882,7 +1077,7 @@ def render_scene( locs=locs, info=infos, **kwargs, - ang=(self.angx, self.angy, self.angz), + ang=self._R, contrast=contrast, invert_colors=self.window.dataset_dialog.wbackground.isChecked(), single_channel_colormap=cmap, @@ -1053,9 +1248,7 @@ def draw_rotation(self, image: QtGui.QImage) -> QtGui.QImage: Image with the drawn rotation axes icon. """ if self.window.rotation_action.isChecked(): - image = render.draw_rotation( - image=image, ang=(self.angx, self.angy, self.angz) - ) + image = render.draw_rotation(image=image, ang=self._R) return image def draw_rotation_angles(self, image: QtGui.QImage) -> QtGui.QImage: @@ -1098,40 +1291,28 @@ def draw_points(self, image: QtGui.QImage) -> QtGui.QImage: ) def rotation_input(self) -> None: - """Ask the user to input 3 rotation angles manually.""" - angx, ok = QtWidgets.QInputDialog.getDouble( - self, - "Rotation angle x", - "Angle x (degrees):", - 0, - decimals=2, - ) + """Ask the user to input 3 rotation angles manually. + + The rotations are applied sequentially around the data's own + (rotated) x, y and z axes - the axes shown by the axes icon - + so each displayed angle changes by exactly the entered amount + regardless of the current orientation. Angles beyond +/- 180 + degrees are preserved, e.g. 720 degrees encodes two full turns + (relevant for animations). + """ + angx, angy, angz, ok = RotateByAngleDialog.getParams(self) if ok: - angy, ok2 = QtWidgets.QInputDialog.getDouble( - self, - "Rotation angle y", - "Angle y (degrees):", - 0, - decimals=2, - ) - if ok2: - angz, ok3 = QtWidgets.QInputDialog.getDouble( - self, - "Rotation angle z", - "Angle z (degrees):", - 0, - decimals=2, - ) - if ok3: - self.angx += np.pi * angx / 180 - self.angy += np.pi * angy / 180 - self.angz += np.pi * angz / 180 - - self.update_scene() + # codebase convention: x angle is the negative of + # scipy's right-handed x rotation; "object" frame + # rotates around the data's own axes + self.apply_rotation([-np.radians(angx), 0.0, 0.0], frame="object") + self.apply_rotation([0.0, np.radians(angy), 0.0], frame="object") + self.apply_rotation([0.0, 0.0, np.radians(angz)], frame="object") + self.update_scene() def delete_rotation(self) -> None: """Reset rotation and any accumulated pan offset.""" - self._R = Rotation.identity() + self.set_rotation(Rotation.identity()) self._pan_z = 0.0 self.update_scene() @@ -1182,24 +1363,18 @@ def fit_in_view_rotated(self, get_viewport: bool = False) -> None: self.update_scene() def xy_projection(self) -> None: - """Set angles to 0 to get XY projection.""" - self.angx = 0 - self.angy = 0 - self.angz = 0 + """Reset rotation to get XY projection.""" + self.set_rotation(Rotation.identity()) self.update_scene() def xz_projection(self) -> None: - """Set angles to get XZ projection.""" - self.angx = np.pi / 2 - self.angy = 0 - self.angz = 0 + """Set rotation to get XZ projection.""" + self.set_rotation(render.rotation_matrix(np.pi / 2, 0, 0)) self.update_scene() def yz_projection(self) -> None: - """Set angles to get YZ projection.""" - self.angx = 0 - self.angy = np.pi / 2 - self.angz = 0 + """Set rotation to get YZ projection.""" + self.set_rotation(render.rotation_matrix(0, np.pi / 2, 0)) self.update_scene() def _arrow_pan(self, sx: float, sy: float) -> None: @@ -1332,12 +1507,13 @@ def _rotate_drag(self, event: QtGui.QMouseEvent) -> None: ax = 0.0 # Axis locks: pressing X/Y/Z constrains rotation to the - # corresponding *world* axis. Project the screen-frame rotation - # vector into world frame, zero the components we don't want, and - # rotate it back. + # corresponding *data* axis (the axes shown by the axes icon). + # Project the screen-frame rotation vector into the data frame, + # keep only the locked component and apply it around the data + # axis; the displayed angle then changes only for that axis. if self.block_x or self.block_y or self.block_z: - v_screen = np.array([ax, ay, az], dtype=float) - v_world = self._R.inv().apply(v_screen) + v_screen = render.rotation_matrix(ax, ay, az).as_rotvec() + v_object = self._R.inv().apply(v_screen) keep = np.array( [ 1.0 if self.block_x else 0.0, @@ -1345,20 +1521,10 @@ def _rotate_drag(self, event: QtGui.QMouseEvent) -> None: 1.0 if self.block_z else 0.0, ] ) - v_world = v_world * keep - v_screen = self._R.apply(v_world) - ax, ay, az = ( - float(v_screen[0]), - float(v_screen[1]), - float(v_screen[2]), - ) - - # The codebase's render.rotation_matrix flips the sign of the X - # rotation relative to scipy's right-hand-rule convention; using - # it here keeps screen-aligned tilts feeling identical to the old - # Euler-update behaviour at R = I. - delta_R = render.rotation_matrix(ax, ay, az) - self._R = delta_R * self._R + self.apply_rotation(v_object * keep, frame="object") + else: + delta_R = render.rotation_matrix(ax, ay, az) + self.apply_rotation(delta_R.as_rotvec()) self.update_scene() def _pan_drag(self, event: QtGui.QMouseEvent) -> None: @@ -1949,9 +2115,14 @@ def save_locs_rotated(self) -> None: "Pick": pick, "Pick shape": self.view_rot.pick_shape, "Pick size (nm)": size * pixelsize, - "angx": self.view_rot.angx, - "angy": self.view_rot.angy, - "angz": self.view_rot.angz, + # accumulated rotation angles incl. full turns + # (radians, codebase convention); legacy keys + "angx": float(self.view_rot.angx), + "angy": float(self.view_rot.angy), + "angz": float(self.view_rot.angz), + "Quaternion (x, y, z, w)": [ + float(q) for q in self.view_rot.rotation.as_quat() + ], } ] diff --git a/picasso/gui/spinna.py b/picasso/gui/spinna.py index 5cbf1d75..92646d2f 100644 --- a/picasso/gui/spinna.py +++ b/picasso/gui/spinna.py @@ -3893,7 +3893,9 @@ def update_prop_str_input_spins( def retrieve_results(): for i, ps in enumerate(prop_str): spin = self.prop_str_input_spins[i] + spin.blockSignals(True) spin.setValue(ps) + spin.blockSignals(False) button.released.connect(retrieve_results) self.prop_str_input.add_widget( diff --git a/picasso/lib.py b/picasso/lib.py index 4a1512cc..330ba9e6 100644 --- a/picasso/lib.py +++ b/picasso/lib.py @@ -1562,6 +1562,7 @@ def calculate_optimal_bins( bins : FloatArray1D Bins for display. """ + data = np.asarray(data) # positional indexing below; Series use labels n = len(data) if n == 0: return np.array([0.0, 1.0]) diff --git a/picasso/render.py b/picasso/render.py index 22d1bc6e..11146913 100755 --- a/picasso/render.py +++ b/picasso/render.py @@ -43,7 +43,7 @@ def render( Literal["gaussian", "gaussian_iso", "smooth", "convolve"] | None ) = None, min_blur_width: float = 0.0, - ang: tuple | None = None, + ang: tuple | Rotation | None = None, disp_px_size: float | None = None, ) -> tuple[int, lib.FloatArray2D]: """Render localizations given FOV and blur method. @@ -74,9 +74,11 @@ def render( localizations which is the median localization precision. min_blur_width : float, optional Minimum size of blur (camera pixels). - ang : tuple, optional - Rotation angles of locs around x, y and z axes in radians. If - None, locs are not rotated. + ang : tuple or scipy.spatial.transform.Rotation, optional + Rotation of locs; either a scipy Rotation (e.g. built from a + quaternion) or a tuple of 3 rotation angles around the x, y + and z axes in radians (legacy Euler convention, see + ``rotation_matrix``). If None, locs are not rotated. disp_px_size : float, optional Display pixel size in nm. Will replace oversampling in v0.11.0. @@ -637,7 +639,7 @@ def _fill_gaussian_rot( sz: lib.FloatArray1D, n_pixel_x: int, n_pixel_y: int, - ang: tuple[float, float, float], + rot_matrix: lib.Array3x3, ) -> None: """Fill image with rotated gaussian-blurred localizations. @@ -654,38 +656,12 @@ def _fill_gaussian_rot( Localization precision in x, y and z for each localization. n_pixel_x, n_pixel_y : int Number of pixels in x and y. - ang : tuple - Rotation angles of locs around x, y and z axes (radians). + rot_matrix : lib.Array3x3 + Rotation matrix (float32) applied to the localizations. """ n_locs = len(x) if n_locs == 0: return - angx, angy, angz = ang - rot_mat_x = np.array( - [ - [1.0, 0.0, 0.0], - [0.0, np.cos(angx), np.sin(angx)], - [0.0, -np.sin(angx), np.cos(angx)], - ], - dtype=np.float32, - ) - rot_mat_y = np.array( - [ - [np.cos(angy), 0.0, np.sin(angy)], - [0.0, 1.0, 0.0], - [-np.sin(angy), 0.0, np.cos(angy)], - ], - dtype=np.float32, - ) - rot_mat_z = np.array( - [ - [np.cos(angz), -np.sin(angz), 0.0], - [np.sin(angz), np.cos(angz), 0.0], - [0.0, 0.0, 1.0], - ], - dtype=np.float32, - ) - rot_matrix = (rot_mat_x @ rot_mat_y @ rot_mat_z).astype(np.float32) rot_matrixT = np.ascontiguousarray(rot_matrix.T) for i in range(n_locs): @@ -804,7 +780,7 @@ def render_hist( x_min: float, y_max: float, x_max: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Alias for _render_hist which will be a private function in v0.11.0. Kept for backward compatibility but will be removed in @@ -826,7 +802,7 @@ def _render_hist( x_min: float, y_max: float, x_max: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Render localizations with no blur by assigning them to pixels. @@ -840,9 +816,11 @@ def _render_hist( Minimum y and x coordinates to be rendered (camera pixels) y_max, x_max : float Maximum y and x coordinates to be rendered (camera pixels) - ang : tuple (default=None) - Rotation angles of locs around x, y and z axes in radians. If - None, locs are not rotated. + ang : tuple or scipy.spatial.transform.Rotation, optional + Rotation of locs; either a scipy Rotation (e.g. built from a + quaternion) or a tuple of 3 rotation angles around the x, y + and z axes in radians (legacy Euler convention, see + ``rotation_matrix``). If None, locs are not rotated. Returns ------- @@ -1017,7 +995,7 @@ def render_gaussian( y_max: float, x_max: float, min_blur_width: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Alias for _render_gaussian which will be a private function in v0.11.0. Kept for backward compatibility but will be removed in v0.11.0. Use @@ -1047,7 +1025,7 @@ def _render_gaussian( y_max: float, x_max: float, min_blur_width: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Render localizations with with individual localization precision which differs in x and y. @@ -1064,9 +1042,11 @@ def _render_gaussian( Minimum and maximum x coordinates to be rendered (camera pixels). min_blur_width : float Minimum localization precision (camera pixels). - ang : tuple, optional - Rotation angles of localizations around x, y and z axes in - radians. If None, localizations are not rotated. + ang : tuple or scipy.spatial.transform.Rotation, optional + Rotation of localizations; either a scipy Rotation (e.g. built + from a quaternion) or a tuple of 3 rotation angles around the + x, y and z axes in radians (legacy Euler convention, see + ``rotation_matrix``). If None, localizations are not rotated. Returns ------- @@ -1124,7 +1104,12 @@ def _render_gaussian( sx = blur_width[in_view] sz = blur_depth[in_view] - _fill_gaussian_rot(image, x, y, sx, sy, sz, n_pixel_x, n_pixel_y, ang) + rot_matrix = np.ascontiguousarray( + to_rotation(ang).as_matrix(), dtype=np.float32 + ) + _fill_gaussian_rot( + image, x, y, sx, sy, sz, n_pixel_x, n_pixel_y, rot_matrix + ) n = len(x) return n, image @@ -1138,7 +1123,7 @@ def render_gaussian_iso( y_max: float, x_max: float, min_blur_width: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Alias for _render_gaussian_iso which will be a private function in v0.11.0. Kept for backward compatibility but will be removed in v0.11.0. Use @@ -1168,7 +1153,7 @@ def _render_gaussian_iso( y_max: float, x_max: float, min_blur_width: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Same as ``_render_gaussian``, but uses the same localization precision in x and y.""" @@ -1221,7 +1206,12 @@ def _render_gaussian_iso( sx = sy sz = blur_depth[in_view] - _fill_gaussian_rot(image, x, y, sx, sy, sz, n_pixel_x, n_pixel_y, ang) + rot_matrix = np.ascontiguousarray( + to_rotation(ang).as_matrix(), dtype=np.float32 + ) + _fill_gaussian_rot( + image, x, y, sx, sy, sz, n_pixel_x, n_pixel_y, rot_matrix + ) return len(x), image @@ -1234,7 +1224,7 @@ def render_convolve( y_max: float, x_max: float, min_blur_width: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Alias for _render_convolve which will be a private function in v0.11.0. Kept for backward compatibility but will be removed in v0.11.0. Use @@ -1264,7 +1254,7 @@ def _render_convolve( y_max: float, x_max: float, min_blur_width: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Render localizations with with global localization precision, i.e. each localization is blurred by the median localization @@ -1282,9 +1272,11 @@ def _render_convolve( Maximum y and x coordinates to be rendered (camera pixels). min_blur_width : float Minimum localization precision (camera pixels). - ang : tuple, optional - Rotation angles of localizations around x, y and z axes in - radians. If None, localizations are not rotated. + ang : tuple or scipy.spatial.transform.Rotation, optional + Rotation of localizations; either a scipy Rotation (e.g. built + from a quaternion) or a tuple of 3 rotation angles around the + x, y and z axes in radians (legacy Euler convention, see + ``rotation_matrix``). If None, localizations are not rotated. Returns ------- @@ -1334,7 +1326,7 @@ def render_smooth( x_min: float, y_max: float, x_max: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Alias for _render_smooth which will be a private function in v0.11.0. Kept for backward compatibility but will be removed in v0.11.0. Use _render_smooth @@ -1361,7 +1353,7 @@ def _render_smooth( x_min: float, y_max: float, x_max: float, - ang: tuple[float, float, float] | None = None, + ang: tuple[float, float, float] | Rotation | None = None, ) -> tuple[int, lib.FloatArray2D]: """Render localizations with with blur of one display pixel (set by oversampling). @@ -1376,9 +1368,11 @@ def _render_smooth( Minimum y and x coordinates to be rendered (camera pixels). y_max, x_max : float Maximum y and x coordinates to be rendered (camera pixels). - ang : tuple, optional - Rotation angles of localizations around x, y and z axes in - radians. If None, localizations are not rotated. + ang : tuple or scipy.spatial.transform.Rotation, optional + Rotation of localizations; either a scipy Rotation (e.g. built + from a quaternion) or a tuple of 3 rotation angles around the + x, y and z axes in radians (legacy Euler convention, see + ``rotation_matrix``). If None, localizations are not rotated. Returns ------- @@ -1504,6 +1498,76 @@ def rotation_matrix(angx: float, angy: float, angz: float) -> Rotation: return Rotation.from_matrix(rot_mat) +def to_rotation( + ang: tuple[float, float, float] | Rotation | None, +) -> Rotation | None: + """Normalize a rotation input to a scipy Rotation. + + Parameters + ---------- + ang : tuple, Rotation or None + Rotation to be normalized. Either a scipy Rotation (e.g. built + from a quaternion; used as is), a tuple of 3 rotation angles + around the x, y and z axes in radians (legacy Euler convention, + see ``rotation_matrix``), or None. + + Returns + ------- + scipy.spatial.transform.Rotation or None + The rotation as a scipy Rotation, or None if ``ang`` is None. + """ + if ang is None: + return None + if isinstance(ang, Rotation): + return ang + return rotation_matrix(*ang) + + +def closest_rotvec( + rotation: Rotation, + reference: lib.FloatArray1D, +) -> lib.FloatArray1D: + """Find the rotation vector representing ``rotation`` that is + closest to ``reference``. + + A rotation vector (axis * angle, radians) is not unique: adding + full turns (2 pi) around the same axis yields the same rotation. + This function picks the representation closest to ``reference``, + which allows keeping track of rotations beyond +/- 180 degrees + (e.g. unwrapping a continuously updated rotation, or encoding + multiple full turns in an animation segment). + + Parameters + ---------- + rotation : scipy.spatial.transform.Rotation + The rotation to be represented. + reference : lib.FloatArray1D + Rotation vector (radians) to which the result should be + closest. + + Returns + ------- + rotvec : lib.FloatArray1D + Rotation vector (radians) such that + ``Rotation.from_rotvec(rotvec) == rotation``, with the number + of full turns chosen to match ``reference``. Its magnitude may + exceed pi. + """ + reference = np.asarray(reference, dtype=float) + base = rotation.as_rotvec() # magnitude <= pi + theta = np.linalg.norm(base) + if theta < 1e-9: + # identity rotation; keep the full turns along reference's axis + ref_magnitude = np.linalg.norm(reference) + if ref_magnitude < 1e-9: + return np.zeros(3) + n_turns = np.round(ref_magnitude / (2 * np.pi)) + return reference * (2 * np.pi * n_turns / ref_magnitude) + axis = base / theta + n_turns = np.round((axis @ reference - theta) / (2 * np.pi)) + return axis * (theta + 2 * np.pi * n_turns) + + def locs_rotation( locs: pd.DataFrame, oversampling: float, @@ -1511,7 +1575,7 @@ def locs_rotation( x_max: float, y_min: float, y_max: float, - ang: tuple[float, float, float], + ang: tuple[float, float, float] | Rotation, ) -> tuple[ lib.FloatArray1D, lib.FloatArray1D, lib.BoolArray1D, lib.FloatArray1D ]: @@ -1527,9 +1591,10 @@ def locs_rotation( Minimum y and x coordinate to be rendered (camera pixels). y_max, x_max : float Maximum y and x coordinate to be rendered (camera pixels). - ang : tuple - Rotation angles of localizations around x, y and z axes in - radians. + ang : tuple or scipy.spatial.transform.Rotation + Rotation of localizations; either a scipy Rotation or a tuple + of 3 rotation angles around the x, y and z axes in radians + (legacy Euler convention, see ``rotation_matrix``). Returns ------- @@ -1551,7 +1616,7 @@ def locs_rotation( locs_coord[:, 1] -= y_min + (y_max - y_min) / 2 # rotate locs - R = rotation_matrix(ang[0], ang[1], ang[2]) + R = to_rotation(ang) locs_coord = R.apply(locs_coord) # unshift locs @@ -2538,7 +2603,7 @@ def draw_minimap( def draw_rotation( image: QtGui.QImage, - ang: tuple[float, float, float], + ang: tuple[float, float, float] | Rotation, axis_length: int = 30, axis_center: tuple[int, int] = (50, -50), # bottom left ) -> QtGui.QImage: @@ -2548,8 +2613,10 @@ def draw_rotation( ---------- image : QImage Image containing rendered localizations. - ang : tuple of float - Rotation angles around x, y, and z axes in radians. + ang : tuple of float or scipy.spatial.transform.Rotation + Rotation of the localizations; either a scipy Rotation or a + tuple of 3 rotation angles around the x, y and z axes in + radians (legacy Euler convention, see ``rotation_matrix``). axis_length : int, optional Length of the rotation axes in display pixels. Default is 30. axis_center : tuple of int, optional @@ -2592,7 +2659,7 @@ def draw_rotation( # rotate these points coordinates = [[xx, xy, xz], [yx, yy, yz], [zx, zy, zz]] - R = rotation_matrix(*ang) + R = to_rotation(ang) coordinates = R.apply(coordinates).astype(int) (xx, xy, xz) = coordinates[0] (yx, yy, yz) = coordinates[1] @@ -2635,7 +2702,8 @@ def draw_rotation_angles( image : QImage Image containing rendered localizations. ang : tuple of float - Rotation angles around x, y, and z axes in radians. + Rotation angles (or rotation-vector components) around x, y, + and z axes in radians, used only for the displayed text. color : QColor, optional Color of the text. Default is white. @@ -2667,7 +2735,7 @@ def render_scene( Literal["gaussian", "gaussian_iso", "smooth", "convolve"] | None ) = None, min_blur_width: float = 0.0, - ang: tuple | None = None, + ang: tuple | Rotation | None = None, contrast: tuple[float, float] | None = None, invert_colors: bool = False, single_channel_colormap: str | lib.FloatArray2D = "magma", @@ -2738,9 +2806,11 @@ def render_scene( localizations which is the median localization precision. min_blur_width : float, optional Minimum size of blur (camera pixels). - ang : tuple, optional - Rotation angles of locs around x, y and z axes in radians. If - None, locs are not rotated. + ang : tuple or scipy.spatial.transform.Rotation, optional + Rotation of locs; either a scipy Rotation (e.g. built from a + quaternion) or a tuple of 3 rotation angles around the x, y + and z axes in radians (legacy Euler convention, see + ``rotation_matrix``). If None, locs are not rotated. contrast : tuple of float, optional Contrast limits for scaling. If None, contrast is automatically determined. @@ -2751,10 +2821,17 @@ def render_scene( corresponding pyplot colormap is selected. If a 2D array, a 256x4 array is expected with values between 0 and 1. Default is 'magma'. - colors : list of tuples, optional - List of RGB tuples corresponding to the colors of the channels. - Only needs to be specified for multi-channel data. Must range - between 0 and 1. Default is None. + colors : list of tuples or list of lib.FloatArray2D, optional + Colors of the channels, one entry per channel. Each entry is + either an ``(r, g, b)`` tuple with values between 0 and 1 + (the channel is rendered as ``intensity * rgb``) or a + ``(256, 3)`` lookup table with values between 0 and 1 (the + channel intensity is indexed into the LUT, allowing + per-channel colormaps; see ``solid_to_lut`` and + ``stops_to_lut``). The two forms cannot be mixed. Channels are + blended additively. Only needs to be specified for + multi-channel data. Default is None, in which case colors are + taken from ``lib.get_colors``. relative_intensities : list of float, optional List of relative intensities for each channel. Only needs to be specified for multi-channel data. Default is None, in which @@ -2854,7 +2931,7 @@ def _render_multi_channel( Literal["gaussian", "gaussian_iso", "smooth", "convolve"] | None ) = None, min_blur_width: float = 0.0, - ang: tuple | None = None, + ang: tuple | Rotation | None = None, contrast: tuple[float, float] | None = None, relative_intensities: list[float] | None = None, invert_colors: bool = False, @@ -2933,7 +3010,7 @@ def _render_single_channel( Literal["gaussian", "gaussian_iso", "smooth", "convolve"] | None ) = None, min_blur_width: float = 0.0, - ang: tuple | None = None, + ang: tuple | Rotation | None = None, contrast: tuple[float, float] | None = None, invert_colors: bool = False, single_channel_colormap: str = "magma", @@ -3244,41 +3321,91 @@ def optimal_scalebar_length(pixelsize: int | float, width: int | float) -> int: return scalebar +def _normalize_animation_positions( + positions: list, +) -> list[tuple[Rotation, tuple]]: + """Normalize animation checkpoints to (Rotation, viewport) tuples. + + Accepts two formats: each position is (rotation, viewport) with + rotation being a scipy Rotation, as well as the legacy format + (angle_x, angle_y, angle_z, viewport) with Euler angles in radians + (see ``rotation_matrix``), which is deprecated. + """ + normalized = [] + legacy = False + for p in positions: + if len(p) == 2 and isinstance(p[0], Rotation): + normalized.append((p[0], p[1])) + elif len(p) == 4: + legacy = True + normalized.append((rotation_matrix(p[0], p[1], p[2]), p[3])) + else: + raise ValueError( + "Each position must be a tuple (rotation, viewport) with " + "rotation being a scipy Rotation, or the deprecated " + "4-element form (angle_x, angle_y, angle_z, viewport)." + ) + if legacy: + lib.deprecation_warning( + "Deprecation warning: passing animation positions as Euler " + "angles (angle_x, angle_y, angle_z, viewport) is deprecated " + "and will be removed in v0.12.0. Pass (rotation, viewport) " + "with rotation being a scipy.spatial.transform.Rotation " + "instead." + ) + return normalized + + def _animation_sequence( - positions: list[list[float, float, float, tuple]], + positions: list[tuple[Rotation, tuple]], durations: list[float], fps: int, -) -> tuple[list, list]: - """Calculate the sequence of angles and viewports for the animation. - See ``build_animation`` for more details.""" - n_frames = [0] - for i in range(len(positions) - 1): - n_frames.append(int(fps * durations[i])) - - # find rotation angles and viewport for each frame - angles = [] + segment_rotations: list | None = None, +) -> tuple[list[Rotation], list]: + """Calculate the sequence of rotations and viewports for the + animation. See ``build_animation`` for more details. + + Each segment is interpolated along the geodesic between the two + checkpoint rotations at constant angular velocity (slerp). If + ``segment_rotations`` is given, the corresponding rotation vector + defines the rotation path of each segment and may include full + turns (magnitude beyond pi), e.g. 4 pi for two full turns. The + vector is snapped to the true relative rotation between the two + checkpoints (see ``closest_rotvec``) so that each segment always + ends exactly at the next checkpoint. This is equivalent to + splitting a segment with a rotation larger than 180 degrees into + sub-180-degree pieces and applying slerp to each piece.""" + rotations = [] viewports = [] for i in range(len(positions) - 1): - # angles - x1, y1, z1 = positions[i][:3] - x2, y2, z2 = positions[i + 1][:3] - current_angles = np.linspace( - [x1, y1, z1], [x2, y2, z2], n_frames[i + 1] + n_frames = int(fps * durations[i]) + + # rotations + R1, vp1 = positions[i] + R2, vp2 = positions[i + 1] + relative = R2 * R1.inv() + if segment_rotations is not None: + rotvec = closest_rotvec( + relative, np.asarray(segment_rotations[i], dtype=float) + ) + else: + rotvec = relative.as_rotvec() + fractions = np.linspace(0, 1, n_frames) + rotations.extend( + Rotation.from_rotvec(fraction * rotvec) * R1 + for fraction in fractions ) - angles.extend(current_angles) # viewports - vp1 = positions[i][3] - vp2 = positions[i + 1][3] - ymin = np.linspace(vp1[0][0], vp2[0][0], n_frames[i + 1]) - xmin = np.linspace(vp1[0][1], vp2[0][1], n_frames[i + 1]) - ymax = np.linspace(vp1[1][0], vp2[1][0], n_frames[i + 1]) - xmax = np.linspace(vp1[1][1], vp2[1][1], n_frames[i + 1]) + ymin = np.linspace(vp1[0][0], vp2[0][0], n_frames) + xmin = np.linspace(vp1[0][1], vp2[0][1], n_frames) + ymax = np.linspace(vp1[1][0], vp2[1][0], n_frames) + xmax = np.linspace(vp1[1][1], vp2[1][1], n_frames) current_viewports = [ ((ymin[j], xmin[j]), (ymax[j], xmax[j])) for j in range(len(ymin)) ] viewports.extend(current_viewports) - return angles, viewports + return rotations, viewports def build_animation( @@ -3286,10 +3413,14 @@ def build_animation( locs: pd.DataFrame | list[pd.DataFrame], info: list[dict] | list[list[dict]], *, - positions: list[tuple[float, float, float, tuple]], + positions: ( + list[tuple[Rotation, tuple]] + | list[tuple[tuple[float, float, float], tuple]] + ), durations: list[float], disp_px_size: int | float, # nm image_size: tuple[int, int], + segment_rotations: list | None = None, blur_method: ( Literal["gaussian", "gaussian_iso", "smooth", "convolve"] | None ) = None, @@ -3306,7 +3437,7 @@ def build_animation( ) = None, ) -> None: """Build an animation of rendered localizations given the - checkpoints (angle, viewport, etc) and the time between them. + checkpoints (rotation, viewport, etc) and the time between them. Parameters ---------- @@ -3328,13 +3459,27 @@ def build_animation( image_size : tuple of int Size of the rendered image in pixels, given as (width, height). positions : list - Each element determines the checkpoint of the animation, which - is a tuple of 4 elements: (angle_x, angle_y, angle_z, viewport). - Angles are in radians. Viewport is given as ((y_min, x_min), - (y_max, x_max)) in camera pixels. + Each element determines a checkpoint of the animation, which + is a tuple of 2 elements: (rotation, viewport). Rotation is a + ``scipy.spatial.transform.Rotation`` defining the orientation + of the localizations at the checkpoint. Viewport is given as + ((y_min, x_min), (y_max, x_max)) in camera pixels. The + deprecated legacy format (angle_x, angle_y, angle_z, viewport) + with Euler angles in radians (see ``rotation_matrix``) is also + accepted and will be removed in v0.12.0. durations : list List of durations in seconds between the checkpoints. Must have the same length as positions - 1. + segment_rotations : list, optional + One rotation vector (3 floats, radians, scipy convention) per + segment, i.e., of length len(positions) - 1, describing the + full rotation path from one checkpoint to the next. The + magnitude may exceed pi to encode rotations larger than 180 + degrees, e.g. (0, 0, 4 * pi) for two full turns around the z + axis. Each vector is snapped to the true relative rotation + between its two checkpoints, so the segment always ends + exactly at the next checkpoint. If None, each segment follows + the shortest path (slerp). Default is None. blur_method : {"gaussian", "gaussian_iso", "smooth", "convolve"} or None, \ optional Defines localizations' blur. The string has to be one of @@ -3359,10 +3504,17 @@ def build_animation( corresponding pyplot colormap is selected. If a 2D array, a 256x4 array is expected with values between 0 and 1. Default is 'magma'. - colors : list of tuples, optional - List of RGB tuples corresponding to the colors of the channels. - Only needs to be specified for multi-channel data. Must range - between 0 and 1. Default is None. + colors : list of tuples or list of lib.FloatArray2D, optional + Colors of the channels, one entry per channel. Each entry is + either an ``(r, g, b)`` tuple with values between 0 and 1 + (the channel is rendered as ``intensity * rgb``) or a + ``(256, 3)`` lookup table with values between 0 and 1 (the + channel intensity is indexed into the LUT, allowing + per-channel colormaps; see ``solid_to_lut`` and + ``stops_to_lut``). The two forms cannot be mixed. Channels are + blended additively. Only needs to be specified for + multi-channel data. Default is None, in which case colors are + taken from ``lib.get_colors``. relative_intensities : list of float, optional List of relative intensities for each channel. Only needs to be specified for multi-channel data. Default is None, in which @@ -3397,10 +3549,16 @@ def build_animation( assert ( isinstance(positions, list) and len(positions) >= 2 ), "positions must be a list with at least 2 elements." - assert all(len(p) == 4 for p in positions), ( - "Each position must be a tuple/list of 4 elements: " - "(angle_x, angle_y, angle_z, viewport)." - ) + positions = _normalize_animation_positions(positions) + if segment_rotations is not None: + assert ( + isinstance(segment_rotations, list) + and len(segment_rotations) == len(positions) - 1 + ), "segment_rotations must be a list of length len(positions) - 1." + assert all( + np.asarray(rotvec, dtype=float).shape == (3,) + for rotvec in segment_rotations + ), "Each segment rotation must be a rotation vector of 3 floats." assert ( isinstance(durations, list) and len(durations) == len(positions) - 1 ), "durations must be a list of length len(positions) - 1." @@ -3448,9 +3606,6 @@ def build_animation( assert ( len(colors) == n_channels ), "colors must have one entry per channel." - assert all( - len(c) == 3 and all(0.0 <= v <= 1.0 for v in c) for c in colors - ), "Each color must be an RGB tuple with values between 0 and 1." if relative_intensities is not None: n_channels = len(locs) if isinstance(locs, list) else 1 assert ( @@ -3475,6 +3630,7 @@ def build_animation( info=info, positions=positions, durations=durations, + segment_rotations=segment_rotations, disp_px_size=disp_px_size, image_size=image_size, blur_method=blur_method, @@ -3494,8 +3650,9 @@ def _build_animation( path: str, locs: pd.DataFrame | list[pd.DataFrame], info: list[dict] | list[list[dict]], - positions: list[tuple[float, float, float, tuple]], + positions: list[tuple[Rotation, tuple]], durations: list[float], + segment_rotations: list | None, disp_px_size: int | float, image_size: tuple[int, int], blur_method: ( @@ -3513,7 +3670,9 @@ def _build_animation( ) -> None: """Internal function to build an animation of rendered localizations given the checkpoints. See ``build_animation`` for more details.""" - angles, viewports = _animation_sequence(positions, durations, fps) + rotations, viewports = _animation_sequence( + positions, durations, fps, segment_rotations=segment_rotations + ) # width and height for building the animation; must be divisible by 16 # as ffmpeg codecs require this for proper encoding @@ -3526,10 +3685,10 @@ def _build_animation( use_tqdm = progress_callback == "console" if use_tqdm: iter_range = tqdm( - range(len(angles)), desc="Building animation", unit="frame" + range(len(rotations)), desc="Building animation", unit="frame" ) else: - iter_range = range(len(angles)) + iter_range = range(len(rotations)) for i in iter_range: if callable(progress_callback): @@ -3546,7 +3705,7 @@ def _build_animation( info=info, disp_px_size=disp_px_size_, viewport=viewports[i], - ang=angles[i], + ang=rotations[i], blur_method=blur_method, min_blur_width=min_blur_width, contrast=contrast_, @@ -3570,30 +3729,32 @@ def _build_animation( video_writer.append_data(frame) if callable(progress_callback): - progress_callback(len(angles)) + progress_callback(len(rotations)) video_writer.close() # save a yaml with animation settings, note that yaml does not support # numpy types and arrays - angles_yaml = [ - ( - float(np.degrees(p[0])), - float(np.degrees(p[1])), - float(np.degrees(p[2])), - ) - for p in positions - ] + quaternions_yaml = [[float(q) for q in R.as_quat()] for R, _ in positions] viewports_yaml = [ ( - (float(p[3][0][0]), float(p[3][0][1])), - (float(p[3][1][0]), float(p[3][1][1])), + (float(vp[0][0]), float(vp[0][1])), + (float(vp[1][0]), float(vp[1][1])), ) - for p in positions + for _, vp in positions + ] + if segment_rotations is None: + segment_rotations = [ + (positions[i + 1][0] * positions[i][0].inv()).as_rotvec() + for i in range(len(positions) - 1) + ] + segments_yaml = [ + [float(np.degrees(v)) for v in rotvec] for rotvec in segment_rotations ] anim_settings = { "Generated by": f"Picasso v{__version__} Render 3D Animation", "FPS": fps, - "Angles at checkpoints (x, y, z) (deg)": angles_yaml, + "Quaternions at checkpoints (x, y, z, w)": quaternions_yaml, + "Rotations between checkpoints (x, y, z) (deg)": segments_yaml, "Viewports at checkpoints (camera pixels)": viewports_yaml, "Durations (s)": durations, } diff --git a/picasso/version.py b/picasso/version.py index 1f4c4d43..17c1a626 100644 --- a/picasso/version.py +++ b/picasso/version.py @@ -1 +1 @@ -__version__ = "0.10.1" +__version__ = "0.10.2" diff --git a/picasso/zfit.py b/picasso/zfit.py index 0dd6a8b8..27576db7 100644 --- a/picasso/zfit.py +++ b/picasso/zfit.py @@ -49,6 +49,7 @@ def calibrate_z( d: float, magnification_factor: float, path: str | None = None, + frame_bounds: tuple[int, int] | None = None, ) -> dict: """Given localizations of a calibration sample (e.g., gold beads at different z positions), calibrate the z-axis by fitting a polynomial @@ -72,6 +73,12 @@ def calibrate_z( path : str, optional Path to save the calibration data as a YAML file. If None, the calibration data will not be saved. Default is None. + frame_bounds : tuple, optional + Minimum and maximum frame numbers to consider for the + calibration. If None, all frames are used. If only min or max + is to be specified, the other is to be set to None, for example, + ``(5, None)`` sets minimum frame to 5 without maximum frame. + Default is None. Returns ------- @@ -85,6 +92,16 @@ def calibrate_z( frame_range = np.arange(n_frames) z_range = -(frame_range * d - range / 2) # negative so that the # first frames of a bottom-to-up scan are positive z coordinates. + if frame_bounds is not None: + frame_min, frame_max = frame_bounds + frame_min = frame_min or 0 + frame_max = frame_max or (n_frames - 1) + # frame bounds are inclusive, like in picasso.localize + frame_range = frame_range[frame_min : frame_max + 1] + z_range = z_range[frame_min : frame_max + 1] + locs = locs[ + (locs["frame"] >= frame_min) & (locs["frame"] <= frame_max) + ] mean_sx = np.array( [locs["sx"][locs["frame"] == _].mean() for _ in frame_range] @@ -99,8 +116,11 @@ def calibrate_z( [locs["sy"][locs["frame"] == _].var() for _ in frame_range] ) - keep_x = (locs["sx"] - mean_sx[locs["frame"]]) ** 2 < var_sx[locs["frame"]] - keep_y = (locs["sy"] - mean_sy[locs["frame"]]) ** 2 < var_sy[locs["frame"]] + # index of each localization's frame in the (possibly bounded) + # frame_range + frame_idx = locs["frame"] - frame_range[0] + keep_x = (locs["sx"] - mean_sx[frame_idx]) ** 2 < var_sx[frame_idx] + keep_y = (locs["sy"] - mean_sy[frame_idx]) ** 2 < var_sy[frame_idx] keep = keep_x & keep_y locs = locs[keep] @@ -134,6 +154,7 @@ def calibrate_z( "Step size in nm": float(d), "Magnification factor": float(magnification_factor), "Path": path if path is not None else "N/A", + "Frame bounds": frame_bounds, } if path is not None: with open(path, "w") as f: @@ -142,6 +163,7 @@ def calibrate_z( # pixelsize does not matter here anyway locs = _fit_z(locs, info, calibration, magnification_factor, pixelsize=130) locs["z"] /= magnification_factor + frame_idx = locs["frame"] - frame_range[0] plt.figure(figsize=(18, 10)) @@ -186,7 +208,7 @@ def calibrate_z( plt.legend(loc="best") ax = plt.subplot(234) - plt.plot(z_range[locs["frame"]], locs["z"], ".k", alpha=0.1) + plt.plot(z_range[frame_idx], locs["z"], ".k", alpha=0.1) plt.plot( [z_range.min(), z_range.max()], [z_range.min(), z_range.max()], @@ -201,7 +223,7 @@ def calibrate_z( plt.legend(loc="best") ax = plt.subplot(235) - deviation = locs["z"] - z_range[locs["frame"]] + deviation = locs["z"] - z_range[frame_idx] bins = lib.calculate_optimal_bins(deviation, max_n_bins=1000) plt.hist(deviation, bins) plt.xlabel("Deviation to true position") diff --git a/pyproject.toml b/pyproject.toml index 1e033b1b..4be6e658 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,6 +54,15 @@ dependencies = [ ] [project.optional-dependencies] +cuda11 = [ + "numba-cuda[cu11]" +] +cuda12 = [ + "numba-cuda[cu12]" +] +cuda13 = [ + "numba-cuda[cu13]" +] dev = [ "build", "black", @@ -86,6 +95,33 @@ installer = [ "PyImarisWriter==0.7.0; sys_platform=='win32'", "hdf5plugin==6.0.0; sys_platform=='win32'" ] +installer_gpu = [ + "PyQt6==6.11.0", + "numpy==2.4.4", + "matplotlib==3.10.8", + "h5py==3.16.0", + "numba==0.65.0", + "scipy==1.17.1", + "pandas==2.3.3", + "tables==3.11.1", + "pyyaml==6.0.3", + "scikit-learn==1.7.2", + "tqdm==4.67.3", + "streamlit==1.56.0", + "nd2==0.11.3", + "sqlalchemy==2.0.49", + "plotly-express==0.4.1", + "watchdog==6.0.0", + "psutil==7.2.2", + "playsound3==3.3.1", + "imageio==2.37.3", + "imageio-ffmpeg==0.6.0", + "tifffile==2026.4.11", + "pyinstaller==6.19.0", + "PyImarisWriter==0.7.0; sys_platform=='win32'", + "hdf5plugin==6.0.0; sys_platform=='win32'", + "numba-cuda[cu12]==0.30.2" +] [project.urls] Homepage = "https://github.com/jungmannlab/picasso" diff --git a/tests/test_render.py b/tests/test_render.py index 6669c40f..b342df12 100644 --- a/tests/test_render.py +++ b/tests/test_render.py @@ -11,6 +11,7 @@ import pandas as pd import pytest from PyQt6 import QtCore, QtGui +from scipy.spatial.transform import Rotation from picasso import io, masking, render @@ -540,6 +541,99 @@ def test_locs_rotation_in_view_consistency(self, locs_3d): assert len(y_out) == n_in_view assert len(z_out) == n_in_view + def test_to_rotation_none(self): + assert render.to_rotation(None) is None + + def test_to_rotation_passes_rotation_through(self): + R = Rotation.from_rotvec([0.1, -0.2, 0.3]) + assert render.to_rotation(R) is R + + def test_to_rotation_legacy_euler_equivalence(self): + ang = (0.4, -0.7, 1.1) + R = render.to_rotation(ang) + expected = render.rotation_matrix(*ang) + assert np.allclose(R.as_matrix(), expected.as_matrix()) + + def test_locs_rotation_accepts_rotation_object(self, locs_3d): + ang = (0.1, 0.2, 0.3) + out_tuple = render.locs_rotation( + locs_3d, 5.0, 0.0, 32.0, 0.0, 32.0, ang + ) + out_rotation = render.locs_rotation( + locs_3d, 5.0, 0.0, 32.0, 0.0, 32.0, render.rotation_matrix(*ang) + ) + for a, b in zip(out_tuple, out_rotation): + assert np.allclose(a, b) + + def test_render_accepts_rotation_object(self, locs_3d, info): + """Rendering with a scipy Rotation matches the legacy Euler + tuple input for all rotation-aware code paths.""" + ang = (0.5, 0.3, 0.2) + for blur_method in (None, "gaussian"): + _, im_tuple = render.render( + locs_3d, + info, + oversampling=5, + viewport=FULL_VIEWPORT, + blur_method=blur_method, + ang=ang, + ) + _, im_rotation = render.render( + locs_3d, + info, + oversampling=5, + viewport=FULL_VIEWPORT, + blur_method=blur_method, + ang=render.rotation_matrix(*ang), + ) + assert np.allclose(im_tuple, im_rotation) + + +class TestClosestRotvec: + def test_zero_reference_returns_base(self): + R = Rotation.from_rotvec([0.3, 0.0, 0.0]) + out = render.closest_rotvec(R, np.zeros(3)) + assert np.allclose(out, [0.3, 0.0, 0.0]) + + def test_unwraps_full_turns(self): + """A 10-degree rotation with a reference near 370 degrees must + unwrap to 370 degrees.""" + axis = np.array([0.0, 0.0, 1.0]) + R = Rotation.from_rotvec(np.radians(10) * axis) + reference = np.radians(365) * axis + out = render.closest_rotvec(R, reference) + assert np.allclose(out, np.radians(370) * axis) + + def test_unwraps_across_pi(self): + """Crossing 180 degrees must continue counting up instead of + flipping the axis.""" + axis = np.array([1.0, 0.0, 0.0]) + R = Rotation.from_rotvec(np.radians(181) * axis) # wraps to -179 + reference = np.radians(180) * axis + out = render.closest_rotvec(R, reference) + assert np.allclose(out, np.radians(181) * axis) + + def test_identity_keeps_turns_of_reference(self): + """The identity rotation with a reference of ~2 turns must keep + the full turns along the reference axis.""" + axis = np.array([0.0, 1.0, 0.0]) + out = render.closest_rotvec( + Rotation.identity(), np.radians(719) * axis + ) + assert np.allclose(out, np.radians(720) * axis) + + def test_identity_zero_reference(self): + out = render.closest_rotvec(Rotation.identity(), np.zeros(3)) + assert np.allclose(out, np.zeros(3)) + + def test_result_represents_same_rotation(self): + R = Rotation.from_rotvec([0.2, -0.4, 0.6]) + reference = np.array([2.5, -4.0, 6.5]) + out = render.closest_rotvec(R, reference) + assert np.allclose( + (Rotation.from_rotvec(out) * R.inv()).magnitude(), 0.0, atol=1e-9 + ) + # --------------------------------------------------------------------------- # 3x3 matrix helpers @@ -1366,20 +1460,137 @@ def test_return_bgra(self): # --------------------------------------------------------------------------- +class TestAnimationSequence: + def test_slerp_endpoints_and_midpoint(self): + """Default interpolation is slerp: endpoints exact, midpoint at + half the rotation angle.""" + R1 = Rotation.identity() + R2 = Rotation.from_rotvec([0.0, 0.0, np.pi / 2]) + positions = [(R1, FULL_VIEWPORT), (R2, FULL_VIEWPORT)] + rotations, viewports = render._animation_sequence( + positions, [1.0], fps=3 + ) + assert len(rotations) == 3 + assert len(viewports) == 3 + assert (rotations[0] * R1.inv()).magnitude() == pytest.approx( + 0.0, abs=1e-9 + ) + assert (rotations[-1] * R2.inv()).magnitude() == pytest.approx( + 0.0, abs=1e-9 + ) + assert np.allclose( + rotations[1].as_rotvec(), [0.0, 0.0, np.pi / 4], atol=1e-9 + ) + + def test_full_turn_segment(self): + """A 360-degree segment rotation spins the full turn even though + both checkpoints have the same orientation.""" + R = Rotation.identity() + positions = [(R, FULL_VIEWPORT), (R, FULL_VIEWPORT)] + rotations, _ = render._animation_sequence( + positions, + [1.0], + fps=5, + segment_rotations=[np.array([0.0, 0.0, 2 * np.pi])], + ) + assert len(rotations) == 5 + # quarter of the way through: 90 degrees around z + assert np.allclose( + rotations[1].as_rotvec(), [0.0, 0.0, np.pi / 2], atol=1e-9 + ) + # halfway through: 180 degrees around z + assert rotations[2].magnitude() == pytest.approx(np.pi, abs=1e-9) + # ends exactly at the checkpoint again + assert rotations[-1].magnitude() == pytest.approx(0.0, abs=1e-9) + + def test_segment_rotation_snaps_to_checkpoints(self): + """The segment rotation vector is snapped so that the segment + ends exactly at the next checkpoint, keeping its turns.""" + R1 = Rotation.identity() + R2 = Rotation.from_rotvec([0.0, 0.0, np.pi / 2]) + positions = [(R1, FULL_VIEWPORT), (R2, FULL_VIEWPORT)] + # 450 degrees = 90 degrees + 1 full turn; pass a slightly off + # target to verify snapping + target = np.array([0.0, 0.0, np.radians(449)]) + rotations, _ = render._animation_sequence( + positions, [1.0], fps=11, segment_rotations=[target] + ) + assert (rotations[-1] * R2.inv()).magnitude() == pytest.approx( + 0.0, abs=1e-9 + ) + # total rotation path is 450 degrees: 1/5 of the way is 90 deg + assert np.allclose( + rotations[2].as_rotvec(), [0.0, 0.0, np.pi / 2], atol=1e-6 + ) + + def test_multi_segment_viewports(self): + """Viewports interpolate linearly per segment.""" + R = Rotation.identity() + vp1 = ((0.0, 0.0), (32.0, 32.0)) + vp2 = ((8.0, 8.0), (24.0, 24.0)) + positions = [(R, vp1), (R, vp2)] + _, viewports = render._animation_sequence(positions, [1.0], fps=3) + assert np.allclose(viewports[0], vp1) + assert np.allclose(viewports[1], ((4.0, 4.0), (28.0, 28.0))) + assert np.allclose(viewports[-1], vp2) + + def test_normalize_legacy_positions(self): + """Legacy Euler positions are converted with a deprecation + warning and match rotation_matrix.""" + legacy = [(0.1, 0.2, 0.3, FULL_VIEWPORT)] + with pytest.warns(DeprecationWarning): + normalized = render._normalize_animation_positions(legacy) + R, vp = normalized[0] + expected = render.rotation_matrix(0.1, 0.2, 0.3) + assert (R * expected.inv()).magnitude() == pytest.approx(0.0, abs=1e-9) + assert vp == FULL_VIEWPORT + + def test_normalize_invalid_position_raises(self): + with pytest.raises(ValueError): + render._normalize_animation_positions([(0.1, FULL_VIEWPORT)]) + + class TestBuildAnimation: - def test_smoke_two_frames(self, locs_3d, info, tmp_path): - """A 2-frame animation writes both an .mp4 and the sidecar .yaml.""" + def test_smoke_two_frames_legacy_euler(self, locs_3d, info, tmp_path): + """A 2-frame animation from legacy Euler positions writes both + an .mp4 and the sidecar .yaml (deprecated input format).""" out_path = tmp_path / "anim.mp4" positions = [ (0.0, 0.0, 0.0, FULL_VIEWPORT), (0.1, 0.0, 0.0, FULL_VIEWPORT), ] + with pytest.warns(DeprecationWarning): + render.build_animation( + str(out_path), + locs_3d, + info, + positions=positions, + durations=[1.0], + disp_px_size=PIXELSIZE, + image_size=(64, 64), + fps=2, + ) + assert out_path.exists() + assert out_path.stat().st_size > 0 + yaml_path = out_path.with_suffix(".yaml") + assert yaml_path.exists() + assert yaml_path.stat().st_size > 0 + + def test_smoke_quaternions_with_turns(self, locs_3d, info, tmp_path): + """An animation from scipy Rotations with a multi-turn segment + writes both an .mp4 and the sidecar .yaml.""" + out_path = tmp_path / "anim.mp4" + positions = [ + (Rotation.identity(), FULL_VIEWPORT), + (Rotation.from_rotvec([0.1, 0.0, 0.0]), FULL_VIEWPORT), + ] render.build_animation( str(out_path), locs_3d, info, positions=positions, durations=[1.0], + segment_rotations=[np.array([0.1 + 2 * np.pi, 0.0, 0.0])], disp_px_size=PIXELSIZE, image_size=(64, 64), fps=2, diff --git a/tests/test_zfit.py b/tests/test_zfit.py index cb0dd808..b6d0d588 100644 --- a/tests/test_zfit.py +++ b/tests/test_zfit.py @@ -497,6 +497,221 @@ def test_writes_yaml_and_png_when_path_given( assert loaded["Step size in nm"] == calib["Step size in nm"] +class TestCalibrateZFrameBounds: + """``frame_bounds`` restricts which frames enter the calibration. + Bounds are inclusive on both ends, matching ``picasso.localize``.""" + + N_FRAMES = 50 + D = 10.0 # nm step + + @pytest.fixture(autouse=True) + def _no_show(self, monkeypatch): + monkeypatch.setattr("matplotlib.pyplot.show", lambda *a, **k: None) + + @pytest.fixture + def bead_stack(self): + """Same synthetic bead stack as in ``TestCalibrateZ``: sx/sy + driven by known polynomials of stage z plus tiny noise.""" + rng = np.random.default_rng(0) + rows = [] + z_total = (self.N_FRAMES - 1) * self.D + for fi in range(self.N_FRAMES): + z = -(fi * self.D - z_total / 2) + sx_mean = 1.5 + 1e-3 * z + 1e-5 * z**2 + sy_mean = 1.5 - 1e-3 * z + 1e-5 * z**2 + for _ in range(40): + rows.append( + { + "frame": fi, + "x": 16.0, + "y": 16.0, + "sx": sx_mean + rng.normal(0, 0.02), + "sy": sy_mean + rng.normal(0, 0.02), + "photons": 5000.0, + "bg": 10.0, + "lpx": 0.01, + "lpy": 0.01, + } + ) + locs = pd.DataFrame(rows) + info = [ + { + "Frames": self.N_FRAMES, + "Pixelsize": 130, + "Width": 32, + "Height": 32, + } + ] + return locs, info + + def _polyfit_lengths(self, monkeypatch): + """Spy on ``np.polyfit`` to record how many frames enter each + calibration fit.""" + lengths = [] + orig = np.polyfit + + def spy(x, y, deg, **kwargs): + lengths.append(len(x)) + return orig(x, y, deg, **kwargs) + + monkeypatch.setattr(np, "polyfit", spy) + return lengths + + def test_none_and_none_none_equivalent(self, bead_stack): + """``frame_bounds=(None, None)`` must behave exactly like + ``frame_bounds=None`` (all frames used).""" + locs, info = bead_stack + calib_none = zfit.calibrate_z( + locs, info, self.D, magnification_factor=0.79, frame_bounds=None + ) + calib_nn = zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(None, None), + ) + np.testing.assert_allclose( + calib_none["X Coefficients"], calib_nn["X Coefficients"] + ) + np.testing.assert_allclose( + calib_none["Y Coefficients"], calib_nn["Y Coefficients"] + ) + + def test_full_range_bounds_equivalent_to_none(self, bead_stack): + """``(0, n_frames - 1)`` covers all frames, so it must match the + unbounded calibration.""" + locs, info = bead_stack + calib_none = zfit.calibrate_z( + locs, info, self.D, magnification_factor=0.79 + ) + calib_full = zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(0, self.N_FRAMES - 1), + ) + np.testing.assert_allclose( + calib_none["X Coefficients"], calib_full["X Coefficients"] + ) + np.testing.assert_allclose( + calib_none["Y Coefficients"], calib_full["Y Coefficients"] + ) + + def test_bounds_are_inclusive(self, bead_stack, monkeypatch): + """``(10, 39)`` keeps frames 10..39 inclusive -> 30 frames enter + every polynomial fit.""" + locs, info = bead_stack + lengths = self._polyfit_lengths(monkeypatch) + zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(10, 39), + ) + assert lengths and all(n == 30 for n in lengths) + + def test_one_sided_bounds(self, bead_stack, monkeypatch): + """``(None, max)`` and ``(min, None)`` each leave the other side + unbounded.""" + locs, info = bead_stack + lengths = self._polyfit_lengths(monkeypatch) + zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(None, 29), # frames 0..29 + ) + assert lengths and all(n == 30 for n in lengths) + lengths.clear() + zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(20, None), # frames 20..49 + ) + assert lengths and all(n == 30 for n in lengths) + + def test_bounds_including_last_frame(self, bead_stack): + """Regression test: a nonzero minimum together with the last + frame as maximum used to raise IndexError (localizations' + frame numbers were used as indices into the sliced per-frame + arrays without offsetting).""" + locs, info = bead_stack + calib = zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(10, self.N_FRAMES - 1), + ) + assert len(calib["X Coefficients"]) == 7 + assert len(calib["Y Coefficients"]) == 7 + + def test_locs_outside_bounds_are_ignored(self, bead_stack): + """Passing the full locs with bounds must equal passing locs + already restricted to those frames.""" + locs, info = bead_stack + bounds = (10, 39) + pre_filtered = locs[ + (locs["frame"] >= bounds[0]) & (locs["frame"] <= bounds[1]) + ] + calib_full = zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=bounds, + ) + calib_pre = zfit.calibrate_z( + pre_filtered, + info, + self.D, + magnification_factor=0.79, + frame_bounds=bounds, + ) + np.testing.assert_allclose( + calib_full["X Coefficients"], calib_pre["X Coefficients"] + ) + np.testing.assert_allclose( + calib_full["Y Coefficients"], calib_pre["Y Coefficients"] + ) + + def test_bounded_calibration_differs_from_unbounded(self, bead_stack): + """Restricting the frames must actually change the fit.""" + locs, info = bead_stack + calib_none = zfit.calibrate_z( + locs, info, self.D, magnification_factor=0.79 + ) + calib_bounded = zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(10, 39), + ) + assert not np.allclose( + calib_none["X Coefficients"], calib_bounded["X Coefficients"] + ) + + def test_frame_bounds_stored_in_calibration(self, bead_stack): + locs, info = bead_stack + calib = zfit.calibrate_z( + locs, + info, + self.D, + magnification_factor=0.79, + frame_bounds=(10, 39), + ) + assert calib["Frame bounds"] == (10, 39) + calib = zfit.calibrate_z(locs, info, self.D, magnification_factor=0.79) + assert calib["Frame bounds"] is None + + # --------------------------------------------------------------------------- # locs_from_futures # ---------------------------------------------------------------------------