From 0a7a4ed6e8ed092b8c31ab8cb517e241f643473e Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Wed, 4 Mar 2026 14:08:07 +0100 Subject: [PATCH 01/11] [feat] add peak-finding algorithm --- src/odemis/acq/align/goffset.py | 199 ++++++++++++++++++ src/odemis/acq/align/goffset_alignment.py | 193 +++++++++++++++++ .../acq/align/test/goffset_alignment_test.py | 125 +++++++++++ src/odemis/acq/align/test/goffset_test.py | 121 +++++++++++ 4 files changed, 638 insertions(+) create mode 100644 src/odemis/acq/align/goffset.py create mode 100644 src/odemis/acq/align/goffset_alignment.py create mode 100644 src/odemis/acq/align/test/goffset_alignment_test.py create mode 100644 src/odemis/acq/align/test/goffset_test.py diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py new file mode 100644 index 0000000000..cfa2350f5d --- /dev/null +++ b/src/odemis/acq/align/goffset.py @@ -0,0 +1,199 @@ +import logging +import threading +import time +from collections.abc import Iterable +from concurrent.futures import TimeoutError, CancelledError +from concurrent.futures._base import CANCELLED, FINISHED, RUNNING +from typing import Any, Dict, List, Optional, Tuple, Union, Callable + +import numpy + +from odemis import model +from odemis.acq.align import light +from odemis.model import InstantaneousFuture +from odemis.util import executeAsyncTask, almost_equal +from odemis.util.driver import guessActuatorMoveDuration +from odemis.util.focus import MeasureSEMFocus, Measure1d, MeasureSpotsFocus, AssessFocus +from odemis.util.img import Subtract +from scipy.optimize import curve_fit + +def Gaussian(x, amplitude, x0, width): + """ + x = input coordinates + amplitude = peak intensity + x0 = peak's center + width = standard deviation + """ + intensity = amplitude*numpy.exp(-0.5*((x-x0)/width)**2) + return intensity + + +def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: + """ + Finds the peak position in the given spectrum data. + It can handle both 1D and 2D data (in which case it averages over the first dimension). + """ + + if data.ndim == 2: + spectrum = data.mean(axis=0) # squash data into a 1D array + else: + spectrum = data + + x = numpy.arange(len(spectrum)) + peak_idx = numpy.argmax(spectrum) # find the absolute highest point + + # create a window around the peak + start = max(0, peak_idx - window_radius) + end = min(len(spectrum), peak_idx + window_radius + 1) + window_data = spectrum[start:end] + window_idx = numpy.arange(start, end) + + # calculate the weighted averages + weights = window_data.clip(min=0) + + if weights.sum() == 0: + weighted_avg = float(peak_idx) + else: + weighted_avg = float(numpy.sum(window_idx * window_data) / numpy.sum(window_data)) + + # try Gaussian fit for better accuracy, but fallback to weighted average if it fails or is out of bounds + + try: + p0 = [window_data.max(), weighted_avg, 2.5] + popt, pcov = curve_fit(Gaussian, window_idx, window_data, p0=p0) + peak = popt[1] + + if start <= peak <= end: + return float(peak) + + except Exception: + pass + + return weighted_avg + + +def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0) -> float: + """ + Estimates how many pixels the peak shifts per 1 unit of goffset. + """ + # get initial state + data0 = detector.data.get(asap=False) + p0 = find_peak_position(data0) + + # check limits before moving + current_pos = spgr.position.value["goffset"] + goffset_max = spgr.axes["goffset"].range[1] + + # ensure that max limit isn't violated + move_direction = 1 if (current_pos + delta < goffset_max) else -1 + actual_delta = delta*move_direction + + # move and measure + spgr.moveRelSync({"goffset": actual_delta}) + data1 = detector.data.get(asap=False) + p1 = find_peak_position(data1) + + # return back to start + spgr.moveRelSync({"goffset": -actual_delta}) + + # calculate goffset scale + scale = (p1-p0)/actual_delta + + logging.info( + f"SCALE TRACKING | p0: {p0:.1f} | p1: {p1:.1f} | Delta: {actual_delta} | Shift: {(p1-p0):.1f} | Result Scale: {scale:.4f}") + + if abs(scale) < 1e-3: + try: + scale = estimate_goffset_scale(spgr, detector) + except RuntimeError: + logging.warning("Scale too small, using default 0.5") + scale = 0.5 + + return scale + +def SparcAutoGratingOffset(spgr: model.Actuator, + detector: model.Detector, + align_grating: bool = True, + tolerance_px: float = 0.4, + max_it: int = 20, + gain: float = 0.4) -> model.ProgressiveFuture: + + est_start = time.time() + 0.05 + est_time = max_it*0.5 # conservative estimate + + f = model.ProgressiveFuture(start=est_start, end=est_start + est_time) + f._centering_state = RUNNING + f._centering_lock = threading.Lock() + f.task_canceller = _CancelSparcAutoGratingOffset + + executeAsyncTask(f, _DoSparcAutoGratingOffset, args=(f, spgr, detector, align_grating, tolerance_px, max_it, gain)) + + return f + +def _DoSparcAutoGratingOffset(future: model.ProgressiveFuture, + spgr: model.Actuator, + detector: model.Detector, + align_grating: bool, + tolerance_px: float, + max_it: int, + gain: float) -> bool: + + """ + Iteratively adjusts the grating offset to align the spectral peak to the center of the detector. + The algorithm estimates the current peak position, calculates the error from the center, and moves the grating offset proportionally to reduce this error. + It continues until the peak is within the specified tolerance or the maximum number of iterations is reached. + """ + + success = False + + logging.info("Running alignment | detector=%s | align_grating=%s",detector.name, align_grating,) + + try: + scale = estimate_goffset_scale(spgr, detector) + center_target = detector.resolution.value[0]/2 # adjust if 0 is not the center + total_goffset_displacement = 0.0 + + for i in range(max_it): + with future._centering_lock: + if future._centering_state == CANCELLED: + raise CancelledError() + + # find peak and calculate error + peak_px = find_peak_position(detector.data.get(asap=False)) + error_px = peak_px - center_target + + if abs(error_px) <= tolerance_px: + logging.info("Spectral peak aligned to 0 px after %d iterations", i + 1) + success = True + return True + + delta_goffset = -gain*(error_px/scale) + + # clamp move to safe fraction of axis range + axis = spgr.axes["goffset"] + minv, maxv = axis.range + max_step = 0.1*(maxv-minv) # max 10% of range + + delta_goffset = max(-max_step, min(max_step, delta_goffset)) + total_goffset_displacement += delta_goffset + + print(f"DEBUG | Iter: {i} | Peak: {peak_px:.1f} | Error: {error_px:.1f} | Move: {delta_goffset:.4f} | Total Change: {total_goffset_displacement:.4f}") + spgr.moveRelSync({"goffset": delta_goffset}) + time.sleep(2) + + future.set_progress( + end=time.time() + (max_it-i-1)*0.5) # update estimated end time + + logging.warning("SparcAutoGratingOffset did not converge") + return False + + except CancelledError: + logging.debug("SparcAutoGratingOffset cancelled") + raise + except Exception as e: + logging.error(f"Alignment error: {e}") + raise + +def _CancelSparcAutoGratingOffset(future: model.ProgressiveFuture): + with future._centering_lock: + future._centering_state = CANCELLED \ No newline at end of file diff --git a/src/odemis/acq/align/goffset_alignment.py b/src/odemis/acq/align/goffset_alignment.py new file mode 100644 index 0000000000..6559912ed0 --- /dev/null +++ b/src/odemis/acq/align/goffset_alignment.py @@ -0,0 +1,193 @@ +import logging +import threading +import time +from collections.abc import Iterable +from concurrent.futures import TimeoutError, CancelledError +from concurrent.futures._base import CANCELLED, FINISHED, RUNNING +from typing import Any, Dict, List, Optional, Tuple, Union, Callable + +from odemis import model +from odemis.model import InstantaneousFuture +from odemis.util import executeAsyncTask, almost_equal + +from odemis.acq.align.goffset import SparcAutoGratingOffset + + +def _mapDetectorToSelector(selector: model.Actuator, + detectors: List[model.Detector]) -> Tuple[str, Dict[str, Any]]: + det_2_sel = {} + sel_axis = None + + for an, ad in selector.axes.items(): + if hasattr(ad, "choices") and isinstance(ad.choices, dict): + sel_axis = an + for pos, value in ad.choices.items(): + for d in detectors: + if d.name in value: + det_2_sel[d] = pos + + if det_2_sel: + break + + if len(det_2_sel) < len(detectors): + raise ValueError("Failed to find all detectors (%s) in positions of selector axes %s" % + (", ".join(d.name for d in detectors), list(selector.axes.keys()))) + + return sel_axis, det_2_sel + +def _checkCancelled(future: "model.ProgressiveFuture"): + lock_name = "_centering_lock" if hasattr(future, "_centering_lock") else "_align_lock" + state_name = "_centering_state" if hasattr(future, "_centering_state") else "_align_state" + + lock = getattr(future, lock_name) + state = getattr(future, state_name) + + with lock: + if state == CANCELLED: + raise CancelledError() + + +def _totalAlignmentTime(n_gratings: int, + n_detectors: int) -> float: + + """ + Estimate total time for aligning all grating-detector combinations. + """ + + runs = n_detectors + max(0, n_gratings - 1) + move_time = ((n_gratings-1)*MOVE_TIME_GRATING+(n_detectors-1)*MOVE_TIME_DETECTOR) + + # total time = time spent running alignment algorithms + time spent moving hardware + return runs*EST_ALIGN_TIME + move_time + + +def AutoAlignGratingDetectorOffsets(spectrograph: model.Actuator, + detectors: Union[model.Detector, List[model.Detector]], + selector: Optional[model.Actuator] = None, + streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: + + if not isinstance(detectors, Iterable): + detectors = [detectors] + if not detectors: + raise ValueError("At least one detector must be provided") + if len(detectors) > 1 and selector is None: + raise ValueError("No selector provided, but multiple detectors") + + if streams is None: + streams = [] + + est_start = time.time() + 0.1 + n_gratings = len(spectrograph.axes["grating"].choices) + n_detectors = len(detectors) + a_time = _totalAlignmentTime(n_gratings, n_detectors) + f = model.ProgressiveFuture(start=est_start, end=est_start + a_time) + f.task_canceller = _CancelAutoAlignGratingDetectorOffsets + + f._align_state = RUNNING + f._align_lock = threading.Lock() + f._subfuture = InstantaneousFuture() + executeAsyncTask(f, _DoAutoAlignGratingDetectorOffsets, args=(f, spectrograph, detectors, selector, streams)) + return f + +MOVE_TIME_GRATING = 20 #s +MOVE_TIME_DETECTOR = 5 #s +EST_ALIGN_TIME = 30 #s + +def _DoAutoAlignGratingDetectorOffsets(future: model.ProgressiveFuture, + spectrograph: model.Actuator, + detectors: List[model.Detector], + selector: Optional[model.Actuator], + streams: List['Stream'], + stabilization_time: float = 15.0) -> Optional[Dict[Any, Any]]: + + """ + Iterate through each grating and detector combination, adjusting the selector if provided, and run the auto-alignment algorithm. + - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector will be used for all subsequent gratings. + - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will be aligned by adjusting the detector offset with the grating alignment fixed. + """ + + results: dict[tuple, bool] = {} + original_pos = {k: v for k, v in spectrograph.position.value.items() + if k in ("wavelength", "grating")} + + gratings = sorted(list(spectrograph.axes["grating"].choices.keys())) + logging.info(f"Available gratings: {list(spectrograph.axes['grating'].choices.keys())}") + + first_detector = detectors[0] + + if selector: + original_selector = selector.position.value + selector_axes, detector_to_selector = _mapDetectorToSelector(selector, detectors) + + def is_current_detector(d): + if selector is None: + return True + return detector_to_selector[d] == selector.position.value[selector_axes] + + try: + g0 = gratings[0] + logging.info(f"Starting alignment for initial grating: {g0}") + + spectrograph.moveAbsSync({"grating": g0, "wavelength": 0}) + time.sleep(stabilization_time) + + detectors_sorted = sorted(detectors, key=is_current_detector, reverse=True) + + # align each detector for the first grating + for d in detectors_sorted: + _checkCancelled(future) + logging.info(f"Starting alignment | Detector: {d.name} | Grating: {g0}") + + if selector: + selector.moveAbsSync({selector_axes: detector_to_selector[d]}) + align_grating = (d is first_detector) + future._subfuture = SparcAutoGratingOffset(spectrograph, d, align_grating = align_grating,) + success = future._subfuture.result() + results[(g0, d.name)] = success + + logging.info(f"Finished alignment | Detector: {d.name} | Grating: {g0} | Success: {success}") + + if selector: + selector.moveAbsSync({selector_axes: detector_to_selector[first_detector]}) + + # align remaining gratings using the first detector + for g in gratings[1:]: + _checkCancelled(future) + logging.info(f"Switching to grating: {g}") + + spectrograph.moveAbsSync({"grating": g, "wavelength": 0}) + time.sleep(stabilization_time) + logging.info(f"Starting alignment | Detector: {first_detector.name} | Grating: {g}") + + future._subfuture = SparcAutoGratingOffset(spectrograph, first_detector, align_grating=True,) + success = future._subfuture.result() + results[(g, first_detector.name)] = success + + logging.info(f"Finished alignment | Detector: {first_detector.name} | Grating: {g} | Success: {success}") + + return results + + except CancelledError: + logging.info("Auto-alignment cancelled") + raise + + finally: + spectrograph.moveAbsSync(original_pos) + if selector: + selector.moveAbsSync(original_selector) + + with future._align_lock: + future._align_state = FINISHED + + +def _CancelAutoAlignGratingDetectorOffsets(future: model.ProgressiveFuture) -> bool: + logging.debug("Cancelling autoalignment...") + + with future._align_lock: + if future._align_state == FINISHED: + return False + future._align_state = CANCELLED + future._subfuture.cancel() + logging.debug("Auto-alignment cancellation requested") + + return True \ No newline at end of file diff --git a/src/odemis/acq/align/test/goffset_alignment_test.py b/src/odemis/acq/align/test/goffset_alignment_test.py new file mode 100644 index 0000000000..d36e6bc49d --- /dev/null +++ b/src/odemis/acq/align/test/goffset_alignment_test.py @@ -0,0 +1,125 @@ +import logging +import os +import time +import unittest +import numpy +from collections.abc import Iterable +from concurrent.futures import CancelledError +from scipy import ndimage + +from odemis import model, acq +import odemis +from odemis.acq import align, stream, path +from odemis.acq.align.autofocus import Sparc2AutoFocus, MTD_BINARY +from odemis.acq.align.goffset_alignment import AutoAlignGratingDetectorOffsets +from odemis.dataio import tiff, hdf5 +from odemis.util import testing, timeout, img +import odemis.util.focus +from unittest.mock import patch + + +CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" +SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" + +logging.getLogger().setLevel(logging.DEBUG) + +class TestAutoAlignGratingDetectorOffsets(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.spgr = model.getComponent(role="spectrograph") + cls.ccd = model.getComponent(role="ccd") + #cls.spccd = model.getComponent(role="sp-ccd") + cls.selector = model.getComponent(role="spec-det-selector") + + def setUp(self): + # Speed up acquisition + self.ccd.exposureTime.value = self.ccd.exposureTime.range[0] + + @timeout(100) + def test_cancel(self): + f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=[self.ccd],) + time.sleep(1) + + cancelled = f.cancel() + self.assertTrue(cancelled) + self.assertTrue(f.cancelled()) + + with self.assertRaises(CancelledError): + f.result(timeout=900) + + + @timeout(1000) + def test_single_detector_iteration(self): + f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + res = f.result(timeout=900) + + n_gratings = len(self.spgr.axes["grating"].choices) + n_detectors = 1 + expected = n_detectors + (n_gratings - 1) + + self.assertEqual(len(res), expected) + + first_grating = list(self.spgr.axes["grating"].choices.keys())[0] + dets_first = [d for (g, d) in res.keys() if g == first_grating] + + self.assertEqual(len(dets_first), n_detectors) + + @timeout(100) + def test_cancel(self): + f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=[self.ccd],) + + time.sleep(1) + + cancelled = f.cancel() + self.assertTrue(cancelled) + self.assertTrue(f.cancelled()) + + with self.assertRaises(CancelledError): + f.result(timeout=900) + + def test_multi_detector_iteration(self): + spccd = model.getComponent(role="sp-ccd") + spccd.exposureTime.value = spccd.exposureTime.range[0] + + detectors = [self.ccd, spccd] + + # run alignment + f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) + res = f.result(timeout=900) + + # calculate expected results + n_gratings = len(self.spgr.axes["grating"].choices) + n_detectors = len(detectors) + expected_count = n_detectors + (n_gratings - 1) + + self.assertEqual(len(res), expected_count, f"Expected {expected_count} results, got {len(res)}") + + # verify that every detector was used for the first grating + gratings_list = list(self.spgr.axes["grating"].choices.keys()) + first_grating = gratings_list[0] + + dets_for_first_grating = [d for (g, d) in res.keys() if g == first_grating] + self.assertEqual(len(dets_for_first_grating), n_detectors) + self.assertIn(self.ccd.name, dets_for_first_grating) + self.assertIn(spccd.name, dets_for_first_grating) + + # verify that only first detector is used for remaining gratings + for g in gratings_list[1:]: + # Ensure only the first detector (index 0) is present for these gratings + dets_for_this_grating = [d for (grating, d) in res.keys() if grating == g] + self.assertEqual(len(dets_for_this_grating), 1) + self.assertEqual(dets_for_this_grating[0], detectors[0].name) + + # move to spectral camera + self.selector.moveAbsSync({"rx": 1.5707963267948966}) + data = spccd.data.get(asap=False) + + # check data is not flat + if data.max() == data.min(): + print("WARNING: sp-ccd is returning a flat image!") + else: + print(f"sp-ccd signal range: {data.min()} to {data.max()}") + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py new file mode 100644 index 0000000000..958d44b1bb --- /dev/null +++ b/src/odemis/acq/align/test/goffset_test.py @@ -0,0 +1,121 @@ +# -*- coding: utf-8 -*- +""" +Test Sparc auto grating offset alignment +""" + +import os +import time +import unittest +import logging +import numpy as np + +from concurrent.futures._base import CancelledError + +from odemis import model +from odemis.util import testing, timeout +from odemis.acq.align.goffset import ( + find_peak_position, + estimate_goffset_scale, + SparcAutoGratingOffset, +) + +import odemis + +logging.getLogger().setLevel(logging.DEBUG) + +CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" +SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" + + +class TestSparcAutoGratingOffset(unittest.TestCase): + """ + Test automatic grating offset alignment. + """ + + @classmethod + def setUpClass(cls): + #testing.start_backend(SPARC_CONFIG) + + cls.detector = model.getComponent(role="ccd") + cls.spgr = model.getComponent(role="spectrograph") + + cls._original_goffset = cls.spgr.goffset + cls._original_position = cls.spgr.position.value.copy() + + # @classmethod + # def tearDownClass(cls): + # # restore original position + # try: + # cls.spgr.moveAbsSync(cls._original_position) + # except Exception: + # logging.exception("Failed restoring spectrograph position") + + def setUp(self): + # speed up detector + self.detector.exposureTime.value = self.detector.exposureTime.range[0] + + def test_find_peak_position_synthetic(self): + """ + Test peak detection on synthetic Gaussian data. + """ + x = np.arange(200) + true_center = 83.4 + spectrum = np.exp(-0.5*((x-true_center)/3.0)**2) + + peak = find_peak_position(spectrum) + self.assertAlmostEqual(peak, true_center, places=1) + + def test_find_peak_position_2d(self): + """ + Test peak detection on 2D data (mean over axis). + """ + x = np.arange(200) + true_center = 120.0 + line = np.exp(-0.5*((x-true_center)/4.0)**2) + image = np.tile(line, (50, 1)) + + peak = find_peak_position(image) + self.assertAlmostEqual(peak, true_center, places=1) + + @timeout(100) + def test_estimate_goffset_scale(self): + """ + Test that goffset scale is non-zero and finite. + """ + scale = estimate_goffset_scale(self.spgr, self.detector) + + self.assertIsInstance(scale, float) + self.assertNotEqual(scale, 0.0) + self.assertTrue(np.isfinite(scale)) + + @timeout(300) + def test_auto_grating_offset(self): + """ + Test automatic centering of spectral peak. + """ + delta = 0 # intentionally misalign + current = self.spgr.position.value["goffset"] + goffset_max = self.spgr.axes["goffset"].range[1] + direction = 1 if (current + delta < goffset_max) else -1 + + self.spgr.moveRelSync({"goffset": delta * direction}) + f = SparcAutoGratingOffset(self.spgr, self.detector, max_it=20) + + result = f.result(timeout=200) + self.assertTrue(result) + + @timeout(100) + def test_cancel(self): + """ + Test cancelling alignment. + """ + f = SparcAutoGratingOffset(self.spgr, self.detector) + # Wait for the result or a timeout + try: + f.result(timeout=5) + except: + pass + self.assertTrue(f.done()) + +if __name__ == "__main__": + unittest.main() \ No newline at end of file From bd6d06f26477fe351fd0d3b4e6a0ab7a5c05f43d Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Wed, 4 Mar 2026 14:08:07 +0100 Subject: [PATCH 02/11] [feat] add peak-finding algorithm --- src/odemis/acq/align/test/goffset_test.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index 958d44b1bb..7e20aade80 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -13,11 +13,10 @@ from odemis import model from odemis.util import testing, timeout -from odemis.acq.align.goffset import ( - find_peak_position, + +from odemis.acq.align.goffset import (find_peak_position, estimate_goffset_scale, - SparcAutoGratingOffset, -) + sparc_auto_grating_offset) import odemis From f75f0f036d9633a6e04b7dc1742167046db19d43 Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Mon, 9 Mar 2026 09:44:36 +0100 Subject: [PATCH 03/11] [fix] incorporated feedback --- src/odemis/acq/align/goffset.py | 165 +++++++++++++----- src/odemis/acq/align/goffset_alignment.py | 52 ++---- .../acq/align/test/goffset_alignment_test.py | 10 +- 3 files changed, 139 insertions(+), 88 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index cfa2350f5d..c01aa00c94 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -17,14 +17,17 @@ from odemis.util.img import Subtract from scipy.optimize import curve_fit -def Gaussian(x, amplitude, x0, width): +def gaussian(x, amplitude, x0, width): """ - x = input coordinates - amplitude = peak intensity - x0 = peak's center - width = standard deviation + Gaussian function (for curve fitting). + + :param x: input coordinates + :param amplitude: peak intensity + :param x0 = peak's center + :param width = standard deviation + :return: Gaussian function evaluated at x """ - intensity = amplitude*numpy.exp(-0.5*((x-x0)/width)**2) + intensity = amplitude*numpy.exp(-0.5 * ((x - x0) / width)**2) return intensity @@ -32,6 +35,17 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: """ Finds the peak position in the given spectrum data. It can handle both 1D and 2D data (in which case it averages over the first dimension). + + The function first identifies the absolute peak, then creates a window around it to minimize noise influence. + It then calculates a weighted average of the positions within this window, using the intensity values as weights. + For improved accuracy, it attempts to fit a Gaussian curve to the windowed data, estimating the peak's center, width + and amplitude. + + However, the curve fit can fail to converge or produce unreasonable results if the initial guess is poor, if there are + outliers, baseline trends, or if the data is too noisy. In such cases, the function falls back to using the weighted + average as the peak position estimate. The weighted average is more robust as it does not run an iterative optimizer, + so it has no convergence or numerical-optimization failure modes, but it may be less accurate if the peak is not well-defined + or if there are multiple peaks within the window. """ if data.ndim == 2: @@ -51,30 +65,55 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: # calculate the weighted averages weights = window_data.clip(min=0) + # Corner case: window has no positive signal. Use the maximum index itself + # as the best estimate of the peak, since weighted average would be undefined if weights.sum() == 0: weighted_avg = float(peak_idx) + logging.info("Weighted average fallback: all window data <= 0, using peak_idx=%d as estimate", + peak_idx) else: weighted_avg = float(numpy.sum(window_idx * window_data) / numpy.sum(window_data)) # try Gaussian fit for better accuracy, but fallback to weighted average if it fails or is out of bounds try: - p0 = [window_data.max(), weighted_avg, 2.5] - popt, pcov = curve_fit(Gaussian, window_idx, window_data, p0=p0) + p0 = [window_data.max(), weighted_avg, 2.5] # intial guess: [amplitude, center, width] + popt, pcov = curve_fit(gaussian, window_idx, window_data, p0=p0) peak = popt[1] if start <= peak <= end: return float(peak) - except Exception: - pass + except RuntimeError: + logging.info("Gaussian peak fit did not converge, falling back to weighted average") + + except ValueError: + logging.exception("Gaussian peak fit failed due to invalid input") return weighted_avg def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0) -> float: """ - Estimates how many pixels the peak shifts per 1 unit of goffset. + Estimate the scale factor between a change in the grating offset ('goffset') + and the resulting shift of the spectral peak on the detector. + + The function moves the actuator by a small step (delta) and measures the peak position before and after the move. + It calculates the ratio of pixel shift per unit of goffset. The actuator is returned to its original position + after measurement. + + If the measured scale is unreasonably small or large, the function retries recursively + and falls back to a default value of 0.5 if necessary. + + :param spgr: spectrograph + :param detector: detector + :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). + The actual step may be negated to avoid exceeding hardware limits. + + :return: Tuple (scale, p0, p1) + scale: estimated pixels per unit of goffset + p0: peak position at the initial goffset + p1: peak position after applying the test delta """ # get initial state data0 = detector.data.get(asap=False) @@ -97,92 +136,122 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta spgr.moveRelSync({"goffset": -actual_delta}) # calculate goffset scale - scale = (p1-p0)/actual_delta + scale = (p1 - p0) / actual_delta logging.info( f"SCALE TRACKING | p0: {p0:.1f} | p1: {p1:.1f} | Delta: {actual_delta} | Shift: {(p1-p0):.1f} | Result Scale: {scale:.4f}") - if abs(scale) < 1e-3: + + # If the estimated scale is extremely small, the measurement is likely unreliable + # (e.g., due to noise or a flat spectrum). + # In that case we retry the estimation once by calling the function recursively. + # If the retry still fails (raising a RuntimeError), we fall back to a default + # scale value to ensure the algorithm can continue and avoid infinite recursion. + + if abs(scale) < 1e-3 or abs(scale)>10.0: try: - scale = estimate_goffset_scale(spgr, detector) + scale, p0, p1 = estimate_goffset_scale(spgr, detector) except RuntimeError: logging.warning("Scale too small, using default 0.5") scale = 0.5 - return scale + return scale, p0, p1 -def SparcAutoGratingOffset(spgr: model.Actuator, +def sparc_auto_grating_offset(spgr: model.Actuator, detector: model.Detector, - align_grating: bool = True, tolerance_px: float = 0.4, max_it: int = 20, gain: float = 0.4) -> model.ProgressiveFuture: + """ + Start an asynchronous task that centers the spectral peak by adjusting the + grating offset (goffset). + + :param spgr: spectrograph + :param detector: detector + :param tolerance_px: the acceptable displacement of the peak from the center in pixels (default: 0.4) + :param max_it: maximum number of iterations to attempt (default: 20) + :param gain: proportional gain factor for adjusting the goffset (default: 0.4) + :return: A ``ProgressiveFuture`` representing the asynchronous alignment + task. The future can be used to monitor progress, retrieve the + result, or cancel the alignment. + """ est_start = time.time() + 0.05 est_time = max_it*0.5 # conservative estimate f = model.ProgressiveFuture(start=est_start, end=est_start + est_time) - f._centering_state = RUNNING - f._centering_lock = threading.Lock() - f.task_canceller = _CancelSparcAutoGratingOffset - executeAsyncTask(f, _DoSparcAutoGratingOffset, args=(f, spgr, detector, align_grating, tolerance_px, max_it, gain)) + f._task_lock = threading.Lock() + f._task_state = RUNNING + f.task_canceller = _cancel_sparc_auto_grating_offset + + executeAsyncTask( + f, + _do_sparc_auto_grating_offset, + args=(f, spgr, detector, tolerance_px, max_it, gain), + ) return f -def _DoSparcAutoGratingOffset(future: model.ProgressiveFuture, +def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, spgr: model.Actuator, detector: model.Detector, - align_grating: bool, tolerance_px: float, max_it: int, gain: float) -> bool: """ Iteratively adjusts the grating offset to align the spectral peak to the center of the detector. - The algorithm estimates the current peak position, calculates the error from the center, and moves the grating offset proportionally to reduce this error. + The algorithm estimates the current peak position, calculates the error from the center, and moves the grating offset + proportionally to reduce this error. It continues until the peak is within the specified tolerance or the maximum number of iterations is reached. """ - success = False - - logging.info("Running alignment | detector=%s | align_grating=%s",detector.name, align_grating,) + logging.info("Running alignment | detector=%s |", + detector.name) try: - scale = estimate_goffset_scale(spgr, detector) - center_target = detector.resolution.value[0]/2 # adjust if 0 is not the center + scale, p0, p1 = estimate_goffset_scale(spgr, detector) + center_target = detector.resolution.value[0] / 2 # adjust if 0 is not the center total_goffset_displacement = 0.0 + # clamp move to safe fraction of axis range + axis = spgr.axes["goffset"] + minv, maxv = axis.range + max_step = 0.1 * (maxv - minv) # max 10% of range + for i in range(max_it): - with future._centering_lock: - if future._centering_state == CANCELLED: + with future._task_lock: + if future._task_state == CANCELLED: raise CancelledError() - # find peak and calculate error - peak_px = find_peak_position(detector.data.get(asap=False)) + if i == 0: + peak_px = p0 + else: + data = detector.data.get(asap=False) + peak_px = find_peak_position(data) + error_px = peak_px - center_target if abs(error_px) <= tolerance_px: - logging.info("Spectral peak aligned to 0 px after %d iterations", i + 1) - success = True + logging.info( + "Spectral peak aligned after %d iterations | error_px=%.2f | total_goffset_change=%.4f", + i + 1, + error_px, + total_goffset_displacement + ) return True delta_goffset = -gain*(error_px/scale) - # clamp move to safe fraction of axis range - axis = spgr.axes["goffset"] - minv, maxv = axis.range - max_step = 0.1*(maxv-minv) # max 10% of range - delta_goffset = max(-max_step, min(max_step, delta_goffset)) total_goffset_displacement += delta_goffset - print(f"DEBUG | Iter: {i} | Peak: {peak_px:.1f} | Error: {error_px:.1f} | Move: {delta_goffset:.4f} | Total Change: {total_goffset_displacement:.4f}") + logging.debug(f"DEBUG | Iter: {i} | Peak: {peak_px:.1f} | Error: {error_px:.1f} | Move: {delta_goffset:.4f} | Total Change: {total_goffset_displacement:.4f}") spgr.moveRelSync({"goffset": delta_goffset}) - time.sleep(2) future.set_progress( - end=time.time() + (max_it-i-1)*0.5) # update estimated end time + end=time.time() + (max_it - i - 1) * 0.5) # update estimated end time logging.warning("SparcAutoGratingOffset did not converge") return False @@ -194,6 +263,10 @@ def _DoSparcAutoGratingOffset(future: model.ProgressiveFuture, logging.error(f"Alignment error: {e}") raise -def _CancelSparcAutoGratingOffset(future: model.ProgressiveFuture): - with future._centering_lock: - future._centering_state = CANCELLED \ No newline at end of file +def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): + + """ + Canceller of _do_sparc_auto_grating_offset task. + """ + with future._task_lock: + future._task_state = CANCELLED \ No newline at end of file diff --git a/src/odemis/acq/align/goffset_alignment.py b/src/odemis/acq/align/goffset_alignment.py index 6559912ed0..e48f07976e 100644 --- a/src/odemis/acq/align/goffset_alignment.py +++ b/src/odemis/acq/align/goffset_alignment.py @@ -4,37 +4,16 @@ from collections.abc import Iterable from concurrent.futures import TimeoutError, CancelledError from concurrent.futures._base import CANCELLED, FINISHED, RUNNING -from typing import Any, Dict, List, Optional, Tuple, Union, Callable +from typing import Any, Dict, List, Optional, Tuple, Union from odemis import model from odemis.model import InstantaneousFuture -from odemis.util import executeAsyncTask, almost_equal +from odemis.util import executeAsyncTask -from odemis.acq.align.goffset import SparcAutoGratingOffset +from odemis.acq.align.autofocus import _mapDetectorToSelector +from odemis.acq.align.goffset import sparc_auto_grating_offset -def _mapDetectorToSelector(selector: model.Actuator, - detectors: List[model.Detector]) -> Tuple[str, Dict[str, Any]]: - det_2_sel = {} - sel_axis = None - - for an, ad in selector.axes.items(): - if hasattr(ad, "choices") and isinstance(ad.choices, dict): - sel_axis = an - for pos, value in ad.choices.items(): - for d in detectors: - if d.name in value: - det_2_sel[d] = pos - - if det_2_sel: - break - - if len(det_2_sel) < len(detectors): - raise ValueError("Failed to find all detectors (%s) in positions of selector axes %s" % - (", ".join(d.name for d in detectors), list(selector.axes.keys()))) - - return sel_axis, det_2_sel - def _checkCancelled(future: "model.ProgressiveFuture"): lock_name = "_centering_lock" if hasattr(future, "_centering_lock") else "_align_lock" state_name = "_centering_state" if hasattr(future, "_centering_state") else "_align_state" @@ -47,7 +26,7 @@ def _checkCancelled(future: "model.ProgressiveFuture"): raise CancelledError() -def _totalAlignmentTime(n_gratings: int, +def _total_alignment_time(n_gratings: int, n_detectors: int) -> float: """ @@ -55,13 +34,13 @@ def _totalAlignmentTime(n_gratings: int, """ runs = n_detectors + max(0, n_gratings - 1) - move_time = ((n_gratings-1)*MOVE_TIME_GRATING+(n_detectors-1)*MOVE_TIME_DETECTOR) + move_time = ((n_gratings-1) * MOVE_TIME_GRATING + (n_detectors - 1) * MOVE_TIME_DETECTOR) # total time = time spent running alignment algorithms + time spent moving hardware - return runs*EST_ALIGN_TIME + move_time + return runs * EST_ALIGN_TIME + move_time -def AutoAlignGratingDetectorOffsets(spectrograph: model.Actuator, +def auto_align_grating_detector_offsets(spectrograph: model.Actuator, detectors: Union[model.Detector, List[model.Detector]], selector: Optional[model.Actuator] = None, streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: @@ -79,21 +58,21 @@ def AutoAlignGratingDetectorOffsets(spectrograph: model.Actuator, est_start = time.time() + 0.1 n_gratings = len(spectrograph.axes["grating"].choices) n_detectors = len(detectors) - a_time = _totalAlignmentTime(n_gratings, n_detectors) + a_time = _total_alignment_time(n_gratings, n_detectors) f = model.ProgressiveFuture(start=est_start, end=est_start + a_time) - f.task_canceller = _CancelAutoAlignGratingDetectorOffsets + f.task_canceller = _cancel_auto_align_grating_detector_offsets f._align_state = RUNNING f._align_lock = threading.Lock() f._subfuture = InstantaneousFuture() - executeAsyncTask(f, _DoAutoAlignGratingDetectorOffsets, args=(f, spectrograph, detectors, selector, streams)) + executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams)) return f MOVE_TIME_GRATING = 20 #s MOVE_TIME_DETECTOR = 5 #s EST_ALIGN_TIME = 30 #s -def _DoAutoAlignGratingDetectorOffsets(future: model.ProgressiveFuture, +def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, spectrograph: model.Actuator, detectors: List[model.Detector], selector: Optional[model.Actuator], @@ -140,8 +119,7 @@ def is_current_detector(d): if selector: selector.moveAbsSync({selector_axes: detector_to_selector[d]}) - align_grating = (d is first_detector) - future._subfuture = SparcAutoGratingOffset(spectrograph, d, align_grating = align_grating,) + future._subfuture = sparc_auto_grating_offset(spectrograph, d) success = future._subfuture.result() results[(g0, d.name)] = success @@ -159,7 +137,7 @@ def is_current_detector(d): time.sleep(stabilization_time) logging.info(f"Starting alignment | Detector: {first_detector.name} | Grating: {g}") - future._subfuture = SparcAutoGratingOffset(spectrograph, first_detector, align_grating=True,) + future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector) success = future._subfuture.result() results[(g, first_detector.name)] = success @@ -180,7 +158,7 @@ def is_current_detector(d): future._align_state = FINISHED -def _CancelAutoAlignGratingDetectorOffsets(future: model.ProgressiveFuture) -> bool: +def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) -> bool: logging.debug("Cancelling autoalignment...") with future._align_lock: diff --git a/src/odemis/acq/align/test/goffset_alignment_test.py b/src/odemis/acq/align/test/goffset_alignment_test.py index d36e6bc49d..5957cfa158 100644 --- a/src/odemis/acq/align/test/goffset_alignment_test.py +++ b/src/odemis/acq/align/test/goffset_alignment_test.py @@ -11,7 +11,7 @@ import odemis from odemis.acq import align, stream, path from odemis.acq.align.autofocus import Sparc2AutoFocus, MTD_BINARY -from odemis.acq.align.goffset_alignment import AutoAlignGratingDetectorOffsets +from odemis.acq.align.goffset_alignment import auto_align_grating_detector_offsets from odemis.dataio import tiff, hdf5 from odemis.util import testing, timeout, img import odemis.util.focus @@ -38,7 +38,7 @@ def setUp(self): @timeout(100) def test_cancel(self): - f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=[self.ccd],) + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd],) time.sleep(1) cancelled = f.cancel() @@ -51,7 +51,7 @@ def test_cancel(self): @timeout(1000) def test_single_detector_iteration(self): - f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) res = f.result(timeout=900) n_gratings = len(self.spgr.axes["grating"].choices) @@ -67,7 +67,7 @@ def test_single_detector_iteration(self): @timeout(100) def test_cancel(self): - f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=[self.ccd],) + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd],) time.sleep(1) @@ -85,7 +85,7 @@ def test_multi_detector_iteration(self): detectors = [self.ccd, spccd] # run alignment - f = AutoAlignGratingDetectorOffsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) res = f.result(timeout=900) # calculate expected results From 2f65bf7499474ed507f738f70354e5a47eb42d4d Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Mon, 9 Mar 2026 09:44:36 +0100 Subject: [PATCH 04/11] [fix] incorporated feedback --- src/odemis/acq/align/test/goffset_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index 7e20aade80..93a5fb99c2 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -98,7 +98,7 @@ def test_auto_grating_offset(self): direction = 1 if (current + delta < goffset_max) else -1 self.spgr.moveRelSync({"goffset": delta * direction}) - f = SparcAutoGratingOffset(self.spgr, self.detector, max_it=20) + f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=20) result = f.result(timeout=200) self.assertTrue(result) @@ -108,7 +108,7 @@ def test_cancel(self): """ Test cancelling alignment. """ - f = SparcAutoGratingOffset(self.spgr, self.detector) + f = sparc_auto_grating_offset(self.spgr, self.detector) # Wait for the result or a timeout try: f.result(timeout=5) From 1d430b91b7e39f2e3829fa3ee9d0f686bdf1ded3 Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Mon, 9 Mar 2026 14:19:48 +0100 Subject: [PATCH 05/11] [feat] merged goffset + goffset_alignment --- src/odemis/acq/align/goffset.py | 243 +++++++++++++++++++--- src/odemis/acq/align/goffset_alignment.py | 61 ++++-- 2 files changed, 264 insertions(+), 40 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index c01aa00c94..22a71e9072 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -2,19 +2,16 @@ import threading import time from collections.abc import Iterable -from concurrent.futures import TimeoutError, CancelledError +from concurrent.futures import CancelledError from concurrent.futures._base import CANCELLED, FINISHED, RUNNING -from typing import Any, Dict, List, Optional, Tuple, Union, Callable +from typing import Any, Dict, List, Optional, Tuple, Union import numpy from odemis import model -from odemis.acq.align import light +from odemis.acq.align.autofocus import _mapDetectorToSelector from odemis.model import InstantaneousFuture -from odemis.util import executeAsyncTask, almost_equal -from odemis.util.driver import guessActuatorMoveDuration -from odemis.util.focus import MeasureSEMFocus, Measure1d, MeasureSpotsFocus, AssessFocus -from odemis.util.img import Subtract +from odemis.util import executeAsyncTask from scipy.optimize import curve_fit def gaussian(x, amplitude, x0, width): @@ -53,7 +50,6 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: else: spectrum = data - x = numpy.arange(len(spectrum)) peak_idx = numpy.argmax(spectrum) # find the absolute highest point # create a window around the peak @@ -93,7 +89,7 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: return weighted_avg -def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0) -> float: +def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0, retries = 1) -> Tuple[float, float, float]: """ Estimate the scale factor between a change in the grating offset ('goffset') and the resulting shift of the spectral peak on the detector. @@ -109,12 +105,14 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta :param detector: detector :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). The actual step may be negated to avoid exceeding hardware limits. + :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). :return: Tuple (scale, p0, p1) scale: estimated pixels per unit of goffset p0: peak position at the initial goffset p1: peak position after applying the test delta """ + # get initial state data0 = detector.data.get(asap=False) p0 = find_peak_position(data0) @@ -139,8 +137,9 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta scale = (p1 - p0) / actual_delta logging.info( - f"SCALE TRACKING | p0: {p0:.1f} | p1: {p1:.1f} | Delta: {actual_delta} | Shift: {(p1-p0):.1f} | Result Scale: {scale:.4f}") - + "SCALE TRACKING | p0: %.1f | p1: %.1f | Delta: %.1f | Shift: %.1f | Result Scale: %.4f", + p0, p1, actual_delta, (p1 - p0), scale, + ) # If the estimated scale is extremely small, the measurement is likely unreliable # (e.g., due to noise or a flat spectrum). @@ -148,12 +147,18 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta # If the retry still fails (raising a RuntimeError), we fall back to a default # scale value to ensure the algorithm can continue and avoid infinite recursion. - if abs(scale) < 1e-3 or abs(scale)>10.0: - try: - scale, p0, p1 = estimate_goffset_scale(spgr, detector) - except RuntimeError: - logging.warning("Scale too small, using default 0.5") - scale = 0.5 + if abs(scale) < 1e-3 or abs(scale) > 10.0: + logging.warning( + "Unreliable scale estimate (%.4f). Retries left: %d", + scale, + retries + ) + + if retries > 0: + return estimate_goffset_scale(spgr, detector, delta, retries - 1) + + logging.warning("Scale estimation failed after retries, using default 0.5") + scale = 0.5 return scale, p0, p1 @@ -221,9 +226,7 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, max_step = 0.1 * (maxv - minv) # max 10% of range for i in range(max_it): - with future._task_lock: - if future._task_state == CANCELLED: - raise CancelledError() + _checkCancelled() if i == 0: peak_px = p0 @@ -247,8 +250,10 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, delta_goffset = max(-max_step, min(max_step, delta_goffset)) total_goffset_displacement += delta_goffset - logging.debug(f"DEBUG | Iter: {i} | Peak: {peak_px:.1f} | Error: {error_px:.1f} | Move: {delta_goffset:.4f} | Total Change: {total_goffset_displacement:.4f}") - spgr.moveRelSync({"goffset": delta_goffset}) + logging.debug( + "DEBUG | Iter: %d | Peak: %.1f | Error: %.1f | Move: %.4f | Total Change: %.4f", + i, peak_px, error_px, delta_goffset, total_goffset_displacement + ) spgr.moveRelSync({"goffset": delta_goffset}) future.set_progress( end=time.time() + (max_it - i - 1) * 0.5) # update estimated end time @@ -260,7 +265,7 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.debug("SparcAutoGratingOffset cancelled") raise except Exception as e: - logging.error(f"Alignment error: {e}") + logging.error("Alignment error: %s", e) raise def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): @@ -269,4 +274,194 @@ def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): Canceller of _do_sparc_auto_grating_offset task. """ with future._task_lock: - future._task_state = CANCELLED \ No newline at end of file + future._task_state = CANCELLED + + +def _checkCancelled(future: "model.ProgressiveFuture"): + + """ + Check if the future has been cancelled, and if so raise CancelledError. + """ + + with future._task_lock: + if future._task_state == CANCELLED: + raise CancelledError() + +def _total_alignment_time(n_gratings: int, + n_detectors: int) -> float: + + """ + Estimate total time for aligning all grating-detector combinations. + + :param n_gratings: number of gratings to align + :param n_detectors: number of detectors to align + :return: estimated total time in seconds + """ + + runs = n_detectors + max(0, n_gratings - 1) + move_time = ((n_gratings-1) * MOVE_TIME_GRATING + (n_detectors - 1) * MOVE_TIME_DETECTOR) + + # total time = time spent running alignment algorithms + time spent moving hardware + return runs * EST_ALIGN_TIME + move_time + + +def auto_align_grating_detector_offsets(spectrograph: model.Actuator, + detectors: Union[model.Detector, List[model.Detector]], + selector: Optional[model.Actuator] = None, + streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: + + """ + Automatically align grating-detector offsets for all combinations of gratings and detectors. + - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector + will be used for all subsequent gratings. + - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will + be aligned by adjusting the detector offset with the grating alignment fixed. + + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) + :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector + :raises CancelledError: if the operation is cancelled + """ + + if not isinstance(detectors, Iterable): + detectors = [detectors] + if not detectors: + raise ValueError("At least one detector must be provided") + if len(detectors) > 1 and selector is None: + raise ValueError("No selector provided, but multiple detectors") + + if streams is None: + streams = [] + + est_start = time.time() + 0.1 + n_gratings = len(spectrograph.axes["grating"].choices) + n_detectors = len(detectors) + a_time = _total_alignment_time(n_gratings, n_detectors) + f = model.ProgressiveFuture(start=est_start, end=est_start + a_time) + f.task_canceller = _cancel_auto_align_grating_detector_offsets + + f._task_lock = threading.Lock() + f._task_state = RUNNING + f._subfuture = InstantaneousFuture() + executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams)) + return f + +MOVE_TIME_GRATING = 20 #s +MOVE_TIME_DETECTOR = 5 #s +EST_ALIGN_TIME = 30 #s + +def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, + spectrograph: model.Actuator, + detectors: List[model.Detector], + selector: Optional[model.Actuator], + streams: List['Stream'], + stabilization_time: float = 15.0) -> Optional[Dict[Any, Any]]: + + """ + Iterate through each grating and detector combination, adjusting the selector if provided, and run the auto-alignment algorithm. + - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector + will be used for all subsequent gratings. + - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will + be aligned by adjusting the detector offset with the grating alignment fixed. + + :param future: ProgressiveFuture to update with progress and results + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :param stabilization_time: time to wait after moving hardware before starting alignment (default: 15s) + + :return: dict mapping (grating, detector) to alignment success boolean + :raises CancelledError: if the operation is cancelled + """ + + results: dict[tuple, bool] = {} + original_pos = {k: v for k, v in spectrograph.position.value.items() + if k in ("wavelength", "grating")} + + gratings = sorted(list(spectrograph.axes["grating"].choices.keys())) + logging.info(f"Available gratings: {list(spectrograph.axes['grating'].choices.keys())}") + + first_detector = detectors[0] + + if selector: + original_selector = selector.position.value + selector_axes, detector_to_selector = _mapDetectorToSelector(selector, detectors) + + def is_current_detector(d): + if selector is None: + return True + return detector_to_selector[d] == selector.position.value[selector_axes] + + try: + g0 = gratings[0] + logging.info("Starting alignment for initial grating: %s", g0) + + spectrograph.moveAbsSync({"grating": g0, "wavelength": 0}) + time.sleep(stabilization_time) + + detectors_sorted = sorted(detectors, key=is_current_detector, reverse=True) + + # align each detector for the first grating + for d in detectors_sorted: + _checkCancelled(future) + logging.info("Starting alignment | Detector: %s | Grating: %s", d.name, g0) + + if selector: + selector.moveAbsSync({selector_axes: detector_to_selector[d]}) + future._subfuture = sparc_auto_grating_offset(spectrograph, d) + success = future._subfuture.result() + results[(g0, d.name)] = success + + logging.info("Finished alignment | Detector: %s | Grating: %s | Success: %s", d.name, g0) + + if selector: + selector.moveAbsSync({selector_axes: detector_to_selector[first_detector]}) + + # align remaining gratings using the first detector + for g in gratings[1:]: + _checkCancelled(future) + logging.info("Switching to grating: %s", g) + + spectrograph.moveAbsSync({"grating": g, "wavelength": 0}) + time.sleep(stabilization_time) + logging.info("Starting alignment | Detector: %s | Grating: %s", first_detector.name, g) + + future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector) + success = future._subfuture.result() + results[(g, first_detector.name)] = success + + logging.info("Finished alignment | Detector: %s | Grating: %s", first_detector.name, g) + + return results + + except CancelledError: + logging.info("Auto-alignment cancelled") + raise + + finally: + spectrograph.moveAbsSync(original_pos) + if selector: + selector.moveAbsSync(original_selector) + + with future._task_lock: + future._task_state = FINISHED + +def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) -> bool: + + """ + Canceller for _do_auto_align_grating_detector_offsets task. + """ + logging.debug("Cancelling autoalignment...") + + with future._task_lock: + if future._task_state == FINISHED: + return False + future._task_state = CANCELLED + future._subfuture.cancel() + logging.debug("Auto-alignment cancellation requested") + + return True \ No newline at end of file diff --git a/src/odemis/acq/align/goffset_alignment.py b/src/odemis/acq/align/goffset_alignment.py index e48f07976e..431fc5dd48 100644 --- a/src/odemis/acq/align/goffset_alignment.py +++ b/src/odemis/acq/align/goffset_alignment.py @@ -15,22 +15,24 @@ def _checkCancelled(future: "model.ProgressiveFuture"): - lock_name = "_centering_lock" if hasattr(future, "_centering_lock") else "_align_lock" - state_name = "_centering_state" if hasattr(future, "_centering_state") else "_align_state" - lock = getattr(future, lock_name) - state = getattr(future, state_name) + """ + Check if the future has been cancelled, and if so raise CancelledError. + """ - with lock: - if state == CANCELLED: + with future._function_lock: + if future._function_state == CANCELLED: raise CancelledError() - def _total_alignment_time(n_gratings: int, n_detectors: int) -> float: """ Estimate total time for aligning all grating-detector combinations. + + :param n_gratings: number of gratings to align + :param n_detectors: number of detectors to align + :return: estimated total time in seconds """ runs = n_detectors + max(0, n_gratings - 1) @@ -45,6 +47,22 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, selector: Optional[model.Actuator] = None, streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: + """ + Automatically align grating-detector offsets for all combinations of gratings and detectors. + - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector + will be used for all subsequent gratings. + - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will + be aligned by adjusting the detector offset with the grating alignment fixed. + + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) + :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector + :raises CancelledError: if the operation is cancelled + """ + if not isinstance(detectors, Iterable): detectors = [detectors] if not detectors: @@ -62,8 +80,8 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, f = model.ProgressiveFuture(start=est_start, end=est_start + a_time) f.task_canceller = _cancel_auto_align_grating_detector_offsets - f._align_state = RUNNING - f._align_lock = threading.Lock() + f._task_lock = threading.Lock() + f._task_state = RUNNING f._subfuture = InstantaneousFuture() executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams)) return f @@ -81,8 +99,19 @@ def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, """ Iterate through each grating and detector combination, adjusting the selector if provided, and run the auto-alignment algorithm. - - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector will be used for all subsequent gratings. - - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will be aligned by adjusting the detector offset with the grating alignment fixed. + - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector + will be used for all subsequent gratings. + - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will + be aligned by adjusting the detector offset with the grating alignment fixed. + + :param future: ProgressiveFuture to update with progress and results + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :param stabilization_time: time to wait after moving hardware before starting alignment (default: 15s) + :return: dict mapping (grating, detector) to alignment success boolean + :raises CancelledError: if the operation is cancelled """ results: dict[tuple, bool] = {} @@ -154,17 +183,17 @@ def is_current_detector(d): if selector: selector.moveAbsSync(original_selector) - with future._align_lock: - future._align_state = FINISHED + with future._task_lock: + future._task_state = FINISHED def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) -> bool: logging.debug("Cancelling autoalignment...") - with future._align_lock: - if future._align_state == FINISHED: + with future._task_lock: + if future._task_state == FINISHED: return False - future._align_state = CANCELLED + future._task_state = CANCELLED future._subfuture.cancel() logging.debug("Auto-alignment cancellation requested") From 37c53b39ce4651f6a5da0ceb3f85669d4f6b0f2c Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Mon, 9 Mar 2026 15:22:47 +0100 Subject: [PATCH 06/11] [fix] fixed _checkcancelled and merged unittests goffset + goffset_alignment --- src/odemis/acq/align/goffset.py | 27 +++-- src/odemis/acq/align/goffset_alignment.py | 2 +- .../acq/align/test/goffset_alignment_test.py | 16 +-- src/odemis/acq/align/test/goffset_test.py | 100 ++++++++++++++++-- 4 files changed, 112 insertions(+), 33 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index 22a71e9072..b593c1700b 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -1,6 +1,7 @@ import logging import threading import time +from odemis.acq.stream import Stream from collections.abc import Iterable from concurrent.futures import CancelledError from concurrent.futures._base import CANCELLED, FINISHED, RUNNING @@ -91,7 +92,7 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0, retries = 1) -> Tuple[float, float, float]: """ - Estimate the scale factor between a change in the grating offset ('goffset') + Estimate the scale factor between a change in the grating offset ('goffset') and the resulting shift of the spectral peak on the detector. The function moves the actuator by a small step (delta) and measures the peak position before and after the move. @@ -106,7 +107,6 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). The actual step may be negated to avoid exceeding hardware limits. :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). - :return: Tuple (scale, p0, p1) scale: estimated pixels per unit of goffset p0: peak position at the initial goffset @@ -226,7 +226,7 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, max_step = 0.1 * (maxv - minv) # max 10% of range for i in range(max_it): - _checkCancelled() + _checkCancelled(future) if i == 0: peak_px = p0 @@ -253,7 +253,8 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.debug( "DEBUG | Iter: %d | Peak: %.1f | Error: %.1f | Move: %.4f | Total Change: %.4f", i, peak_px, error_px, delta_goffset, total_goffset_displacement - ) spgr.moveRelSync({"goffset": delta_goffset}) + ) + spgr.moveRelSync({"goffset": delta_goffset}) future.set_progress( end=time.time() + (max_it - i - 1) * 0.5) # update estimated end time @@ -268,16 +269,20 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.error("Alignment error: %s", e) raise -def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): +def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture) -> Bool: """ Canceller of _do_sparc_auto_grating_offset task. """ with future._task_lock: + if future._task_state == FINISHED: + return False future._task_state = CANCELLED + logging.debug("Cancelling alignment") + return True -def _checkCancelled(future: "model.ProgressiveFuture"): +def _checkCancelled(future: "model.ProgressiveFuture") -> None: """ Check if the future has been cancelled, and if so raise CancelledError. @@ -412,11 +417,11 @@ def is_current_detector(d): if selector: selector.moveAbsSync({selector_axes: detector_to_selector[d]}) - future._subfuture = sparc_auto_grating_offset(spectrograph, d) - success = future._subfuture.result() - results[(g0, d.name)] = success + future._subfuture = sparc_auto_grating_offset(spectrograph, d) + success = future._subfuture.result() + results[(g0, d.name)] = success - logging.info("Finished alignment | Detector: %s | Grating: %s | Success: %s", d.name, g0) + logging.info("Finished alignment | Detector: %s | Grating: %s", d.name, g0) if selector: selector.moveAbsSync({selector_axes: detector_to_selector[first_detector]}) @@ -464,4 +469,4 @@ def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) future._subfuture.cancel() logging.debug("Auto-alignment cancellation requested") - return True \ No newline at end of file + return True diff --git a/src/odemis/acq/align/goffset_alignment.py b/src/odemis/acq/align/goffset_alignment.py index 431fc5dd48..ee841cd0cc 100644 --- a/src/odemis/acq/align/goffset_alignment.py +++ b/src/odemis/acq/align/goffset_alignment.py @@ -197,4 +197,4 @@ def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) future._subfuture.cancel() logging.debug("Auto-alignment cancellation requested") - return True \ No newline at end of file + return True diff --git a/src/odemis/acq/align/test/goffset_alignment_test.py b/src/odemis/acq/align/test/goffset_alignment_test.py index 5957cfa158..cc5ad75caf 100644 --- a/src/odemis/acq/align/test/goffset_alignment_test.py +++ b/src/odemis/acq/align/test/goffset_alignment_test.py @@ -2,21 +2,12 @@ import os import time import unittest -import numpy -from collections.abc import Iterable from concurrent.futures import CancelledError -from scipy import ndimage -from odemis import model, acq -import odemis -from odemis.acq import align, stream, path -from odemis.acq.align.autofocus import Sparc2AutoFocus, MTD_BINARY +from odemis import model from odemis.acq.align.goffset_alignment import auto_align_grating_detector_offsets -from odemis.dataio import tiff, hdf5 -from odemis.util import testing, timeout, img +from odemis.util import timeout import odemis.util.focus -from unittest.mock import patch - CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" @@ -114,7 +105,6 @@ def test_multi_detector_iteration(self): # move to spectral camera self.selector.moveAbsSync({"rx": 1.5707963267948966}) data = spccd.data.get(asap=False) - # check data is not flat if data.max() == data.min(): print("WARNING: sp-ccd is returning a flat image!") @@ -122,4 +112,4 @@ def test_multi_detector_iteration(self): print(f"sp-ccd signal range: {data.min()} to {data.max()}") if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index 93a5fb99c2..8b69803c8e 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -8,17 +8,16 @@ import unittest import logging import numpy as np +import odemis from concurrent.futures._base import CancelledError - from odemis import model -from odemis.util import testing, timeout +from odemis.util import timeout from odemis.acq.align.goffset import (find_peak_position, estimate_goffset_scale, - sparc_auto_grating_offset) - -import odemis + sparc_auto_grating_offset, + auto_align_grating_detector_offsets) logging.getLogger().setLevel(logging.DEBUG) @@ -59,7 +58,7 @@ def test_find_peak_position_synthetic(self): """ x = np.arange(200) true_center = 83.4 - spectrum = np.exp(-0.5*((x-true_center)/3.0)**2) + spectrum = np.exp(-0.5 * ((x - true_center) / 3.0)**2) peak = find_peak_position(spectrum) self.assertAlmostEqual(peak, true_center, places=1) @@ -70,7 +69,7 @@ def test_find_peak_position_2d(self): """ x = np.arange(200) true_center = 120.0 - line = np.exp(-0.5*((x-true_center)/4.0)**2) + line = np.exp(-0.5 * ((x - true_center) / 4.0)**2) image = np.tile(line, (50, 1)) peak = find_peak_position(image) @@ -116,5 +115,90 @@ def test_cancel(self): pass self.assertTrue(f.done()) +class TestAutoAlignGratingDetectorOffsets(unittest.TestCase): + + @classmethod + def setUpClass(cls): + cls.spgr = model.getComponent(role="spectrograph") + cls.ccd = model.getComponent(role="ccd") + #cls.spccd = model.getComponent(role="sp-ccd") + cls.selector = model.getComponent(role="spec-det-selector") + + def setUp(self): + # Speed up acquisition + self.ccd.exposureTime.value = self.ccd.exposureTime.range[0] + + @timeout(1000) + def test_single_detector_iteration(self): + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + res = f.result(timeout=900) + + n_gratings = len(self.spgr.axes["grating"].choices) + n_detectors = 1 + expected = n_detectors + (n_gratings - 1) + + self.assertEqual(len(res), expected) + + first_grating = list(self.spgr.axes["grating"].choices.keys())[0] + dets_first = [d for (g, d) in res.keys() if g == first_grating] + + self.assertEqual(len(dets_first), n_detectors) + + def test_multi_detector_iteration(self): + spccd = model.getComponent(role="sp-ccd") + spccd.exposureTime.value = spccd.exposureTime.range[0] + + detectors = [self.ccd, spccd] + + # run alignment + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) + res = f.result(timeout=900) + + # calculate expected results + n_gratings = len(self.spgr.axes["grating"].choices) + n_detectors = len(detectors) + expected_count = n_detectors + (n_gratings - 1) + + self.assertEqual(len(res), expected_count, f"Expected {expected_count} results, got {len(res)}") + + # verify that every detector was used for the first grating + gratings_list = list(self.spgr.axes["grating"].choices.keys()) + first_grating = gratings_list[0] + + dets_for_first_grating = [d for (g, d) in res.keys() if g == first_grating] + self.assertEqual(len(dets_for_first_grating), n_detectors) + self.assertIn(self.ccd.name, dets_for_first_grating) + self.assertIn(spccd.name, dets_for_first_grating) + + # verify that only first detector is used for remaining gratings + for g in gratings_list[1:]: + # Ensure only the first detector (index 0) is present for these gratings + dets_for_this_grating = [d for (grating, d) in res.keys() if grating == g] + self.assertEqual(len(dets_for_this_grating), 1) + self.assertEqual(dets_for_this_grating[0], detectors[0].name) + + # move to spectral camera + self.selector.moveAbsSync({"rx": 1.5707963267948966}) + data = spccd.data.get(asap=False) + + # check data is not flat + if data.max() == data.min(): + print("WARNING: sp-ccd is returning a flat image!") + else: + print(f"sp-ccd signal range: {data.min()} to {data.max()}") + + @timeout(100) + def test_cancel(self): + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd],) + + time.sleep(1) + + cancelled = f.cancel() + self.assertTrue(cancelled) + self.assertTrue(f.cancelled()) + + with self.assertRaises(CancelledError): + f.result(timeout=900) + if __name__ == "__main__": - unittest.main() \ No newline at end of file + unittest.main() From 75e74351f774af2ab0baeb64a138195039849771 Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Tue, 17 Mar 2026 09:46:50 +0100 Subject: [PATCH 07/11] [feat] added cornercase --- src/odemis/acq/align/goffset.py | 375 +++++++++++++----- .../acq/align/test/goffset_alignment_test.py | 44 +- src/odemis/acq/align/test/goffset_test.py | 250 +++++++++--- 3 files changed, 508 insertions(+), 161 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index b593c1700b..047dd7fd78 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -1,7 +1,6 @@ import logging import threading import time -from odemis.acq.stream import Stream from collections.abc import Iterable from concurrent.futures import CancelledError from concurrent.futures._base import CANCELLED, FINISHED, RUNNING @@ -15,17 +14,19 @@ from odemis.util import executeAsyncTask from scipy.optimize import curve_fit -def gaussian(x, amplitude, x0, width): + +def gaussian(x, amplitude, x0, width) -> numpy.ndarray: """ Gaussian function (for curve fitting). - :param x: input coordinates - :param amplitude: peak intensity - :param x0 = peak's center - :param width = standard deviation - :return: Gaussian function evaluated at x + :param x: input coordinates + :param amplitude: peak intensity + :param x0 = peak's center + :param width = standard deviation + :return: Gaussian function evaluated at x (numpy array) """ - intensity = amplitude*numpy.exp(-0.5 * ((x - x0) / width)**2) + + intensity = amplitude * numpy.exp(-0.5 * ((x - x0) / width) ** 2) return intensity @@ -44,6 +45,11 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: average as the peak position estimate. The weighted average is more robust as it does not run an iterative optimizer, so it has no convergence or numerical-optimization failure modes, but it may be less accurate if the peak is not well-defined or if there are multiple peaks within the window. + + :param data: 1D or 2D array containing the spectrum (if 2D, it will be averaged to 1D) + :param window_radius: number of pixels on either side of the peak to include in the window for fitting (default: 15) + :return: estimated peak position in pixels (float) + :raises RuntimeError: if no significant peak is detected (SNR too low) """ if data.ndim == 2: @@ -51,6 +57,17 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: else: spectrum = data + # compute simple SNR by comparing the peak height above the median background to the background magnitude and returns + # False if that SNR is below the threshold. This helps reject cases with no real peak, + # where the Gaussian fit would fail or produce nonsense results. + + peak_value = float(spectrum.max()) + background = float(numpy.median(spectrum)) + snr = (peak_value - background) / (abs(background) + 1e-6) + + if snr < 1: # tune threshold + raise RuntimeError("No peak detected (SNR too low)") + peak_idx = numpy.argmax(spectrum) # find the absolute highest point # create a window around the peak @@ -67,14 +84,14 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: if weights.sum() == 0: weighted_avg = float(peak_idx) logging.info("Weighted average fallback: all window data <= 0, using peak_idx=%d as estimate", - peak_idx) + peak_idx) else: weighted_avg = float(numpy.sum(window_idx * window_data) / numpy.sum(window_data)) # try Gaussian fit for better accuracy, but fallback to weighted average if it fails or is out of bounds try: - p0 = [window_data.max(), weighted_avg, 2.5] # intial guess: [amplitude, center, width] + p0 = [window_data.max(), weighted_avg, 2.5] # intial guess: [amplitude, center, width] popt, pcov = curve_fit(gaussian, window_idx, window_data, p0=p0) peak = popt[1] @@ -89,8 +106,113 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: return weighted_avg +def peak_is_present(spectrum, snr_threshold=1, width_range=(0.5, 12.0)) -> bool: + """ + Test to decide whether spectral peak is present. + + The test uses: + - a simple SNR threshold comparing the maximum to the median background, + - a local width estimate computed from a small window around the peak to + reject hot pixels and extremely broad features. -def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0, retries = 1) -> Tuple[float, float, float]: + :param spectrum: 1D array containing the spectrum (intensity vs pixel) + :param snr_threshold: minimum required signal-to-noise ratio for a peak to be considered present (default: 1) + :param width_range: acceptable range of estimated peak widths in pixels (default: (0.5, 12.0)) + :return: True if a peak meeting the criteria is present, False otherwise + """ + + # basic stats + peak_value = float(spectrum.max()) + background = float(numpy.median(spectrum)) + snr = (peak_value - background) / (abs(background) + 1e-6) + + if snr < snr_threshold: + return False + + # estimate width around the peak + peak_idx = int(numpy.argmax(spectrum)) + if peak_idx < 1 or peak_idx > len(spectrum) - 2: + logging.debug("Peak too close to edge idx=%d len=%d", peak_idx, len(spectrum)) + return False + + window = spectrum[peak_idx-2 : peak_idx+3] + x = numpy.arange(len(window)) + w = window - window.min() + if w.sum() == 0: + return False + + mean = numpy.sum(x * w) / numpy.sum(w) + var = numpy.sum(w * (x - mean)**2) / numpy.sum(w) + width = numpy.sqrt(var) + present = width_range[0] <= width <= width_range[1] + logging.debug("snr=%.2f width=%.2f present=%s", snr, width, present) + + return present + + +def acquire_peak(spgr, detector, step=2000) -> Tuple[float, float]: + """ + Coarse scan across the goffset axis until a real peak becomes visible. + + The function performs absolute moves across the configured goffset axis + from min to max in steps of `step`. After each move it reads the detector, + computes a 1D spectrum (mean across the first axis), and applies + `peak_is_present` to decide whether a real peak is on the detector. + When a peak is found, `find_peak_position` is used to estimate its pixel + position and the function returns the actual goffset (as reported by the + actuator) and the peak pixel. + + :param spgr: spectrograph + :param detector: detector + :param step: step size in goffset units for the coarse scan (default: 2000) + :return: tuple (actual_goffset, peak_pixel) + :raises RuntimeError: if no peak is found across the full goffset range + """ + + current = float(spgr.position.value["goffset"]) + step = abs(step) if step != 0 else 2000.0 + max_span = 20000.0 # limit how far we wander from the current valid position + + logging.debug( + "Coarse local scan around goffset %.1f with step %.1f and max span %.1f", + current, step, max_span + ) + + # Build a sequence of positions: current, +step, -step, +2*step, -2*step, ... + positions = [current] + k = 1 + while k * step <= max_span: + positions.append(current + k * step) + positions.append(current - k * step) + k += 1 + + tried = set() + + for g in positions: + # avoid duplicate moves due to symmetry or float rounding + key = round(g, 3) + if key in tried: + continue + tried.add(key) + + logging.debug("Coarse scan move attempt to goffset %.3f", g) + try: + spgr.moveAbsSync({"goffset": g}) + except Exception as e: + logging.warning("Skipping invalid goffset %.3f (%s)", g, e) + continue + + data = detector.data.get(asap=False) + spectrum = data.mean(axis=0) + + if peak_is_present(spectrum): + peak = find_peak_position(data) + return g, peak + + raise RuntimeError("Peak not found in local goffset scan around current position") + +def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0, retries=1) -> Tuple[ + float, float, float]: """ Estimate the scale factor between a change in the grating offset ('goffset') and the resulting shift of the spectral peak on the detector. @@ -102,15 +224,16 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta If the measured scale is unreasonably small or large, the function retries recursively and falls back to a default value of 0.5 if necessary. - :param spgr: spectrograph - :param detector: detector - :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). - The actual step may be negated to avoid exceeding hardware limits. - :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). - :return: Tuple (scale, p0, p1) - scale: estimated pixels per unit of goffset - p0: peak position at the initial goffset - p1: peak position after applying the test delta + :param spgr: spectrograph + :param detector: detector + :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). + The actual step may be negated to avoid exceeding hardware limits. + :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). + + :return: Tuple (scale, p0, p1) + scale: estimated pixels per unit of goffset + p0: peak position at the initial goffset + p1: peak position after applying the test delta """ # get initial state @@ -123,7 +246,7 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta # ensure that max limit isn't violated move_direction = 1 if (current_pos + delta < goffset_max) else -1 - actual_delta = delta*move_direction + actual_delta = delta * move_direction # move and measure spgr.moveRelSync({"goffset": actual_delta}) @@ -162,27 +285,29 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta return scale, p0, p1 + def sparc_auto_grating_offset(spgr: model.Actuator, - detector: model.Detector, - tolerance_px: float = 0.4, - max_it: int = 20, - gain: float = 0.4) -> model.ProgressiveFuture: + detector: model.Detector, + single_detector_mode: bool = False, + tolerance_px: float = 0.4, + max_it: int = 20, + gain: float = 0.4) -> model.ProgressiveFuture: """ Start an asynchronous task that centers the spectral peak by adjusting the grating offset (goffset). - :param spgr: spectrograph - :param detector: detector - :param tolerance_px: the acceptable displacement of the peak from the center in pixels (default: 0.4) - :param max_it: maximum number of iterations to attempt (default: 20) - :param gain: proportional gain factor for adjusting the goffset (default: 0.4) - :return: A ``ProgressiveFuture`` representing the asynchronous alignment - task. The future can be used to monitor progress, retrieve the - result, or cancel the alignment. + :param spgr: spectrograph + :param detector: detector + :param tolerance_px: the acceptable displacement of the peak from the center in pixels (default: 0.4) + :param max_it: maximum number of iterations to attempt (default: 20) + :param gain: proportional gain factor for adjusting the goffset (default: 0.4) + :return: A ``ProgressiveFuture`` representing the asynchronous alignment + task. The future can be used to monitor progress, retrieve the + result, or cancel the alignment. """ est_start = time.time() + 0.05 - est_time = max_it*0.5 # conservative estimate + est_time = max_it * 0.5 # conservative estimate f = model.ProgressiveFuture(start=est_start, end=est_start + est_time) @@ -193,73 +318,125 @@ def sparc_auto_grating_offset(spgr: model.Actuator, executeAsyncTask( f, _do_sparc_auto_grating_offset, - args=(f, spgr, detector, tolerance_px, max_it, gain), + args=(f, spgr, detector, single_detector_mode, tolerance_px, max_it, gain), ) return f def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, - spgr: model.Actuator, - detector: model.Detector, - tolerance_px: float, - max_it: int, - gain: float) -> bool: + spgr: model.Actuator, + detector: model.Detector, + single_detector_mode, + tolerance_px: float, + max_it: int, + gain: float) -> bool: """ - Iteratively adjusts the grating offset to align the spectral peak to the center of the detector. - The algorithm estimates the current peak position, calculates the error from the center, and moves the grating offset - proportionally to reduce this error. - It continues until the peak is within the specified tolerance or the maximum number of iterations is reached. + Core alignment routine that iteratively adjusts the grating offset to center the peak. + + Behavior: + - Attempts to acquire a peak on the detector (coarse scan) if none is present. + - Measures the local goffset-to-pixel scale only when a peak is present and not centered. + - Runs a centering loop until the peak is within `tolerance_px` or the max #iterations is reached. + - Respects cancellation via the provided ProgressiveFuture. + + :param future: model.ProgressiveFuture + :param spgr: spectrograph + :param detector: detector + :param tolerance_px: pixel tolerance for successful alignment + :param max_it: maximum number of centering iterations + :param gain: proportional gain for converting pixel error to goffset correction + :return: True if alignment succeeded (peak within tolerance), False otherwise (bool) + :raises CancelledError: if the future was cancelled during execution """ - logging.info("Running alignment | detector=%s |", - detector.name) + logging.info("Running alignment | detector=%s |", detector.name) try: + center_target = detector.resolution.value[0] / 2 + + # initial read: try to get a valid peak without moving the grating + try: + data0 = detector.data.get(asap=False) + peak0 = find_peak_position(data0) # raises RuntimeError if no peak present + logging.debug("Initial read: peak0=%.2f px", peak0) + peak_present = True + except RuntimeError: + logging.debug("Initial read: no significant peak detected, scanning the axes") + peak_present = False + peak0 = None + + # if peak present and already centered, do nothing + if peak_present: + initial_error_px = peak0 - center_target + logging.debug("Initial error_px=%.3f px (tolerance=%.3f)", initial_error_px, tolerance_px) + if abs(initial_error_px) <= tolerance_px: + logging.info("Peak already centered | peak=%.2f | center=%.2f | error=%.3f", + peak0, center_target, initial_error_px) + return True + + # if no peak present, run acquisition (this will move the grating) + if not peak_present: + try: + g_acq, p_acq = acquire_peak(spgr, detector, step=2000) + logging.info("Peak acquired at goffset=%d pixel=%.2f", g_acq, p_acq) + peak0 = p_acq + peak_present = True + except RuntimeError: + logging.error("Peak acquisition failed — aborting alignment") + return False + + # after acquisition, check if acquisition already placed peak within tolerance + post_acq_error_px = peak0 - center_target + logging.debug("Post-acquisition error_px=%.3f px", post_acq_error_px) + if abs(post_acq_error_px) <= tolerance_px: + logging.info("Peak centered by acquisition | peak=%.2f | center=%.2f | error=%.3f", + peak0, center_target, post_acq_error_px) + return True + + # peak is present and not centered -> estimate scale and center scale, p0, p1 = estimate_goffset_scale(spgr, detector) - center_target = detector.resolution.value[0] / 2 # adjust if 0 is not the center - total_goffset_displacement = 0.0 + logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) + + # prefer p1 (after probe move) if available, else p0 + start_peak = p1 if p1 is not None else p0 - # clamp move to safe fraction of axis range + total_goffset_displacement = 0.0 axis = spgr.axes["goffset"] minv, maxv = axis.range - max_step = 0.1 * (maxv - minv) # max 10% of range for i in range(max_it): _checkCancelled(future) if i == 0: - peak_px = p0 + peak_px = start_peak else: data = detector.data.get(asap=False) - peak_px = find_peak_position(data) + try: + peak_px = find_peak_position(data) + except RuntimeError: + logging.error("No peak detected during centering iteration %d — aborting", i) + return False error_px = peak_px - center_target if abs(error_px) <= tolerance_px: - logging.info( - "Spectral peak aligned after %d iterations | error_px=%.2f | total_goffset_change=%.4f", - i + 1, - error_px, - total_goffset_displacement - ) return True - delta_goffset = -gain*(error_px/scale) - - delta_goffset = max(-max_step, min(max_step, delta_goffset)) + delta_goffset = -gain * (error_px / scale) + current = spgr.position.value["goffset"] + delta_goffset = max(minv - current, min(maxv - current, delta_goffset)) total_goffset_displacement += delta_goffset logging.debug( - "DEBUG | Iter: %d | Peak: %.1f | Error: %.1f | Move: %.4f | Total Change: %.4f", + "Iter: %d | Peak: %.2f | Error: %.2f | Move: %.6f | Total Change: %.6f", i, peak_px, error_px, delta_goffset, total_goffset_displacement ) - spgr.moveRelSync({"goffset": delta_goffset}) - future.set_progress( - end=time.time() + (max_it - i - 1) * 0.5) # update estimated end time + spgr.moveRelSync({"goffset": delta_goffset}) + future.set_progress(end=time.time() + (max_it - i - 1) * 0.5) - logging.warning("SparcAutoGratingOffset did not converge") + logging.warning("SparcAutoGratingOffset did not converge within max iterations") return False except CancelledError: @@ -269,21 +446,16 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.error("Alignment error: %s", e) raise -def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture) -> Bool: +def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): """ Canceller of _do_sparc_auto_grating_offset task. """ with future._task_lock: - if future._task_state == FINISHED: - return False future._task_state = CANCELLED - logging.debug("Cancelling alignment") - - return True -def _checkCancelled(future: "model.ProgressiveFuture") -> None: +def _checkCancelled(future: "model.ProgressiveFuture"): """ Check if the future has been cancelled, and if so raise CancelledError. """ @@ -292,9 +464,9 @@ def _checkCancelled(future: "model.ProgressiveFuture") -> None: if future._task_state == CANCELLED: raise CancelledError() -def _total_alignment_time(n_gratings: int, - n_detectors: int) -> float: +def _total_alignment_time(n_gratings: int, + n_detectors: int) -> float: """ Estimate total time for aligning all grating-detector combinations. @@ -304,17 +476,16 @@ def _total_alignment_time(n_gratings: int, """ runs = n_detectors + max(0, n_gratings - 1) - move_time = ((n_gratings-1) * MOVE_TIME_GRATING + (n_detectors - 1) * MOVE_TIME_DETECTOR) + move_time = ((n_gratings - 1) * MOVE_TIME_GRATING + (n_detectors - 1) * MOVE_TIME_DETECTOR) # total time = time spent running alignment algorithms + time spent moving hardware return runs * EST_ALIGN_TIME + move_time def auto_align_grating_detector_offsets(spectrograph: model.Actuator, - detectors: Union[model.Detector, List[model.Detector]], - selector: Optional[model.Actuator] = None, - streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: - + detectors: Union[model.Detector, List[model.Detector]], + selector: Optional[model.Actuator] = None, + streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: """ Automatically align grating-detector offsets for all combinations of gratings and detectors. - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector @@ -327,8 +498,8 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, :param selector: optional selector to switch between detectors :param streams: optional list of streams to update with progress :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) - :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector - :raises CancelledError: if the operation is cancelled + :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector + :raises CancelledError: if the operation is cancelled """ if not isinstance(detectors, Iterable): @@ -341,6 +512,8 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, if streams is None: streams = [] + single_detector_mode = len(detectors) == 1 + est_start = time.time() + 0.1 n_gratings = len(spectrograph.axes["grating"].choices) n_detectors = len(detectors) @@ -351,20 +524,22 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, f._task_lock = threading.Lock() f._task_state = RUNNING f._subfuture = InstantaneousFuture() - executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams)) + executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams, single_detector_mode)) return f -MOVE_TIME_GRATING = 20 #s -MOVE_TIME_DETECTOR = 5 #s -EST_ALIGN_TIME = 30 #s -def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, - spectrograph: model.Actuator, - detectors: List[model.Detector], - selector: Optional[model.Actuator], - streams: List['Stream'], - stabilization_time: float = 15.0) -> Optional[Dict[Any, Any]]: +MOVE_TIME_GRATING = 20 # s +MOVE_TIME_DETECTOR = 5 # s +EST_ALIGN_TIME = 30 # s + +def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, + spectrograph: model.Actuator, + detectors: List[model.Detector], + selector: Optional[model.Actuator], + streams: List['Stream'], + single_detector_mode: bool = False, + stabilization_time: float = 10.0) -> Optional[Dict[Any, Any]]: """ Iterate through each grating and detector combination, adjusting the selector if provided, and run the auto-alignment algorithm. - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector @@ -379,8 +554,8 @@ def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, :param streams: optional list of streams to update with progress :param stabilization_time: time to wait after moving hardware before starting alignment (default: 15s) - :return: dict mapping (grating, detector) to alignment success boolean - :raises CancelledError: if the operation is cancelled + :return: dict mapping (grating, detector) to alignment success boolean + :raises CancelledError: if the operation is cancelled """ results: dict[tuple, bool] = {} @@ -417,11 +592,11 @@ def is_current_detector(d): if selector: selector.moveAbsSync({selector_axes: detector_to_selector[d]}) - future._subfuture = sparc_auto_grating_offset(spectrograph, d) - success = future._subfuture.result() - results[(g0, d.name)] = success + future._subfuture = sparc_auto_grating_offset(spectrograph, d, single_detector_mode=single_detector_mode) + success = future._subfuture.result() + results[(g0, d.name)] = success - logging.info("Finished alignment | Detector: %s | Grating: %s", d.name, g0) + logging.info("Finished alignment | Detector: %s | Grating: %s", d.name, g0) if selector: selector.moveAbsSync({selector_axes: detector_to_selector[first_detector]}) @@ -435,7 +610,7 @@ def is_current_detector(d): time.sleep(stabilization_time) logging.info("Starting alignment | Detector: %s | Grating: %s", first_detector.name, g) - future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector) + future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector, single_detector_mode=single_detector_mode) success = future._subfuture.result() results[(g, first_detector.name)] = success @@ -455,8 +630,8 @@ def is_current_detector(d): with future._task_lock: future._task_state = FINISHED -def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) -> bool: +def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) -> bool: """ Canceller for _do_auto_align_grating_detector_offsets task. """ diff --git a/src/odemis/acq/align/test/goffset_alignment_test.py b/src/odemis/acq/align/test/goffset_alignment_test.py index cc5ad75caf..cf32fc3233 100644 --- a/src/odemis/acq/align/test/goffset_alignment_test.py +++ b/src/odemis/acq/align/test/goffset_alignment_test.py @@ -2,13 +2,28 @@ import os import time import unittest +import numpy +from collections.abc import Iterable from concurrent.futures import CancelledError +from scipy import ndimage -from odemis import model -from odemis.acq.align.goffset_alignment import auto_align_grating_detector_offsets -from odemis.util import timeout +from odemis import model, acq + +#from odemis.acq.align.goffset import auto_align_grating_detector_offsets + +from odemis.util import testing, timeout, img import odemis.util.focus +from odemis.acq.align.goffset import( + find_peak_position, + acquire_peak, + estimate_goffset_scale, + sparc_auto_grating_offset, + auto_align_grating_detector_offsets +) + + + CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" @@ -27,17 +42,17 @@ def setUp(self): # Speed up acquisition self.ccd.exposureTime.value = self.ccd.exposureTime.range[0] - @timeout(100) - def test_cancel(self): - f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd],) - time.sleep(1) - - cancelled = f.cancel() - self.assertTrue(cancelled) - self.assertTrue(f.cancelled()) - - with self.assertRaises(CancelledError): - f.result(timeout=900) + # @timeout(100) + # def test_cancel(self): + # f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd],) + # time.sleep(1) + # + # cancelled = f.cancel() + # self.assertTrue(cancelled) + # self.assertTrue(f.cancelled()) + # + # with self.assertRaises(CancelledError): + # f.result(timeout=900) @timeout(1000) @@ -105,6 +120,7 @@ def test_multi_detector_iteration(self): # move to spectral camera self.selector.moveAbsSync({"rx": 1.5707963267948966}) data = spccd.data.get(asap=False) + # check data is not flat if data.max() == data.min(): print("WARNING: sp-ccd is returning a flat image!") diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index 8b69803c8e..a01a72d9b8 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -4,27 +4,25 @@ """ import os -import time import unittest import logging import numpy as np import odemis -from concurrent.futures._base import CancelledError from odemis import model from odemis.util import timeout - -from odemis.acq.align.goffset import (find_peak_position, +from odemis.acq.align.goffset import( + find_peak_position, estimate_goffset_scale, sparc_auto_grating_offset, - auto_align_grating_detector_offsets) + auto_align_grating_detector_offsets +) logging.getLogger().setLevel(logging.DEBUG) CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" - class TestSparcAutoGratingOffset(unittest.TestCase): """ Test automatic grating offset alignment. @@ -36,6 +34,8 @@ def setUpClass(cls): cls.detector = model.getComponent(role="ccd") cls.spgr = model.getComponent(role="spectrograph") + # cls.spccd = model.getComponent(role="sp-ccd") + cls.selector = model.getComponent(role="spec-det-selector") cls._original_goffset = cls.spgr.goffset cls._original_position = cls.spgr.position.value.copy() @@ -58,7 +58,7 @@ def test_find_peak_position_synthetic(self): """ x = np.arange(200) true_center = 83.4 - spectrum = np.exp(-0.5 * ((x - true_center) / 3.0)**2) + spectrum = np.exp(-0.5*((x-true_center)/3.0)**2) peak = find_peak_position(spectrum) self.assertAlmostEqual(peak, true_center, places=1) @@ -69,7 +69,7 @@ def test_find_peak_position_2d(self): """ x = np.arange(200) true_center = 120.0 - line = np.exp(-0.5 * ((x - true_center) / 4.0)**2) + line = np.exp(-0.5*((x-true_center)/4.0)**2) image = np.tile(line, (50, 1)) peak = find_peak_position(image) @@ -86,51 +86,83 @@ def test_estimate_goffset_scale(self): self.assertNotEqual(scale, 0.0) self.assertTrue(np.isfinite(scale)) - @timeout(300) - def test_auto_grating_offset(self): + @timeout(800) + def test_scale_not_misaligned(self): """ - Test automatic centering of spectral peak. + Verify scale estimation only happens when the peak is misaligned. + This is inferred from the probe move performed by estimate_goffset_scale(). """ - delta = 0 # intentionally misalign - current = self.spgr.position.value["goffset"] - goffset_max = self.spgr.axes["goffset"].range[1] - direction = 1 if (current + delta < goffset_max) else -1 + # reset spectrograph to known aligned position + self.spgr.moveAbsSync(self._original_position) + + start_goffset = self.spgr.position.value["goffset"] - self.spgr.moveRelSync({"goffset": delta * direction}) f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=20) + result = f.result(timeout=300) + + end_goffset = self.spgr.position.value["goffset"] - result = f.result(timeout=200) self.assertTrue(result) - @timeout(100) - def test_cancel(self): + # If peak is already centered, the algorithm exits immediately + # so goffset should not change. + self.assertAlmostEqual( + start_goffset, + end_goffset, + places=6, + msg="goffset changed even though peak was already centered (scale estimation likely ran)" + ) + + def test_scale_estimation_misaligned(self): """ - Test cancelling alignment. + Verify that scale estimation runs when the peak is misaligned, which is inferred + from the probe move performed by estimate_goffset_scale(). """ - f = sparc_auto_grating_offset(self.spgr, self.detector) - # Wait for the result or a timeout - try: - f.result(timeout=5) - except: - pass - self.assertTrue(f.done()) -class TestAutoAlignGratingDetectorOffsets(unittest.TestCase): + delta = 500 + current = self.spgr.position.value["goffset"] + maxv = self.spgr.axes["goffset"].range[1] + direction = 1 if (current + delta < maxv) else -1 - @classmethod - def setUpClass(cls): - cls.spgr = model.getComponent(role="spectrograph") - cls.ccd = model.getComponent(role="ccd") - #cls.spccd = model.getComponent(role="sp-ccd") - cls.selector = model.getComponent(role="spec-det-selector") + self.spgr.moveRelSync({"goffset": delta * direction}) - def setUp(self): - # Speed up acquisition - self.ccd.exposureTime.value = self.ccd.exposureTime.range[0] + start_goffset = self.spgr.position.value["goffset"] + + f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=50) + result = f.result(timeout=600) + + end_goffset = self.spgr.position.value["goffset"] + + self.assertTrue(result) + + # If misaligned, centering should move the grating + self.assertNotAlmostEqual( + start_goffset, + end_goffset, + places=3, + msg="goffset did not change during alignment when peak was misaligned" + ) + + @timeout(800) + def test_auto_grating_offset(self): + """ + Test automatic centering of spectral peak. + """ + delta = 0 # intentionally misalign + current = self.spgr.position.value["goffset"] + goffset_max = self.spgr.axes["goffset"].range[1] + direction = 1 if (current + delta < goffset_max) else -1 + + self.spgr.moveRelSync({"goffset": delta * direction}) + logging.info("Test: after misalign move, spgr.position.gooffset = %s", self.spgr.position.value["goffset"]) + f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=50) + + result = f.result(timeout=800) + self.assertTrue(result) @timeout(1000) def test_single_detector_iteration(self): - f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.detector], selector=self.selector) res = f.result(timeout=900) n_gratings = len(self.spgr.axes["grating"].choices) @@ -187,18 +219,142 @@ def test_multi_detector_iteration(self): else: print(f"sp-ccd signal range: {data.min()} to {data.max()}") - @timeout(100) - def test_cancel(self): - f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd],) + @timeout(600) + def test_single_detector_updates_grating(self): + """ + Single-detector mode: aligning a detector on the secondary port must update + the *grating* goffset (not a detector-specific offset). + """ + spccd = model.getComponent(role="sp-ccd") + spccd.exposureTime.value = spccd.exposureTime.range[0] + + # try to move selector to secondary port; ignore if not available + try: + self.selector.moveAbsSync({"rx": 1.5707963267948966}) + except Exception: + logging.debug("Selector move to secondary failed or not present; continuing") + + start_goffset = self.spgr.position.value["goffset"] + + f = sparc_auto_grating_offset(self.spgr, spccd, single_detector_mode=True, max_it=50) + result = f.result(timeout=300) + self.assertTrue(result, "Single-detector alignment failed") + + after_goffset = self.spgr.position.value["goffset"] + self.assertNotAlmostEqual( + start_goffset, + after_goffset, + places=3, + msg="Grating goffset did not change for single-detector alignment on secondary port" + ) + + @timeout(900) + def test_multi_detector_does_not_change_grating(self): + """ + Multi-detector mode: first detector sets the grating; aligning the second + detector should not change the grating goffset (it should adjust detector offset). + """ + spccd = model.getComponent(role="sp-ccd") + spccd.exposureTime.value = spccd.exposureTime.range[0] + + # ensure selector points to primary detector for initial grating-setting alignment + try: + self.selector.moveAbsSync({"rx": 0.0}) + except Exception: + logging.debug("Selector move to primary failed or not present; continuing") + + # align first detector in single-detector mode to set grating + f_first = sparc_auto_grating_offset(self.spgr, self.detector, single_detector_mode=True, max_it=50) + self.assertTrue(f_first.result(timeout=300), "First-detector alignment failed") + grating_after_first = self.spgr.position.value["goffset"] - time.sleep(1) + # switch to secondary detector + try: + self.selector.moveAbsSync({"rx": 1.5707963267948966}) + except Exception: + logging.debug("Selector move to secondary failed or not present; continuing") + + # align second detector in multi-detector mode (should not change grating) + f_second = sparc_auto_grating_offset(self.spgr, spccd, single_detector_mode=False, max_it=50) + self.assertTrue(f_second.result(timeout=300), "Second-detector alignment failed") + + grating_after_second = self.spgr.position.value["goffset"] + self.assertAlmostEqual( + grating_after_first, + grating_after_second, + places=3, + msg="Grating goffset changed when aligning second detector in multi-detector mode" + ) + + @timeout(900) + def test_multi_detector_detector_offset_changes(self): + """ + Multi-detector mode: after the first detector sets the grating, aligning the + second detector should NOT change the grating goffset but should change the + detector-specific response (verified by the peak moving closer to center). + """ + spccd = model.getComponent(role="sp-ccd") + spccd.exposureTime.value = spccd.exposureTime.range[0] - cancelled = f.cancel() - self.assertTrue(cancelled) - self.assertTrue(f.cancelled()) + # Ensure primary detector sets the grating first + try: + self.selector.moveAbsSync({"rx": 0.0}) + except Exception: + logging.debug("Selector move to primary failed or not present; continuing") - with self.assertRaises(CancelledError): - f.result(timeout=900) + f_first = sparc_auto_grating_offset(self.spgr, self.detector, single_detector_mode=True, max_it=50) + self.assertTrue(f_first.result(timeout=300), "First-detector alignment failed") + grating_after_first = self.spgr.position.value["goffset"] + + # Switch to secondary detector + try: + self.selector.moveAbsSync({"rx": 1.5707963267948966}) + except Exception: + logging.debug("Selector move to secondary failed or not present; continuing") + + # Measure peak before alignment on secondary detector + data_before = spccd.data.get(asap=False) + before_peak = float(find_peak_position(data_before)) + + # Run alignment for second detector in multi-detector mode + f_second = sparc_auto_grating_offset(self.spgr, spccd, single_detector_mode=False, max_it=50) + self.assertTrue(f_second.result(timeout=600), "Second-detector alignment failed") + + # Measure peak after alignment + data_after = spccd.data.get(asap=False) + after_peak = float(find_peak_position(data_after)) + + # Assert grating did not change + grating_after_second = self.spgr.position.value["goffset"] + self.assertAlmostEqual( + grating_after_first, + grating_after_second, + places=3, + msg="Grating goffset changed when aligning second detector in multi-detector mode" + ) + + # Assert peak moved closer to center on the secondary detector + center = spccd.resolution.value[0] / 2 + before_dist = abs(before_peak - center) + after_dist = abs(after_peak - center) + self.assertLess( + after_dist, + before_dist, + msg="Peak did not move closer to center on second detector after alignment" + ) + + @timeout(100) + def test_cancel(self): + """ + Test cancelling alignment. + """ + f = sparc_auto_grating_offset(self.spgr, self.detector) + # Wait for the result or a timeout + try: + f.result(timeout=5) + except: + pass + self.assertTrue(f.done()) if __name__ == "__main__": unittest.main() From 9670b75eeb97d29dbd4e656c7bfa852b2a3870dc Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Thu, 2 Apr 2026 10:49:42 +0200 Subject: [PATCH 08/11] [fix] updated code for GUI and feedback --- src/odemis/acq/align/goffset.py | 292 ++++++++++++++-------- src/odemis/acq/align/test/goffset_test.py | 84 +++++-- 2 files changed, 255 insertions(+), 121 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index 047dd7fd78..6e26124c5c 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -1,29 +1,31 @@ import logging +import numpy import threading import time + from collections.abc import Iterable from concurrent.futures import CancelledError from concurrent.futures._base import CANCELLED, FINISHED, RUNNING -from typing import Any, Dict, List, Optional, Tuple, Union - -import numpy - from odemis import model from odemis.acq.align.autofocus import _mapDetectorToSelector from odemis.model import InstantaneousFuture from odemis.util import executeAsyncTask from scipy.optimize import curve_fit +from typing import Any, Dict, List, Optional, Tuple, Union +MOVE_TIME_GRATING = 20 # s +MOVE_TIME_DETECTOR = 5 # s +EST_ALIGN_TIME = 30 # s -def gaussian(x, amplitude, x0, width) -> numpy.ndarray: +def gaussian(x: numpy.ndarray, amplitude: float, x0: float, width: float) -> numpy.ndarray: """ Gaussian function (for curve fitting). - :param x: input coordinates - :param amplitude: peak intensity - :param x0 = peak's center - :param width = standard deviation - :return: Gaussian function evaluated at x (numpy array) + :param x: input coordinates + :param amplitude: peak intensity + :param x0 = peak's center + :param width = standard deviation + :return: Gaussian function evaluated at x (numpy array) """ intensity = amplitude * numpy.exp(-0.5 * ((x - x0) / width) ** 2) @@ -46,10 +48,10 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: so it has no convergence or numerical-optimization failure modes, but it may be less accurate if the peak is not well-defined or if there are multiple peaks within the window. - :param data: 1D or 2D array containing the spectrum (if 2D, it will be averaged to 1D) - :param window_radius: number of pixels on either side of the peak to include in the window for fitting (default: 15) - :return: estimated peak position in pixels (float) - :raises RuntimeError: if no significant peak is detected (SNR too low) + :param data: 1D or 2D array containing the spectrum (if 2D, it will be averaged to 1D) + :param window_radius: number of pixels on either side of the peak to include in the window for fitting (default: 15) + :return: estimated peak position in pixels (float) + :raises RuntimeError: if no significant peak is detected (SNR too low) """ if data.ndim == 2: @@ -86,7 +88,7 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: logging.info("Weighted average fallback: all window data <= 0, using peak_idx=%d as estimate", peak_idx) else: - weighted_avg = float(numpy.sum(window_idx * window_data) / numpy.sum(window_data)) + weighted_avg = float(numpy.sum(window_idx * window_data) / numpy.sum(weights)) # try Gaussian fit for better accuracy, but fallback to weighted average if it fails or is out of bounds @@ -106,7 +108,9 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: return weighted_avg -def peak_is_present(spectrum, snr_threshold=1, width_range=(0.5, 12.0)) -> bool: +def peak_is_present(spectrum: numpy.ndarray, + snr_threshold: float=1.0, + width_range: Tuple[float, float]=(0.5, 12.0)) -> bool: """ Test to decide whether spectral peak is present. @@ -115,10 +119,10 @@ def peak_is_present(spectrum, snr_threshold=1, width_range=(0.5, 12.0)) -> bool: - a local width estimate computed from a small window around the peak to reject hot pixels and extremely broad features. - :param spectrum: 1D array containing the spectrum (intensity vs pixel) - :param snr_threshold: minimum required signal-to-noise ratio for a peak to be considered present (default: 1) - :param width_range: acceptable range of estimated peak widths in pixels (default: (0.5, 12.0)) - :return: True if a peak meeting the criteria is present, False otherwise + :param spectrum: 1D array containing the spectrum (intensity vs pixel) + :param snr_threshold: minimum required signal-to-noise ratio for a peak to be considered present (default: 1) + :param width_range: acceptable range of estimated peak widths in pixels (default: (0.5, 12.0)) + :return: True if a peak meeting the criteria is present, False otherwise """ # basic stats @@ -150,7 +154,7 @@ def peak_is_present(spectrum, snr_threshold=1, width_range=(0.5, 12.0)) -> bool: return present -def acquire_peak(spgr, detector, step=2000) -> Tuple[float, float]: +def coarse_scan_goffset_for_peak(spgr, detector, step: int=2000) -> Tuple[float, float]: """ Coarse scan across the goffset axis until a real peak becomes visible. @@ -162,11 +166,11 @@ def acquire_peak(spgr, detector, step=2000) -> Tuple[float, float]: position and the function returns the actual goffset (as reported by the actuator) and the peak pixel. - :param spgr: spectrograph - :param detector: detector - :param step: step size in goffset units for the coarse scan (default: 2000) - :return: tuple (actual_goffset, peak_pixel) - :raises RuntimeError: if no peak is found across the full goffset range + :param spgr: spectrograph + :param detector: detector + :param step: step size in goffset units for the coarse scan (default: 2000) + :return: tuple (actual_goffset, peak_pixel) + :raises RuntimeError: if no peak is found across the full goffset range """ current = float(spgr.position.value["goffset"]) @@ -175,8 +179,7 @@ def acquire_peak(spgr, detector, step=2000) -> Tuple[float, float]: logging.debug( "Coarse local scan around goffset %.1f with step %.1f and max span %.1f", - current, step, max_span - ) + current, step, max_span) # Build a sequence of positions: current, +step, -step, +2*step, -2*step, ... positions = [current] @@ -211,8 +214,10 @@ def acquire_peak(spgr, detector, step=2000) -> Tuple[float, float]: raise RuntimeError("Peak not found in local goffset scan around current position") -def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta=5.0, retries=1) -> Tuple[ - float, float, float]: +def estimate_goffset_scale(spgr: model.Actuator, + detector: model.Detector, + delta: float=5.0, + retries: int=1) -> Tuple[float, float, float]: """ Estimate the scale factor between a change in the grating offset ('goffset') and the resulting shift of the spectral peak on the detector. @@ -224,16 +229,16 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta If the measured scale is unreasonably small or large, the function retries recursively and falls back to a default value of 0.5 if necessary. - :param spgr: spectrograph - :param detector: detector - :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). - The actual step may be negated to avoid exceeding hardware limits. - :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). + :param spgr: spectrograph + :param detector: detector + :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). + The actual step may be negated to avoid exceeding hardware limits. + :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). - :return: Tuple (scale, p0, p1) - scale: estimated pixels per unit of goffset - p0: peak position at the initial goffset - p1: peak position after applying the test delta + :return: Tuple (scale, p0, p1) + scale: estimated pixels per unit of goffset + p0: peak position at the initial goffset + p1: peak position after applying the test delta """ # get initial state @@ -261,8 +266,7 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta logging.info( "SCALE TRACKING | p0: %.1f | p1: %.1f | Delta: %.1f | Shift: %.1f | Result Scale: %.4f", - p0, p1, actual_delta, (p1 - p0), scale, - ) + p0, p1, actual_delta, (p1 - p0), scale) # If the estimated scale is extremely small, the measurement is likely unreliable # (e.g., due to noise or a flat spectrum). @@ -273,9 +277,7 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta if abs(scale) < 1e-3 or abs(scale) > 10.0: logging.warning( "Unreliable scale estimate (%.4f). Retries left: %d", - scale, - retries - ) + scale, retries) if retries > 0: return estimate_goffset_scale(spgr, detector, delta, retries - 1) @@ -288,26 +290,25 @@ def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, delta def sparc_auto_grating_offset(spgr: model.Actuator, detector: model.Detector, - single_detector_mode: bool = False, tolerance_px: float = 0.4, - max_it: int = 20, + max_it: int = 50, gain: float = 0.4) -> model.ProgressiveFuture: """ Start an asynchronous task that centers the spectral peak by adjusting the grating offset (goffset). - :param spgr: spectrograph - :param detector: detector - :param tolerance_px: the acceptable displacement of the peak from the center in pixels (default: 0.4) - :param max_it: maximum number of iterations to attempt (default: 20) - :param gain: proportional gain factor for adjusting the goffset (default: 0.4) - :return: A ``ProgressiveFuture`` representing the asynchronous alignment - task. The future can be used to monitor progress, retrieve the - result, or cancel the alignment. + :param spgr: spectrograph + :param detector: detector + :param tolerance_px: the acceptable displacement of the peak from the center in pixels (default: 0.4) + :param max_it: maximum number of iterations to attempt (default: 20) + :param gain: proportional gain factor for adjusting the goffset (default: 0.4) + :return: A ``ProgressiveFuture`` representing the asynchronous alignment + task. The future can be used to monitor progress, retrieve the + result, or cancel the alignment. """ est_start = time.time() + 0.05 - est_time = max_it * 0.5 # conservative estimate + est_time = max_it * 0.5 # rough estimated time f = model.ProgressiveFuture(start=est_start, end=est_start + est_time) @@ -315,18 +316,13 @@ def sparc_auto_grating_offset(spgr: model.Actuator, f._task_state = RUNNING f.task_canceller = _cancel_sparc_auto_grating_offset - executeAsyncTask( - f, - _do_sparc_auto_grating_offset, - args=(f, spgr, detector, single_detector_mode, tolerance_px, max_it, gain), - ) + executeAsyncTask(f, _do_sparc_auto_grating_offset, args=(f, spgr, detector, tolerance_px, max_it, gain)) return f def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, spgr: model.Actuator, detector: model.Detector, - single_detector_mode, tolerance_px: float, max_it: int, gain: float) -> bool: @@ -340,14 +336,14 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, - Runs a centering loop until the peak is within `tolerance_px` or the max #iterations is reached. - Respects cancellation via the provided ProgressiveFuture. - :param future: model.ProgressiveFuture - :param spgr: spectrograph - :param detector: detector - :param tolerance_px: pixel tolerance for successful alignment - :param max_it: maximum number of centering iterations - :param gain: proportional gain for converting pixel error to goffset correction - :return: True if alignment succeeded (peak within tolerance), False otherwise (bool) - :raises CancelledError: if the future was cancelled during execution + :param future: model.ProgressiveFuture + :param spgr: spectrograph + :param detector: detector + :param tolerance_px: pixel tolerance for successful alignment + :param max_it: maximum number of centering iterations + :param gain: proportional gain for converting pixel error to goffset correction + :return: True if alignment succeeded (peak within tolerance), False otherwise (bool) + :raises CancelledError: if the future was cancelled during execution """ logging.info("Running alignment | detector=%s |", detector.name) @@ -378,10 +374,9 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, # if no peak present, run acquisition (this will move the grating) if not peak_present: try: - g_acq, p_acq = acquire_peak(spgr, detector, step=2000) + g_acq, p_acq = coarse_scan_goffset_for_peak(spgr, detector, step=2000) logging.info("Peak acquired at goffset=%d pixel=%.2f", g_acq, p_acq) peak0 = p_acq - peak_present = True except RuntimeError: logging.error("Peak acquisition failed — aborting alignment") return False @@ -399,6 +394,8 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) # prefer p1 (after probe move) if available, else p0 + # p1 is preferred because it reflects the peak position after a controlled probe move, + # making it more up‑to‑date and reliable than the initial p0 measurement. start_peak = p1 if p1 is not None else p0 total_goffset_displacement = 0.0 @@ -423,15 +420,20 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, if abs(error_px) <= tolerance_px: return True + # Calculate required adjustment based on pixel error and scaling factor delta_goffset = -gain * (error_px / scale) current = spgr.position.value["goffset"] + + # Clamp move to stay within the allowed goffset range. + # Ensures that we don't command a step that would exceed axis limits. delta_goffset = max(minv - current, min(maxv - current, delta_goffset)) + + # Accumulate the actual displacement for total movement tracking total_goffset_displacement += delta_goffset logging.debug( "Iter: %d | Peak: %.2f | Error: %.2f | Move: %.6f | Total Change: %.6f", - i, peak_px, error_px, delta_goffset, total_goffset_displacement - ) + i, peak_px, error_px, delta_goffset, total_goffset_displacement) spgr.moveRelSync({"goffset": delta_goffset}) future.set_progress(end=time.time() + (max_it - i - 1) * 0.5) @@ -447,7 +449,7 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, raise -def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): +def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture) -> bool: """ Canceller of _do_sparc_auto_grating_offset task. """ @@ -455,7 +457,7 @@ def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture): future._task_state = CANCELLED -def _checkCancelled(future: "model.ProgressiveFuture"): +def _checkCancelled(future: "model.ProgressiveFuture") -> None: """ Check if the future has been cancelled, and if so raise CancelledError. """ @@ -481,7 +483,6 @@ def _total_alignment_time(n_gratings: int, # total time = time spent running alignment algorithms + time spent moving hardware return runs * EST_ALIGN_TIME + move_time - def auto_align_grating_detector_offsets(spectrograph: model.Actuator, detectors: Union[model.Detector, List[model.Detector]], selector: Optional[model.Actuator] = None, @@ -493,13 +494,13 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will be aligned by adjusting the detector offset with the grating alignment fixed. - :param spectrograph: spectrograph - :param detectors: list of detectors - :param selector: optional selector to switch between detectors - :param streams: optional list of streams to update with progress - :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) - :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector - :raises CancelledError: if the operation is cancelled + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) + :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector + :raises CancelledError: if the operation is cancelled """ if not isinstance(detectors, Iterable): @@ -512,33 +513,25 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, if streams is None: streams = [] - single_detector_mode = len(detectors) == 1 - est_start = time.time() + 0.1 n_gratings = len(spectrograph.axes["grating"].choices) n_detectors = len(detectors) a_time = _total_alignment_time(n_gratings, n_detectors) f = model.ProgressiveFuture(start=est_start, end=est_start + a_time) + f._progress = 0.0 f.task_canceller = _cancel_auto_align_grating_detector_offsets f._task_lock = threading.Lock() f._task_state = RUNNING f._subfuture = InstantaneousFuture() - executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams, single_detector_mode)) + executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams)) return f - -MOVE_TIME_GRATING = 20 # s -MOVE_TIME_DETECTOR = 5 # s -EST_ALIGN_TIME = 30 # s - - def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, spectrograph: model.Actuator, detectors: List[model.Detector], selector: Optional[model.Actuator], streams: List['Stream'], - single_detector_mode: bool = False, stabilization_time: float = 10.0) -> Optional[Dict[Any, Any]]: """ Iterate through each grating and detector combination, adjusting the selector if provided, and run the auto-alignment algorithm. @@ -547,15 +540,15 @@ def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will be aligned by adjusting the detector offset with the grating alignment fixed. - :param future: ProgressiveFuture to update with progress and results - :param spectrograph: spectrograph - :param detectors: list of detectors - :param selector: optional selector to switch between detectors - :param streams: optional list of streams to update with progress - :param stabilization_time: time to wait after moving hardware before starting alignment (default: 15s) + :param future: ProgressiveFuture to update with progress and results + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :param stabilization_time: time to wait after moving hardware before starting alignment (default: 15s) - :return: dict mapping (grating, detector) to alignment success boolean - :raises CancelledError: if the operation is cancelled + :return: dict mapping (grating, detector) to alignment success boolean + :raises CancelledError: if the operation is cancelled """ results: dict[tuple, bool] = {} @@ -565,23 +558,72 @@ def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, gratings = sorted(list(spectrograph.axes["grating"].choices.keys())) logging.info(f"Available gratings: {list(spectrograph.axes['grating'].choices.keys())}") + total_time = _total_alignment_time(len(gratings), len(detectors)) + start_time = time.time() + + # helper functions for progress tracking in GUI + def update_progress(): + elapsed = time.time() - start_time + future._progress = min(1.0, elapsed / total_time) + + def set_step(duration): + future._step_start_time = time.time() + future._step_duration = duration + + # main alignment loop first_detector = detectors[0] if selector: original_selector = selector.position.value selector_axes, detector_to_selector = _mapDetectorToSelector(selector, detectors) + logging.debug("selector_axes=%s detector_to_selector keys=%s", + selector_axes, list(detector_to_selector.keys())) def is_current_detector(d): if selector is None: return True return detector_to_selector[d] == selector.position.value[selector_axes] + # Store original states of detectors to restore later + # Force GUI to use full resolution + original_detector_states = {} + for d in detectors: + state = {} + if hasattr(d, 'translation'): state['translation'] = d.translation.value + if hasattr(d, 'binning'): state['binning'] = d.binning.value + if hasattr(d, 'resolution'): state['resolution'] = d.resolution.value + original_detector_states[d.name] = state + + logging.info("Forcing detector %s to full horizontal sensor", d.name) + try: + # 1. Reset ROI translation to top-left corner + if hasattr(d, 'translation'): + d.translation.value = (0, 0) + + # 2. Maximize vertical binning (squashes light into a bright 1D line) + if hasattr(d, 'binning'): + max_bin_y = d.binning.range[1][1] + d.binning.value = (1, max_bin_y) + + # 3. Read the full horizontal width of the sensor + if hasattr(d, 'resolution'): + max_res_x = d.resolution.range[1][0] + d.resolution.value = (max_res_x, 1) + except Exception as e: + logging.warning("Failed to enforce full-sensor state for %s: %s", d.name, e) + + # start alignment for the first grating try: g0 = gratings[0] logging.info("Starting alignment for initial grating: %s", g0) + set_step(MOVE_TIME_GRATING) spectrograph.moveAbsSync({"grating": g0, "wavelength": 0}) + update_progress() + + set_step(stabilization_time) time.sleep(stabilization_time) + update_progress() detectors_sorted = sorted(detectors, key=is_current_detector, reverse=True) @@ -591,31 +633,52 @@ def is_current_detector(d): logging.info("Starting alignment | Detector: %s | Grating: %s", d.name, g0) if selector: + set_step(MOVE_TIME_DETECTOR) selector.moveAbsSync({selector_axes: detector_to_selector[d]}) - future._subfuture = sparc_auto_grating_offset(spectrograph, d, single_detector_mode=single_detector_mode) - success = future._subfuture.result() - results[(g0, d.name)] = success + update_progress() - logging.info("Finished alignment | Detector: %s | Grating: %s", d.name, g0) + set_step(stabilization_time) + time.sleep(stabilization_time) + update_progress() + + future._subfuture = sparc_auto_grating_offset(spectrograph, d) + success = future._subfuture.result() + results[(g0, d.name)] = success + + logging.info("Finished alignment | Detector: %s | Grating: %s", d.name, g0) if selector: + set_step(MOVE_TIME_DETECTOR) selector.moveAbsSync({selector_axes: detector_to_selector[first_detector]}) + update_progress() + + set_step(stabilization_time) + time.sleep(stabilization_time) + update_progress() # align remaining gratings using the first detector for g in gratings[1:]: _checkCancelled(future) logging.info("Switching to grating: %s", g) + set_step(MOVE_TIME_GRATING) spectrograph.moveAbsSync({"grating": g, "wavelength": 0}) + update_progress() + + set_step(stabilization_time) time.sleep(stabilization_time) + update_progress() logging.info("Starting alignment | Detector: %s | Grating: %s", first_detector.name, g) - future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector, single_detector_mode=single_detector_mode) + set_step(EST_ALIGN_TIME) + future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector) success = future._subfuture.result() results[(g, first_detector.name)] = success + update_progress() logging.info("Finished alignment | Detector: %s | Grating: %s", first_detector.name, g) + future._progress = 1.0 return results except CancelledError: @@ -627,10 +690,25 @@ def is_current_detector(d): if selector: selector.moveAbsSync(original_selector) + # Restore original GUI settings + for d in detectors: + if d.name in original_detector_states: + state = original_detector_states[d.name] + logging.info("Restoring original state for detector %s", d.name) + try: + # Restore in reverse order: resolution, binning, then offset + if 'resolution' in state and hasattr(d, 'resolution'): + d.resolution.value = state['resolution'] + if 'binning' in state and hasattr(d, 'binning'): + d.binning.value = state['binning'] + if 'offset' in state and hasattr(d, 'offset'): + d.offset.value = state['offset'] + except Exception as e: + logging.warning("Could not fully restore state for %s: %s", d.name, e) + with future._task_lock: future._task_state = FINISHED - def _cancel_auto_align_grating_detector_offsets(future: model.ProgressiveFuture) -> bool: """ Canceller for _do_auto_align_grating_detector_offsets task. diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index a01a72d9b8..e07ddeeb50 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -140,8 +140,7 @@ def test_scale_estimation_misaligned(self): start_goffset, end_goffset, places=3, - msg="goffset did not change during alignment when peak was misaligned" - ) + msg="goffset did not change during alignment when peak was misaligned") @timeout(800) def test_auto_grating_offset(self): @@ -223,7 +222,7 @@ def test_multi_detector_iteration(self): def test_single_detector_updates_grating(self): """ Single-detector mode: aligning a detector on the secondary port must update - the *grating* goffset (not a detector-specific offset). + the grating offset and not the detector offset. """ spccd = model.getComponent(role="sp-ccd") spccd.exposureTime.value = spccd.exposureTime.range[0] @@ -252,7 +251,7 @@ def test_single_detector_updates_grating(self): def test_multi_detector_does_not_change_grating(self): """ Multi-detector mode: first detector sets the grating; aligning the second - detector should not change the grating goffset (it should adjust detector offset). + detector should not change the grating offset (it should adjust detector offset). """ spccd = model.getComponent(role="sp-ccd") spccd.exposureTime.value = spccd.exposureTime.range[0] @@ -264,7 +263,7 @@ def test_multi_detector_does_not_change_grating(self): logging.debug("Selector move to primary failed or not present; continuing") # align first detector in single-detector mode to set grating - f_first = sparc_auto_grating_offset(self.spgr, self.detector, single_detector_mode=True, max_it=50) + f_first = sparc_auto_grating_offset(self.spgr, self.detector, max_it=50) self.assertTrue(f_first.result(timeout=300), "First-detector alignment failed") grating_after_first = self.spgr.position.value["goffset"] @@ -275,7 +274,7 @@ def test_multi_detector_does_not_change_grating(self): logging.debug("Selector move to secondary failed or not present; continuing") # align second detector in multi-detector mode (should not change grating) - f_second = sparc_auto_grating_offset(self.spgr, spccd, single_detector_mode=False, max_it=50) + f_second = sparc_auto_grating_offset(self.spgr, spccd, max_it=50) self.assertTrue(f_second.result(timeout=300), "Second-detector alignment failed") grating_after_second = self.spgr.position.value["goffset"] @@ -283,8 +282,7 @@ def test_multi_detector_does_not_change_grating(self): grating_after_first, grating_after_second, places=3, - msg="Grating goffset changed when aligning second detector in multi-detector mode" - ) + msg="Grating goffset changed when aligning second detector in multi-detector mode") @timeout(900) def test_multi_detector_detector_offset_changes(self): @@ -302,7 +300,7 @@ def test_multi_detector_detector_offset_changes(self): except Exception: logging.debug("Selector move to primary failed or not present; continuing") - f_first = sparc_auto_grating_offset(self.spgr, self.detector, single_detector_mode=True, max_it=50) + f_first = sparc_auto_grating_offset(self.spgr, self.detector, max_it=50) self.assertTrue(f_first.result(timeout=300), "First-detector alignment failed") grating_after_first = self.spgr.position.value["goffset"] @@ -317,7 +315,7 @@ def test_multi_detector_detector_offset_changes(self): before_peak = float(find_peak_position(data_before)) # Run alignment for second detector in multi-detector mode - f_second = sparc_auto_grating_offset(self.spgr, spccd, single_detector_mode=False, max_it=50) + f_second = sparc_auto_grating_offset(self.spgr, spccd, max_it=50) self.assertTrue(f_second.result(timeout=600), "Second-detector alignment failed") # Measure peak after alignment @@ -330,8 +328,7 @@ def test_multi_detector_detector_offset_changes(self): grating_after_first, grating_after_second, places=3, - msg="Grating goffset changed when aligning second detector in multi-detector mode" - ) + msg="Grating goffset changed when aligning second detector in multi-detector mode") # Assert peak moved closer to center on the secondary detector center = spccd.resolution.value[0] / 2 @@ -340,8 +337,67 @@ def test_multi_detector_detector_offset_changes(self): self.assertLess( after_dist, before_dist, - msg="Peak did not move closer to center on second detector after alignment" - ) + msg="Peak did not move closer to center on second detector after alignment") + + def test_single_detector_iteration(self): + """ + Verifies that the auto-alignment algorithm generates the minimum required set + of offsets when operating with a single detector and multiple gratings. + """ + + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + res = f.result(timeout=900) + + n_gratings = len(self.spgr.axes["grating"].choices) + n_detectors = 1 + expected = n_detectors + (n_gratings - 1) + + self.assertEqual(len(res), expected) + + first_grating = list(self.spgr.axes["grating"].choices.keys())[0] + dets_first = [d for (g, d) in res.keys() if g == first_grating] + + self.assertEqual(len(dets_first), n_detectors) + + def test_multi_detector_iteration(self): + spccd = model.getComponent(role="sp-ccd") + spccd.exposureTime.value = spccd.exposureTime.range[0] + + detectors = [self.ccd, spccd] + + # run alignment + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) + res = f.result(timeout=900) + + # calculate expected results + n_gratings = len(self.spgr.axes["grating"].choices) + n_detectors = len(detectors) + expected_count = n_detectors + (n_gratings - 1) + + self.assertEqual(len(res), expected_count, f"Expected {expected_count} results, got {len(res)}") + + # verify that every detector was used for the first grating + gratings_list = list(self.spgr.axes["grating"].choices.keys()) + first_grating = gratings_list[0] + + dets_for_first_grating = [d for (g, d) in res.keys() if g == first_grating] + self.assertEqual(len(dets_for_first_grating), n_detectors) + self.assertIn(self.ccd.name, dets_for_first_grating) + self.assertIn(spccd.name, dets_for_first_grating) + + # verify that only first detector is used for remaining gratings + for g in gratings_list[1:]: + # Ensure only the first detector (index 0) is present for these gratings + dets_for_this_grating = [d for (grating, d) in res.keys() if grating == g] + self.assertEqual(len(dets_for_this_grating), 1) + self.assertEqual(dets_for_this_grating[0], detectors[0].name) + + # move to spectral camera + self.selector.moveAbsSync({"rx": 1.5707963267948966}) + data = spccd.data.get(asap=False) + + # check data is not flat + self.assertNotEqual(data.max(), data.min) @timeout(100) def test_cancel(self): From 07f936e80db73b20336fea7441c91735f3897af0 Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Tue, 14 Apr 2026 13:38:15 +0200 Subject: [PATCH 09/11] [fix] handled ValueError + refined peak detection --- src/odemis/acq/align/goffset.py | 224 +++++++++-------- src/odemis/acq/align/test/goffset_test.py | 249 ++++++++++++------- src/odemis/gui/cont/tabs/sparc2_align_tab.py | 98 ++++++++ 3 files changed, 371 insertions(+), 200 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index 6e26124c5c..6dce28e169 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -12,6 +12,8 @@ from odemis.util import executeAsyncTask from scipy.optimize import curve_fit from typing import Any, Dict, List, Optional, Tuple, Union +from scipy.ndimage import gaussian_filter1d + MOVE_TIME_GRATING = 20 # s MOVE_TIME_DETECTOR = 5 # s @@ -59,18 +61,25 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: else: spectrum = data - # compute simple SNR by comparing the peak height above the median background to the background magnitude and returns + # compute SNR by comparing the peak height above the median background to the background magnitude and returns # False if that SNR is below the threshold. This helps reject cases with no real peak, # where the Gaussian fit would fail or produce nonsense results. - peak_value = float(spectrum.max()) - background = float(numpy.median(spectrum)) - snr = (peak_value - background) / (abs(background) + 1e-6) + # peak_value = float(spectrum.max()) + # background = float(numpy.median(spectrum)) + # snr = (peak_value - background) / (abs(background) + 1e-6) + + spectrum = spectrum - numpy.median(spectrum) + noise_std = numpy.std(spectrum) + peak_height = spectrum.max() + snr = peak_height / (noise_std + 1e-6) - if snr < 1: # tune threshold + if snr < 10.0: # tune threshold raise RuntimeError("No peak detected (SNR too low)") - peak_idx = numpy.argmax(spectrum) # find the absolute highest point + spectrum_smooth = gaussian_filter1d(spectrum, sigma=2) + peak_idx = int(numpy.argmax(spectrum_smooth)) + # peak_idx = numpy.argmax(spectrum) # find the absolute highest point # create a window around the peak start = max(0, peak_idx - window_radius) @@ -109,7 +118,7 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: return weighted_avg def peak_is_present(spectrum: numpy.ndarray, - snr_threshold: float=1.0, + snr_threshold: float=5.0, width_range: Tuple[float, float]=(0.5, 12.0)) -> bool: """ Test to decide whether spectral peak is present. @@ -126,15 +135,22 @@ def peak_is_present(spectrum: numpy.ndarray, """ # basic stats - peak_value = float(spectrum.max()) - background = float(numpy.median(spectrum)) - snr = (peak_value - background) / (abs(background) + 1e-6) + # peak_value = float(spectrum.max()) + # background = float(numpy.median(spectrum)) + # snr = (peak_value - background) / (abs(background) + 1e-6) + + spectrum = spectrum - numpy.median(spectrum) + noise_std = numpy.std(spectrum) + peak_value = spectrum.max() + snr = peak_value / (noise_std + 1e-6) if snr < snr_threshold: return False # estimate width around the peak - peak_idx = int(numpy.argmax(spectrum)) + spectrum_smooth = gaussian_filter1d(spectrum, sigma=2) + peak_idx = int(numpy.argmax(spectrum_smooth)) + # peak_idx = int(numpy.argmax(spectrum)) if peak_idx < 1 or peak_idx > len(spectrum) - 2: logging.debug("Peak too close to edge idx=%d len=%d", peak_idx, len(spectrum)) return False @@ -154,7 +170,7 @@ def peak_is_present(spectrum: numpy.ndarray, return present -def coarse_scan_goffset_for_peak(spgr, detector, step: int=2000) -> Tuple[float, float]: +def coarse_scan_goffset_for_peak(spgr, detector, future: model.ProgressiveFuture, step: int=2000) -> Tuple[float, float]: """ Coarse scan across the goffset axis until a real peak becomes visible. @@ -175,7 +191,7 @@ def coarse_scan_goffset_for_peak(spgr, detector, step: int=2000) -> Tuple[float, current = float(spgr.position.value["goffset"]) step = abs(step) if step != 0 else 2000.0 - max_span = 20000.0 # limit how far we wander from the current valid position + max_span = 200000.0 # limit how far we wander from the current valid position logging.debug( "Coarse local scan around goffset %.1f with step %.1f and max span %.1f", @@ -192,6 +208,8 @@ def coarse_scan_goffset_for_peak(spgr, detector, step: int=2000) -> Tuple[float, tried = set() for g in positions: + _checkCancelled(future) + # avoid duplicate moves due to symmetry or float rounding key = round(g, 3) if key in tried: @@ -201,8 +219,8 @@ def coarse_scan_goffset_for_peak(spgr, detector, step: int=2000) -> Tuple[float, logging.debug("Coarse scan move attempt to goffset %.3f", g) try: spgr.moveAbsSync({"goffset": g}) - except Exception as e: - logging.warning("Skipping invalid goffset %.3f (%s)", g, e) + except ValueError: + logging.warning("Skipping invalid goffset %.3f (%s)", g, ValueError) continue data = detector.data.get(asap=False) @@ -287,11 +305,37 @@ def estimate_goffset_scale(spgr: model.Actuator, return scale, p0, p1 +def log_detector_state(caller: str, stage: str, detector: model.Detector, data: numpy.ndarray): + shape = data.shape + + # Keep consistent with detection pipeline + spectrum = data.mean(axis=0) if data.ndim == 2 else data + + max_val = float(spectrum.max()) + min_val = float(spectrum.min()) + mean_val = float(spectrum.mean()) + + background = float(numpy.median(spectrum)) + signal = spectrum - background + peak_height = float(signal.max()) + + # robust noise estimate + noise_region = signal[signal < 0] + if len(noise_region) > 10: + noise_std = float(numpy.std(noise_region)) + else: + noise_std = float(numpy.std(signal)) + + snr = peak_height / (noise_std + 1e-6) + + logging.warning( + "%s: goffset: [%s] Detector=%s | shape=%s | min=%.2f max=%.2f mean=%.2f | bg=%.2f | noise=%.2f | snr=%.2f", + caller, stage, detector.name, shape, min_val, max_val, mean_val, background, noise_std, snr) def sparc_auto_grating_offset(spgr: model.Actuator, detector: model.Detector, tolerance_px: float = 0.4, - max_it: int = 50, + max_it: int = 60, gain: float = 0.4) -> model.ProgressiveFuture: """ Start an asynchronous task that centers the spectral peak by adjusting the @@ -326,7 +370,6 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, tolerance_px: float, max_it: int, gain: float) -> bool: - """ Core alignment routine that iteratively adjusts the grating offset to center the peak. @@ -350,18 +393,29 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, try: center_target = detector.resolution.value[0] / 2 + data0 = detector.data.get(asap=False) + spectrum0 = data0.mean(axis=0) if data0.ndim == 2 else data0 - # initial read: try to get a valid peak without moving the grating - try: - data0 = detector.data.get(asap=False) - peak0 = find_peak_position(data0) # raises RuntimeError if no peak present + if peak_is_present(spectrum0): + peak0 = find_peak_position(data0) logging.debug("Initial read: peak0=%.2f px", peak0) peak_present = True - except RuntimeError: + else: logging.debug("Initial read: no significant peak detected, scanning the axes") peak_present = False peak0 = None + # # initial read: try to get a valid peak without moving the grating + # try: + # data0 = detector.data.get(asap=False) + # peak0 = find_peak_position(data0) # raises RuntimeError if no peak present + # logging.debug("Initial read: peak0=%.2f px", peak0) + # peak_present = True + # except RuntimeError: + # logging.debug("Initial read: no significant peak detected, scanning the axes") + # peak_present = False + # peak0 = None + # if peak present and already centered, do nothing if peak_present: initial_error_px = peak0 - center_target @@ -372,9 +426,9 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, return True # if no peak present, run acquisition (this will move the grating) - if not peak_present: + else: try: - g_acq, p_acq = coarse_scan_goffset_for_peak(spgr, detector, step=2000) + g_acq, p_acq = coarse_scan_goffset_for_peak(spgr, detector, future, step=2000) logging.info("Peak acquired at goffset=%d pixel=%.2f", g_acq, p_acq) peak0 = p_acq except RuntimeError: @@ -390,8 +444,20 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, return True # peak is present and not centered -> estimate scale and center - scale, p0, p1 = estimate_goffset_scale(spgr, detector) - logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) + # scale, p0, p1 = estimate_goffset_scale(spgr, detector) + # logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) + + try: + scale, p0, p1 = estimate_goffset_scale(spgr, detector) + logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) + except RuntimeError: + logging.warning("Peak lost during scale estimation. Forcing re-acquisition.") + try: + g_acq, p_acq = coarse_scan_goffset_for_peak(spgr, detector, future, step=2000) + scale, p0, p1 = estimate_goffset_scale(spgr, detector) # try one more time + except RuntimeError: + logging.error("Peak re-acquisition failed — aborting alignment") + return False # prefer p1 (after probe move) if available, else p0 # p1 is preferred because it reflects the peak position after a controlled probe move, @@ -412,7 +478,7 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, try: peak_px = find_peak_position(data) except RuntimeError: - logging.error("No peak detected during centering iteration %d — aborting", i) + logging.error("Peak re-acquisition failed during iteration %d — aborting", i) return False error_px = peak_px - center_target @@ -435,7 +501,12 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, "Iter: %d | Peak: %.2f | Error: %.2f | Move: %.6f | Total Change: %.6f", i, peak_px, error_px, delta_goffset, total_goffset_displacement) - spgr.moveRelSync({"goffset": delta_goffset}) + try: + spgr.moveRelSync({"goffset": delta_goffset}) + except ValueError: + logging.warning("Hardware offset limit reached: %s. Keeping max allowed move.", ValueError) + break + future.set_progress(end=time.time() + (max_it - i - 1) * 0.5) logging.warning("SparcAutoGratingOffset did not converge within max iterations") @@ -558,19 +629,6 @@ def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, gratings = sorted(list(spectrograph.axes["grating"].choices.keys())) logging.info(f"Available gratings: {list(spectrograph.axes['grating'].choices.keys())}") - total_time = _total_alignment_time(len(gratings), len(detectors)) - start_time = time.time() - - # helper functions for progress tracking in GUI - def update_progress(): - elapsed = time.time() - start_time - future._progress = min(1.0, elapsed / total_time) - - def set_step(duration): - future._step_start_time = time.time() - future._step_duration = duration - - # main alignment loop first_detector = detectors[0] if selector: @@ -584,46 +642,18 @@ def is_current_detector(d): return True return detector_to_selector[d] == selector.position.value[selector_axes] - # Store original states of detectors to restore later - # Force GUI to use full resolution - original_detector_states = {} - for d in detectors: - state = {} - if hasattr(d, 'translation'): state['translation'] = d.translation.value - if hasattr(d, 'binning'): state['binning'] = d.binning.value - if hasattr(d, 'resolution'): state['resolution'] = d.resolution.value - original_detector_states[d.name] = state - - logging.info("Forcing detector %s to full horizontal sensor", d.name) - try: - # 1. Reset ROI translation to top-left corner - if hasattr(d, 'translation'): - d.translation.value = (0, 0) - - # 2. Maximize vertical binning (squashes light into a bright 1D line) - if hasattr(d, 'binning'): - max_bin_y = d.binning.range[1][1] - d.binning.value = (1, max_bin_y) - - # 3. Read the full horizontal width of the sensor - if hasattr(d, 'resolution'): - max_res_x = d.resolution.range[1][0] - d.resolution.value = (max_res_x, 1) - except Exception as e: - logging.warning("Failed to enforce full-sensor state for %s: %s", d.name, e) - - # start alignment for the first grating + # Calculate total steps for a simple progress bar + total_steps = len(detectors) + (len(gratings) - 1) + current_step = 0 + future._progress = 0.0 + + # Start alignment for the first grating try: g0 = gratings[0] logging.info("Starting alignment for initial grating: %s", g0) - set_step(MOVE_TIME_GRATING) spectrograph.moveAbsSync({"grating": g0, "wavelength": 0}) - update_progress() - - set_step(stabilization_time) time.sleep(stabilization_time) - update_progress() detectors_sorted = sorted(detectors, key=is_current_detector, reverse=True) @@ -633,48 +663,42 @@ def is_current_detector(d): logging.info("Starting alignment | Detector: %s | Grating: %s", d.name, g0) if selector: - set_step(MOVE_TIME_DETECTOR) selector.moveAbsSync({selector_axes: detector_to_selector[d]}) - update_progress() - - set_step(stabilization_time) time.sleep(stabilization_time) - update_progress() + + gui_data = d.data.get(asap=False) + log_detector_state("GUI", "INITIAL", d, gui_data) future._subfuture = sparc_auto_grating_offset(spectrograph, d) success = future._subfuture.result() results[(g0, d.name)] = success + # Update progress + current_step += 1 + future._progress = current_step / total_steps + logging.info("Finished alignment | Detector: %s | Grating: %s", d.name, g0) if selector: - set_step(MOVE_TIME_DETECTOR) selector.moveAbsSync({selector_axes: detector_to_selector[first_detector]}) - update_progress() - - set_step(stabilization_time) time.sleep(stabilization_time) - update_progress() # align remaining gratings using the first detector for g in gratings[1:]: _checkCancelled(future) logging.info("Switching to grating: %s", g) - set_step(MOVE_TIME_GRATING) spectrograph.moveAbsSync({"grating": g, "wavelength": 0}) - update_progress() - - set_step(stabilization_time) time.sleep(stabilization_time) - update_progress() logging.info("Starting alignment | Detector: %s | Grating: %s", first_detector.name, g) - set_step(EST_ALIGN_TIME) future._subfuture = sparc_auto_grating_offset(spectrograph, first_detector) success = future._subfuture.result() results[(g, first_detector.name)] = success - update_progress() + + # Update progress + current_step += 1 + future._progress = current_step / total_steps logging.info("Finished alignment | Detector: %s | Grating: %s", first_detector.name, g) @@ -690,22 +714,6 @@ def is_current_detector(d): if selector: selector.moveAbsSync(original_selector) - # Restore original GUI settings - for d in detectors: - if d.name in original_detector_states: - state = original_detector_states[d.name] - logging.info("Restoring original state for detector %s", d.name) - try: - # Restore in reverse order: resolution, binning, then offset - if 'resolution' in state and hasattr(d, 'resolution'): - d.resolution.value = state['resolution'] - if 'binning' in state and hasattr(d, 'binning'): - d.binning.value = state['binning'] - if 'offset' in state and hasattr(d, 'offset'): - d.offset.value = state['offset'] - except Exception as e: - logging.warning("Could not fully restore state for %s: %s", d.name, e) - with future._task_lock: future._task_state = FINISHED diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index e07ddeeb50..044e1769ba 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -11,18 +11,28 @@ from odemis import model from odemis.util import timeout -from odemis.acq.align.goffset import( +from odemis.acq.align.goffset_ext import( find_peak_position, + peak_is_present, estimate_goffset_scale, sparc_auto_grating_offset, - auto_align_grating_detector_offsets + auto_align_grating_detector_offsets, + log_detector_state ) + + +from odemis.dataio import hdf5 + + logging.getLogger().setLevel(logging.DEBUG) CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" +HOME_PATH = os.path.expanduser("~") + "/" +H5_FILE_2d_NO_PEAK = HOME_PATH + "development/odemis/grating 2 1024x256/" + class TestSparcAutoGratingOffset(unittest.TestCase): """ Test automatic grating offset alignment. @@ -106,12 +116,8 @@ def test_scale_not_misaligned(self): # If peak is already centered, the algorithm exits immediately # so goffset should not change. - self.assertAlmostEqual( - start_goffset, - end_goffset, - places=6, - msg="goffset changed even though peak was already centered (scale estimation likely ran)" - ) + self.assertAlmostEqual(start_goffset, end_goffset, places=6, + msg="goffset changed even though peak was already centered (scale estimation likely ran)") def test_scale_estimation_misaligned(self): """ @@ -136,94 +142,37 @@ def test_scale_estimation_misaligned(self): self.assertTrue(result) # If misaligned, centering should move the grating - self.assertNotAlmostEqual( - start_goffset, - end_goffset, - places=3, + self.assertNotAlmostEqual(start_goffset, end_goffset, places=3, msg="goffset did not change during alignment when peak was misaligned") @timeout(800) def test_auto_grating_offset(self): """ - Test automatic centering of spectral peak. + Test automatic centering of spectral peak (one grating/detector only). """ + delta = 0 # intentionally misalign current = self.spgr.position.value["goffset"] goffset_max = self.spgr.axes["goffset"].range[1] direction = 1 if (current + delta < goffset_max) else -1 + test_data = self.detector.data.get(asap=False) + log_detector_state("TEST", "INITIAL", self.detector, test_data) + self.spgr.moveRelSync({"goffset": delta * direction}) logging.info("Test: after misalign move, spgr.position.gooffset = %s", self.spgr.position.value["goffset"]) - f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=50) + f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=100) result = f.result(timeout=800) self.assertTrue(result) - @timeout(1000) - def test_single_detector_iteration(self): - f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.detector], selector=self.selector) - res = f.result(timeout=900) - - n_gratings = len(self.spgr.axes["grating"].choices) - n_detectors = 1 - expected = n_detectors + (n_gratings - 1) - - self.assertEqual(len(res), expected) - - first_grating = list(self.spgr.axes["grating"].choices.keys())[0] - dets_first = [d for (g, d) in res.keys() if g == first_grating] - - self.assertEqual(len(dets_first), n_detectors) - - def test_multi_detector_iteration(self): - spccd = model.getComponent(role="sp-ccd") - spccd.exposureTime.value = spccd.exposureTime.range[0] - - detectors = [self.ccd, spccd] - - # run alignment - f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) - res = f.result(timeout=900) - - # calculate expected results - n_gratings = len(self.spgr.axes["grating"].choices) - n_detectors = len(detectors) - expected_count = n_detectors + (n_gratings - 1) - - self.assertEqual(len(res), expected_count, f"Expected {expected_count} results, got {len(res)}") - - # verify that every detector was used for the first grating - gratings_list = list(self.spgr.axes["grating"].choices.keys()) - first_grating = gratings_list[0] - - dets_for_first_grating = [d for (g, d) in res.keys() if g == first_grating] - self.assertEqual(len(dets_for_first_grating), n_detectors) - self.assertIn(self.ccd.name, dets_for_first_grating) - self.assertIn(spccd.name, dets_for_first_grating) - - # verify that only first detector is used for remaining gratings - for g in gratings_list[1:]: - # Ensure only the first detector (index 0) is present for these gratings - dets_for_this_grating = [d for (grating, d) in res.keys() if grating == g] - self.assertEqual(len(dets_for_this_grating), 1) - self.assertEqual(dets_for_this_grating[0], detectors[0].name) - - # move to spectral camera - self.selector.moveAbsSync({"rx": 1.5707963267948966}) - data = spccd.data.get(asap=False) - - # check data is not flat - if data.max() == data.min(): - print("WARNING: sp-ccd is returning a flat image!") - else: - print(f"sp-ccd signal range: {data.min()} to {data.max()}") - @timeout(600) def test_single_detector_updates_grating(self): """ Single-detector mode: aligning a detector on the secondary port must update the grating offset and not the detector offset. """ + spccd = model.getComponent(role="sp-ccd") spccd.exposureTime.value = spccd.exposureTime.range[0] @@ -240,12 +189,8 @@ def test_single_detector_updates_grating(self): self.assertTrue(result, "Single-detector alignment failed") after_goffset = self.spgr.position.value["goffset"] - self.assertNotAlmostEqual( - start_goffset, - after_goffset, - places=3, - msg="Grating goffset did not change for single-detector alignment on secondary port" - ) + self.assertNotAlmostEqual(start_goffset, after_goffset, places=3, + msg="Grating goffset did not change for single-detector alignment on secondary port") @timeout(900) def test_multi_detector_does_not_change_grating(self): @@ -253,6 +198,7 @@ def test_multi_detector_does_not_change_grating(self): Multi-detector mode: first detector sets the grating; aligning the second detector should not change the grating offset (it should adjust detector offset). """ + spccd = model.getComponent(role="sp-ccd") spccd.exposureTime.value = spccd.exposureTime.range[0] @@ -262,7 +208,7 @@ def test_multi_detector_does_not_change_grating(self): except Exception: logging.debug("Selector move to primary failed or not present; continuing") - # align first detector in single-detector mode to set grating + # align first detector to set grating offset f_first = sparc_auto_grating_offset(self.spgr, self.detector, max_it=50) self.assertTrue(f_first.result(timeout=300), "First-detector alignment failed") grating_after_first = self.spgr.position.value["goffset"] @@ -278,10 +224,7 @@ def test_multi_detector_does_not_change_grating(self): self.assertTrue(f_second.result(timeout=300), "Second-detector alignment failed") grating_after_second = self.spgr.position.value["goffset"] - self.assertAlmostEqual( - grating_after_first, - grating_after_second, - places=3, + self.assertAlmostEqual(grating_after_first, grating_after_second, places=3, msg="Grating goffset changed when aligning second detector in multi-detector mode") @timeout(900) @@ -324,19 +267,14 @@ def test_multi_detector_detector_offset_changes(self): # Assert grating did not change grating_after_second = self.spgr.position.value["goffset"] - self.assertAlmostEqual( - grating_after_first, - grating_after_second, - places=3, + self.assertAlmostEqual(grating_after_first, grating_after_second, places=3, msg="Grating goffset changed when aligning second detector in multi-detector mode") # Assert peak moved closer to center on the secondary detector center = spccd.resolution.value[0] / 2 before_dist = abs(before_peak - center) after_dist = abs(after_peak - center) - self.assertLess( - after_dist, - before_dist, + self.assertLess(after_dist, before_dist, msg="Peak did not move closer to center on second detector after alignment") def test_single_detector_iteration(self): @@ -345,7 +283,9 @@ def test_single_detector_iteration(self): of offsets when operating with a single detector and multiple gratings. """ - f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + #f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=self.detector, selector=self.selector) + res = f.result(timeout=900) n_gratings = len(self.spgr.axes["grating"].choices) @@ -399,6 +339,131 @@ def test_multi_detector_iteration(self): # check data is not flat self.assertNotEqual(data.max(), data.min) + def test_driver_raises_valueerror_on_hardware_limit(self): + """ + Verifies the driver successfully raises a ValueError when computing + an offset outside of the hardware limits. + """ + # A value guaranteed to exceed DET_OFFSET_MAX + GRAT_OFFSET_MAX + out_of_bounds_target = 9999999.0 + + # Test method directly to ensure the new bounds-checking logic works + with self.assertRaises(ValueError, msg="Driver should raise ValueError for out-of-bounds offset"): + self.spgr._doSetGoffsetAbs(out_of_bounds_target) + + @timeout(150) + def test_hardware_limit_valueerror_handled(self): + """ + Verifies that when a goffset move exceeds hardware limits and throws a ValueError, + the alignment algorithm catches it gracefully instead of crashing. + """ + # Slightly misalign the peak so the algorithm wants to move + self.spgr.moveRelSync({"goffset": 50}) + + # Temporarily bypass the software axis range clamp to ensure the raw command + # reaches the hardware driver to trigger the ValueError + original_range = self.spgr.axes["goffset"].range + + try: + self.spgr.axes["goffset"].range = (-9999999.0, 9999999.0) + + # Execute auto-alignment procedure + f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=10, gain=99999.0) + result = f.result(timeout=100) + + # Driver should raise ValueError for out-of-bounds move, which should be caught by the alignment code, + self.assertFalse(result, "Alignment should gracefully fail (return False) when blocked by hardware limits.") + + finally: + # Restore original software limits and position + self.spgr.axes["goffset"].range = original_range + self.spgr.moveAbsSync({"goffset": self._original_position["goffset"]}) + + def test_find_peak_position_realistic_2d(self): + """ + Test peak detection on realistic 2D detector data (like GUI), + including noise and baseline offset. + """ + np.random.seed(0) + + x = np.arange(1024) + true_center = 512.3 + + # base gaussian peak + peak = 50 * np.exp(-0.5 * ((x - true_center) / 6.0) ** 2) + + # simulate camera baseline + baseline = 300 + + # build 2D image with row variations + rows = 256 + image = [] + + for i in range(rows): + row_noise = np.random.normal(0, 3, size=x.shape) # noise + row_gain = 1 + np.random.normal(0, 0.02) # slight variation + row = baseline + row_gain * peak + row_noise + image.append(row) + + image = np.array(image) + + # run your actual function + peak_pos = find_peak_position(image) + + self.assertAlmostEqual(peak_pos, true_center, places=1) + + def test_no_peak_present_image(self): + im_no_peak = hdf5.read_data(H5_FILE_2d_NO_PEAK + "no-peak-2d-1.h5") + data = (im_no_peak[0].data) + + spectrum = data.mean(axis=0) if data.ndim == 2 else data + + present = peak_is_present(spectrum, snr_threshold=10.0, width_range=(0.5, 12.0)) + + self.assertFalse(present, f"Peak was not detected. Check SNR and width limits.") + + def test_peak_present_image(self): + im_peak = hdf5.read_data(H5_FILE_2d_NO_PEAK + "with-peak-2d-1.h5") + raw_data = np.asarray(im_peak[0].data) + + spectrum = np.squeeze(raw_data) + if spectrum.ndim > 1: + # Average over the row dimension (256) to get the 1D spectrum (1024) + spectrum = spectrum.mean(axis=0) + + clean_spec = spectrum - np.median(spectrum) + noise_std = np.std(clean_spec) + peak_val = clean_spec.max() + snr = peak_val / (noise_std + 1e-6) + peak_idx = np.argmax(clean_spec) + + # Simple width estimation (standard deviation around peak) + window = clean_spec[max(0, peak_idx - 2): min(len(clean_spec), peak_idx + 3)] + x = np.arange(len(window)) + w = window - window.min() + width = 0.0 + if w.sum() > 0: + mean = np.sum(x * w) / np.sum(w) + var = np.sum(w * (x - mean) ** 2) / np.sum(w) + width = np.sqrt(var) + + # --- THE DEBUG LOGS --- + logging.info("=" * 40) + logging.info(f"IMAGE DEBUG SCORECARD") + logging.info(f"Final Spectrum Shape: {spectrum.shape}") + logging.info(f"Peak Location: Index {peak_idx}") + logging.info(f"Calculated SNR: {snr:.2f} (Threshold: 10.0)") + logging.info(f"Calculated Width: {width:.2f} (Range: 0.5 - 12.0)") + logging.info("=" * 40) + + logging.info(f"Cleaned Spectrum Shape for Function: {spectrum.shape}") + + # 3. Now pass the cleaned 1D array to your existing function + present = peak_is_present(spectrum, snr_threshold=10.0, width_range=(0.5, 12.0)) + + self.assertTrue(present, "Peak should have been detected in the cleaned spectrum.") + + @timeout(100) def test_cancel(self): """ diff --git a/src/odemis/gui/cont/tabs/sparc2_align_tab.py b/src/odemis/gui/cont/tabs/sparc2_align_tab.py index 17af22d75c..2b0d7e9f5c 100644 --- a/src/odemis/gui/cont/tabs/sparc2_align_tab.py +++ b/src/odemis/gui/cont/tabs/sparc2_align_tab.py @@ -45,6 +45,7 @@ Sparc2AutoFocus, Sparc2ManualFocus, ) +from odemis.acq.align.goffset_ext import auto_align_grating_detector_offsets from odemis.acq.stream_settings import StreamSettingsConfig from odemis.gui.comp import popup from odemis.gui.conf.data import get_hw_config, get_local_vas @@ -864,6 +865,9 @@ def add_axis(axisname, comp, label=None): # TODO: Auto remove the background when the image shape changes? # TODO: Use a toggle button to show the background is in use or not? + # Auto-calibration button + self.panel.btn_auto_calibrate.Bind(wx.EVT_BUTTON, self._on_btn_auto_calibrate) + def _on_btn_auto_align(self, evt): """ Handle the "Auto alignment" button click. @@ -1537,6 +1541,100 @@ def _on_align_mode_done(self, f): else: logging.debug("Optical path was updated.") + + # Auto-Calibration + def _on_btn_auto_calibrate(self, evt): + + # Check if there's a process running, if so, cancel and reset + if hasattr(self, "_auto_calibrate_future") and not self._auto_calibrate_future.done(): + self._auto_calibrate_future.cancel() + + self.panel.btn_auto_calibrate.SetLabel("Auto calibrate") + self.panel.gauge_auto_calibrate.SetValue(0) + return + + # Reset progress bar + self.panel.gauge_auto_calibrate.SetValue(0) + self.panel.btn_auto_calibrate.SetLabel("Cancel") + + wx.CallAfter(self._start_auto_calibration) + + self._progress_timer = wx.Timer(self.panel) + self.panel.Bind(wx.EVT_TIMER, self._on_progress_timer, self._progress_timer) + self._progress_timer.Start(200) + + def _start_auto_calibration(self): + main = self.tab_data_model.main + spectrograph = main.spectrograph + + # combine all detectors + detectors = getattr(main, "ccds", []) + getattr(main, "sp_ccds", []) + if not detectors: + detectors = [main.ccd] + + selector = getattr(main, "detector_selector", None) + + # If multiple detectors but no selector, only use the first detector + if len(detectors) > 1 and selector is None: + logging.warning( + "Multiple detectors detected but no selector; aligning only the first detector") + detectors = [detectors[0]] + + logging.info( + "Starting auto-calibration: detectors=%s selector=%s", + [d.name for d in detectors], selector.position.value if selector else None) + + # Start alignment procedure + self._auto_calibrate_future = auto_align_grating_detector_offsets( + spectrograph, detectors, selector=selector) + + # Bind progress & done callbacks + self._auto_calibrate_future.add_done_callback(self._on_auto_calibrate_done) + self.panel.btn_auto_calibrate.SetLabel("Cancel") + + def _on_progress_timer(self, evt): + if not hasattr(self, "_auto_calibrate_future"): + return + + f = self._auto_calibrate_future + + if f.done(): + self._progress_timer.Stop() + self.panel.gauge_auto_calibrate.SetValue(100) + return + + # Fetch where the backend logic says we are + target_val = int(getattr(f, "_progress", 0.0) * 100) + current_val = self.panel.gauge_auto_calibrate.GetValue() + + # Smooth UI trick: + # If we are behind the target, catch up quickly. + if current_val < target_val: + self.panel.gauge_auto_calibrate.SetValue(current_val + 1) + + # If we are caught up, STILL allow the bar to creep up to 5% ahead of the target. + # This ensures the bar instantly starts moving to ~4% while waiting + # for the hardware's 10-second stabilization time. + elif current_val < (target_val + 4) and current_val < 95: + self.panel.gauge_auto_calibrate.SetValue(current_val + 1) + + def _on_auto_calibrate_done(self, f): + try: + result = f.result() + logging.info("Auto calibration finished: %s", result) + except CancelledError: + logging.info("Auto calibration cancelled") + except Exception: + logging.exception("Auto calibration failed") + finally: + self._auto_calibrate_future = model.InstantaneousFuture() + + if hasattr(self, "_progress_timer"): + self._progress_timer.Stop() + + wx.CallAfter(self.panel.btn_auto_calibrate.SetLabel, "Auto calibrate") + wx.CallAfter(self.panel.gauge_auto_calibrate.SetValue, 0) + @call_in_wx_main def _on_lens_align_done(self, f): # Has no effect now, as OPM future are not cancellable (but it makes the From 39c44f1945e7c6f338db09eb543d082bc87b5a1b Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Thu, 7 May 2026 10:11:25 +0200 Subject: [PATCH 10/11] [fix] changed peak-finding method, added test to unit test --- src/odemis/acq/align/goffset.py | 277 ++++++++++--------- src/odemis/acq/align/test/goffset_test.py | 82 +++--- src/odemis/gui/cont/tabs/sparc2_align_tab.py | 20 +- 3 files changed, 206 insertions(+), 173 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index 6dce28e169..ba6d18df0e 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -10,10 +10,9 @@ from odemis.acq.align.autofocus import _mapDetectorToSelector from odemis.model import InstantaneousFuture from odemis.util import executeAsyncTask +from odemis.util.peak import Smooth, Detect from scipy.optimize import curve_fit from typing import Any, Dict, List, Optional, Tuple, Union -from scipy.ndimage import gaussian_filter1d - MOVE_TIME_GRATING = 20 # s MOVE_TIME_DETECTOR = 5 # s @@ -34,52 +33,43 @@ def gaussian(x: numpy.ndarray, amplitude: float, x0: float, width: float) -> num return intensity -def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: +def find_peak_position(data: numpy.ndarray, window_radius: int = 15, snr_threshold: float = 10.0) -> float: """ Finds the peak position in the given spectrum data. - It can handle both 1D and 2D data (in which case it averages over the first dimension). - - The function first identifies the absolute peak, then creates a window around it to minimize noise influence. - It then calculates a weighted average of the positions within this window, using the intensity values as weights. - For improved accuracy, it attempts to fit a Gaussian curve to the windowed data, estimating the peak's center, width - and amplitude. + It can handle both 1D and 2D data (in which case it uses Maximum Intensity Projection). - However, the curve fit can fail to converge or produce unreasonable results if the initial guess is poor, if there are - outliers, baseline trends, or if the data is too noisy. In such cases, the function falls back to using the weighted - average as the peak position estimate. The weighted average is more robust as it does not run an iterative optimizer, - so it has no convergence or numerical-optimization failure modes, but it may be less accurate if the peak is not well-defined - or if there are multiple peaks within the window. + The function first uses geometric shape-checking (Smooth and Detect) to accurately + find the peak index while ignoring high-amplitude 1-pixel noise spikes. + It then calculates a weighted average of the positions within a window, using the + intensity values as weights. For improved accuracy, it attempts to fit a Gaussian curve. - :param data: 1D or 2D array containing the spectrum (if 2D, it will be averaged to 1D) + :param data: 1D or 2D array containing the spectrum :param window_radius: number of pixels on either side of the peak to include in the window for fitting (default: 15) :return: estimated peak position in pixels (float) :raises RuntimeError: if no significant peak is detected (SNR too low) """ - if data.ndim == 2: - spectrum = data.mean(axis=0) # squash data into a 1D array - else: - spectrum = data + raw_data = numpy.asarray(data, dtype=float) + spectrum = numpy.squeeze(raw_data) - # compute SNR by comparing the peak height above the median background to the background magnitude and returns - # False if that SNR is below the threshold. This helps reject cases with no real peak, - # where the Gaussian fit would fail or produce nonsense results. + # maximum intensity projection + if spectrum.ndim > 1: + spectrum = spectrum.max(axis=0) - # peak_value = float(spectrum.max()) - # background = float(numpy.median(spectrum)) - # snr = (peak_value - background) / (abs(background) + 1e-6) + # geometric smoothing and noise calculation + smoothed_data = Smooth(spectrum, window_len=11) + clean_data = smoothed_data - numpy.median(smoothed_data) + noise_std = numpy.std(clean_data) - spectrum = spectrum - numpy.median(spectrum) - noise_std = numpy.std(spectrum) - peak_height = spectrum.max() - snr = peak_height / (noise_std + 1e-6) + # detect peak + maxtab, _ = Detect(smoothed_data, lookahead=10, delta=(noise_std + 1e-6) * snr_threshold) - if snr < 10.0: # tune threshold - raise RuntimeError("No peak detected (SNR too low)") + if not maxtab: + raise RuntimeError("No peak detected (SNR too low or peak geometry invalid)") - spectrum_smooth = gaussian_filter1d(spectrum, sigma=2) - peak_idx = int(numpy.argmax(spectrum_smooth)) - # peak_idx = numpy.argmax(spectrum) # find the absolute highest point + # sort peaks by intensity and grab the highest valid one + maxtab.sort(key=lambda x: x[1], reverse=True) + peak_idx = int(maxtab[0][0]) # create a window around the peak start = max(0, peak_idx - window_radius) @@ -94,15 +84,13 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: # as the best estimate of the peak, since weighted average would be undefined if weights.sum() == 0: weighted_avg = float(peak_idx) - logging.info("Weighted average fallback: all window data <= 0, using peak_idx=%d as estimate", - peak_idx) + logging.info("Weighted average fallback: all window data <= 0, using peak_idx=%d as estimate", peak_idx) else: weighted_avg = float(numpy.sum(window_idx * window_data) / numpy.sum(weights)) # try Gaussian fit for better accuracy, but fallback to weighted average if it fails or is out of bounds - try: - p0 = [window_data.max(), weighted_avg, 2.5] # intial guess: [amplitude, center, width] + p0 = [window_data.max(), weighted_avg, 2.5] # initial guess: [amplitude, center, width] popt, pcov = curve_fit(gaussian, window_idx, window_data, p0=p0) peak = popt[1] @@ -110,7 +98,7 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: return float(peak) except RuntimeError: - logging.info("Gaussian peak fit did not converge, falling back to weighted average") + logging.debug("Gaussian peak fit did not converge, falling back to weighted average") except ValueError: logging.exception("Gaussian peak fit failed due to invalid input") @@ -118,60 +106,69 @@ def find_peak_position(data: numpy.ndarray, window_radius: int = 15) -> float: return weighted_avg def peak_is_present(spectrum: numpy.ndarray, - snr_threshold: float=5.0, - width_range: Tuple[float, float]=(0.5, 12.0)) -> bool: + snr_threshold: float = 10.0, + width_range: Tuple[float, float] = (0.5, 12.0)) -> bool: """ Test to decide whether spectral peak is present. The test uses: - - a simple SNR threshold comparing the maximum to the median background, - - a local width estimate computed from a small window around the peak to - reject hot pixels and extremely broad features. + - The proven Odemis Detect() algorithm for geometric shape verification. + - A dynamic SNR threshold to reject background noise. + - A local width estimate computed from a small window around the peak. - :param spectrum: 1D array containing the spectrum (intensity vs pixel) - :param snr_threshold: minimum required signal-to-noise ratio for a peak to be considered present (default: 1) + :param spectrum: 1D or 2D array containing the spectrum + :param snr_threshold: minimum required signal-to-noise ratio for a peak (default: 10.0) :param width_range: acceptable range of estimated peak widths in pixels (default: (0.5, 12.0)) :return: True if a peak meeting the criteria is present, False otherwise """ - # basic stats - # peak_value = float(spectrum.max()) - # background = float(numpy.median(spectrum)) - # snr = (peak_value - background) / (abs(background) + 1e-6) + raw_data = numpy.asarray(spectrum, dtype=float) # convert to float + spectrum_1d = numpy.squeeze(raw_data) # remove any dummy dimensions - spectrum = spectrum - numpy.median(spectrum) - noise_std = numpy.std(spectrum) - peak_value = spectrum.max() - snr = peak_value / (noise_std + 1e-6) + # convert 2D data to 1D via MIP if it wasn't done by the caller + if spectrum_1d.ndim > 1: + spectrum_1d = spectrum_1d.max(axis=0) - if snr < snr_threshold: + smoothed_data = Smooth(spectrum_1d, window_len=11) # remove any hot pixels + clean_data = smoothed_data - numpy.median(smoothed_data) + noise_std = numpy.std(clean_data) + + # use Detect to reject high-amplitude noise spikes that lack Gaussian geometry + maxtab, _ = Detect(smoothed_data, lookahead=10, delta=(noise_std + 1e-6) * snr_threshold) + + if not maxtab: return False - # estimate width around the peak - spectrum_smooth = gaussian_filter1d(spectrum, sigma=2) - peak_idx = int(numpy.argmax(spectrum_smooth)) - # peak_idx = int(numpy.argmax(spectrum)) - if peak_idx < 1 or peak_idx > len(spectrum) - 2: - logging.debug("Peak too close to edge idx=%d len=%d", peak_idx, len(spectrum)) + # get the best peak found + maxtab.sort(key=lambda x: x[1], reverse=True) + peak_idx = int(maxtab[0][0]) + peak_value = maxtab[0][1] + + if peak_idx < 1 or peak_idx > len(spectrum_1d) - 2: + logging.debug("Peak too close to edge idx=%d len=%d", peak_idx, len(spectrum_1d)) return False - window = spectrum[peak_idx-2 : peak_idx+3] + # estimate width around the verified peak + window = spectrum_1d[peak_idx - 2: peak_idx + 3] # create 5-value window (2 on the left, 2 on the right of the peak) x = numpy.arange(len(window)) w = window - window.min() if w.sum() == 0: return False mean = numpy.sum(x * w) / numpy.sum(w) - var = numpy.sum(w * (x - mean)**2) / numpy.sum(w) + var = numpy.sum(w * (x - mean) ** 2) / numpy.sum(w) width = numpy.sqrt(var) + + snr = (peak_value - numpy.median(smoothed_data)) / (noise_std + 1e-6) present = width_range[0] <= width <= width_range[1] + logging.debug("snr=%.2f width=%.2f present=%s", snr, width, present) return present - -def coarse_scan_goffset_for_peak(spgr, detector, future: model.ProgressiveFuture, step: int=2000) -> Tuple[float, float]: - """ +def coarse_scan_goffset_for_peak(spgr, detector, future: model.ProgressiveFuture, + step: int = 2000) -> Tuple[float, float]: + """" Coarse scan across the goffset axis until a real peak becomes visible. The function performs absolute moves across the configured goffset axis @@ -197,7 +194,7 @@ def coarse_scan_goffset_for_peak(spgr, detector, future: model.ProgressiveFuture "Coarse local scan around goffset %.1f with step %.1f and max span %.1f", current, step, max_span) - # Build a sequence of positions: current, +step, -step, +2*step, -2*step, ... + # build a sequence of positions: current, +step, -step, +2*step, -2*step, ... positions = [current] k = 1 while k * step <= max_span: @@ -219,12 +216,13 @@ def coarse_scan_goffset_for_peak(spgr, detector, future: model.ProgressiveFuture logging.debug("Coarse scan move attempt to goffset %.3f", g) try: spgr.moveAbsSync({"goffset": g}) + time.sleep(2) # give hardware time to stabilize except ValueError: logging.warning("Skipping invalid goffset %.3f (%s)", g, ValueError) continue data = detector.data.get(asap=False) - spectrum = data.mean(axis=0) + spectrum = data.max(axis=0) if data.ndim == 2 else data if peak_is_present(spectrum): peak = find_peak_position(data) @@ -234,30 +232,30 @@ def coarse_scan_goffset_for_peak(spgr, detector, future: model.ProgressiveFuture def estimate_goffset_scale(spgr: model.Actuator, detector: model.Detector, - delta: float=5.0, - retries: int=1) -> Tuple[float, float, float]: + delta: float = 5.0, + retries: int = 1) -> Tuple[float, float, float]: """ - Estimate the scale factor between a change in the grating offset ('goffset') - and the resulting shift of the spectral peak on the detector. + Estimate the scale factor between a change in the grating offset ('goffset') + and the resulting shift of the spectral peak on the detector. - The function moves the actuator by a small step (delta) and measures the peak position before and after the move. - It calculates the ratio of pixel shift per unit of goffset. The actuator is returned to its original position - after measurement. + The function moves the actuator by a small step (delta) and measures the peak position before and after the move. + It calculates the ratio of pixel shift per unit of goffset. The actuator is returned to its original position + after measurement. - If the measured scale is unreasonably small or large, the function retries recursively - and falls back to a default value of 0.5 if necessary. + If the measured scale is unreasonably small or large, the function retries recursively + and falls back to a default value of 0.5 if necessary. - :param spgr: spectrograph - :param detector: detector - :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). - The actual step may be negated to avoid exceeding hardware limits. - :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). - - :return: Tuple (scale, p0, p1) - scale: estimated pixels per unit of goffset - p0: peak position at the initial goffset - p1: peak position after applying the test delta - """ + :param spgr: spectrograph + :param detector: detector + :param delta: The relative goffset step size to apply when measuring the scale (default: 5.0). + The actual step may be negated to avoid exceeding hardware limits. + :param retries: number of retries allowed if the estimated scale is unreliable (default: 1). + + :return: Tuple (scale, p0, p1) + scale: estimated pixels per unit of goffset + p0: peak position at the initial goffset + p1: peak position after applying the test delta + """ # get initial state data0 = detector.data.get(asap=False) @@ -292,7 +290,7 @@ def estimate_goffset_scale(spgr: model.Actuator, # If the retry still fails (raising a RuntimeError), we fall back to a default # scale value to ensure the algorithm can continue and avoid infinite recursion. - if abs(scale) < 1e-3 or abs(scale) > 10.0: + if abs(scale) < 1e-3 or abs(scale) > 2.0: logging.warning( "Unreliable scale estimate (%.4f). Retries left: %d", scale, retries) @@ -306,10 +304,20 @@ def estimate_goffset_scale(spgr: model.Actuator, return scale, p0, p1 def log_detector_state(caller: str, stage: str, detector: model.Detector, data: numpy.ndarray): + """ + Log a compact summary of the detector data: min, max, mean, background, noise, peak and SNR. + + :param caller: origin of the log message + :param stage: stage of the log message + :param detector: detector + :param data: data + :return: None + """ + shape = data.shape - # Keep consistent with detection pipeline - spectrum = data.mean(axis=0) if data.ndim == 2 else data + # convert 2D data to 1D if necessary + spectrum = data.max(axis=0) if data.ndim == 2 else data max_val = float(spectrum.max()) min_val = float(spectrum.min()) @@ -392,9 +400,9 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.info("Running alignment | detector=%s |", detector.name) try: - center_target = detector.resolution.value[0] / 2 + center_target = (detector.resolution.value[0] - 1) / 2 data0 = detector.data.get(asap=False) - spectrum0 = data0.mean(axis=0) if data0.ndim == 2 else data0 + spectrum0 = data0.max(axis=0) if data0.ndim == 2 else data0 if peak_is_present(spectrum0): peak0 = find_peak_position(data0) @@ -405,17 +413,6 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, peak_present = False peak0 = None - # # initial read: try to get a valid peak without moving the grating - # try: - # data0 = detector.data.get(asap=False) - # peak0 = find_peak_position(data0) # raises RuntimeError if no peak present - # logging.debug("Initial read: peak0=%.2f px", peak0) - # peak_present = True - # except RuntimeError: - # logging.debug("Initial read: no significant peak detected, scanning the axes") - # peak_present = False - # peak0 = None - # if peak present and already centered, do nothing if peak_present: initial_error_px = peak0 - center_target @@ -444,9 +441,6 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, return True # peak is present and not centered -> estimate scale and center - # scale, p0, p1 = estimate_goffset_scale(spgr, detector) - # logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) - try: scale, p0, p1 = estimate_goffset_scale(spgr, detector) logging.info("Scale estimated: %.4f px/goffset | p0=%.2f p1=%.2f", scale, p0, p1) @@ -471,14 +465,29 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, for i in range(max_it): _checkCancelled(future) + # if the peak is lost during the iteration, retry twice before aborting. if i == 0: peak_px = start_peak else: - data = detector.data.get(asap=False) - try: - peak_px = find_peak_position(data) - except RuntimeError: - logging.error("Peak re-acquisition failed during iteration %d — aborting", i) + max_retries = 2 + peak_found = False + + for attempt in range(max_retries + 1): + data = detector.data.get(asap=False) + try: + peak_px = find_peak_position(data) + peak_found = True + break + except RuntimeError: + if attempt < max_retries: + logging.warning("Peak lost in iteration %d. Settle and retry (%d/%d)...", + i, attempt + 1, max_retries) + time.sleep(2) + else: + pass + + if not peak_found: + logging.error(f"Peak re-acquisition permanently failed during iteration {i} — aborting") return False error_px = peak_px - center_target @@ -486,15 +495,15 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, if abs(error_px) <= tolerance_px: return True - # Calculate required adjustment based on pixel error and scaling factor + # calculate required adjustment based on pixel error and scaling factor delta_goffset = -gain * (error_px / scale) current = spgr.position.value["goffset"] - # Clamp move to stay within the allowed goffset range. - # Ensures that we don't command a step that would exceed axis limits. + # clamp move to stay within the allowed goffset range + # ensures that we don't command a step that would exceed axis limits. delta_goffset = max(minv - current, min(maxv - current, delta_goffset)) - # Accumulate the actual displacement for total movement tracking + # accumulate the actual displacement for total movement tracking total_goffset_displacement += delta_goffset logging.debug( @@ -503,6 +512,7 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, try: spgr.moveRelSync({"goffset": delta_goffset}) + time.sleep(0.2) # settle time for hardware to stabilize except ValueError: logging.warning("Hardware offset limit reached: %s. Keeping max allowed move.", ValueError) break @@ -519,7 +529,6 @@ def _do_sparc_auto_grating_offset(future: model.ProgressiveFuture, logging.error("Alignment error: %s", e) raise - def _cancel_sparc_auto_grating_offset(future: model.ProgressiveFuture) -> bool: """ Canceller of _do_sparc_auto_grating_offset task. @@ -559,19 +568,19 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, selector: Optional[model.Actuator] = None, streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: """ - Automatically align grating-detector offsets for all combinations of gratings and detectors. - - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector - will be used for all subsequent gratings. - - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will - be aligned by adjusting the detector offset with the grating alignment fixed. - - :param spectrograph: spectrograph - :param detectors: list of detectors - :param selector: optional selector to switch between detectors - :param streams: optional list of streams to update with progress - :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) - :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector - :raises CancelledError: if the operation is cancelled + Automatically align grating-detector offsets for all combinations of gratings and detectors. + - If a selector is provided, it will be used to switch between detectors for the first grating, then the first detector + will be used for all subsequent gratings. + - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will + be aligned by adjusting the detector offset with the grating alignment fixed. + + :param spectrograph: spectrograph + :param detectors: list of detectors + :param selector: optional selector to switch between detectors + :param streams: optional list of streams to update with progress + :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) + :raises ValueError: if no detectors provided, or if multiple detectors provided without a selector + :raises CancelledError: if the operation is cancelled """ if not isinstance(detectors, Iterable): @@ -642,12 +651,12 @@ def is_current_detector(d): return True return detector_to_selector[d] == selector.position.value[selector_axes] - # Calculate total steps for a simple progress bar + # calculate total steps for progress bar total_steps = len(detectors) + (len(gratings) - 1) current_step = 0 future._progress = 0.0 - # Start alignment for the first grating + # start alignment for the first grating try: g0 = gratings[0] logging.info("Starting alignment for initial grating: %s", g0) @@ -673,7 +682,7 @@ def is_current_detector(d): success = future._subfuture.result() results[(g0, d.name)] = success - # Update progress + # update progress current_step += 1 future._progress = current_step / total_steps @@ -696,7 +705,7 @@ def is_current_detector(d): success = future._subfuture.result() results[(g, first_detector.name)] = success - # Update progress + # update progress current_step += 1 future._progress = current_step / total_steps diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index 044e1769ba..1e2a694689 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -7,7 +7,6 @@ import unittest import logging import numpy as np -import odemis from odemis import model from odemis.util import timeout @@ -20,16 +19,9 @@ log_detector_state ) - - from odemis.dataio import hdf5 - - logging.getLogger().setLevel(logging.DEBUG) -CONFIG_PATH = os.path.dirname(odemis.__file__) + "/../../install/linux/usr/share/odemis/" -SPARC_CONFIG = CONFIG_PATH + "sim/sparc2-focus-test.odm.yaml" - HOME_PATH = os.path.expanduser("~") + "/" H5_FILE_2d_NO_PEAK = HOME_PATH + "development/odemis/grating 2 1024x256/" @@ -40,7 +32,6 @@ class TestSparcAutoGratingOffset(unittest.TestCase): @classmethod def setUpClass(cls): - #testing.start_backend(SPARC_CONFIG) cls.detector = model.getComponent(role="ccd") cls.spgr = model.getComponent(role="spectrograph") @@ -66,6 +57,7 @@ def test_find_peak_position_synthetic(self): """ Test peak detection on synthetic Gaussian data. """ + x = np.arange(200) true_center = 83.4 spectrum = np.exp(-0.5*((x-true_center)/3.0)**2) @@ -77,6 +69,7 @@ def test_find_peak_position_2d(self): """ Test peak detection on 2D data (mean over axis). """ + x = np.arange(200) true_center = 120.0 line = np.exp(-0.5*((x-true_center)/4.0)**2) @@ -184,7 +177,7 @@ def test_single_detector_updates_grating(self): start_goffset = self.spgr.position.value["goffset"] - f = sparc_auto_grating_offset(self.spgr, spccd, single_detector_mode=True, max_it=50) + f = sparc_auto_grating_offset(self.spgr, spccd, max_it=50) result = f.result(timeout=300) self.assertTrue(result, "Single-detector alignment failed") @@ -237,7 +230,7 @@ def test_multi_detector_detector_offset_changes(self): spccd = model.getComponent(role="sp-ccd") spccd.exposureTime.value = spccd.exposureTime.range[0] - # Ensure primary detector sets the grating first + # ensure primary detector sets the grating first try: self.selector.moveAbsSync({"rx": 0.0}) except Exception: @@ -247,13 +240,13 @@ def test_multi_detector_detector_offset_changes(self): self.assertTrue(f_first.result(timeout=300), "First-detector alignment failed") grating_after_first = self.spgr.position.value["goffset"] - # Switch to secondary detector + # switch to secondary detector try: self.selector.moveAbsSync({"rx": 1.5707963267948966}) except Exception: logging.debug("Selector move to secondary failed or not present; continuing") - # Measure peak before alignment on secondary detector + # measure peak before alignment on secondary detector data_before = spccd.data.get(asap=False) before_peak = float(find_peak_position(data_before)) @@ -261,16 +254,16 @@ def test_multi_detector_detector_offset_changes(self): f_second = sparc_auto_grating_offset(self.spgr, spccd, max_it=50) self.assertTrue(f_second.result(timeout=600), "Second-detector alignment failed") - # Measure peak after alignment + # measure peak after alignment data_after = spccd.data.get(asap=False) after_peak = float(find_peak_position(data_after)) - # Assert grating did not change + # assert grating did not change grating_after_second = self.spgr.position.value["goffset"] self.assertAlmostEqual(grating_after_first, grating_after_second, places=3, msg="Grating goffset changed when aligning second detector in multi-detector mode") - # Assert peak moved closer to center on the secondary detector + # assert peak moved closer to center on the secondary detector center = spccd.resolution.value[0] / 2 before_dist = abs(before_peak - center) after_dist = abs(after_peak - center) @@ -283,9 +276,7 @@ def test_single_detector_iteration(self): of offsets when operating with a single detector and multiple gratings. """ - #f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=[self.ccd], selector=self.selector) f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=self.detector, selector=self.selector) - res = f.result(timeout=900) n_gratings = len(self.spgr.axes["grating"].choices) @@ -299,11 +290,32 @@ def test_single_detector_iteration(self): self.assertEqual(len(dets_first), n_detectors) + def test_single_detector_alignment_algorithm(self): + + """ + Test automatic alignment of the gratings using a single detector. + """ + + f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=self.detector, selector=self.selector) + result = f.result(timeout=800) + + self.assertTrue(result) + def test_multi_detector_iteration(self): + """ + Tests the auto-alignment sequence for a multi-detector setup. + + Verifies that: + - Grating 1 aligns ALL detectors (calibrating the detector offsets). + - Subsequent gratings align ONLY the primary detector (calibrating the grating offsets). + - The total number of alignment runs is exactly as expected. + - The secondary detector receives valid optical data post-alignment. + """ + spccd = model.getComponent(role="sp-ccd") spccd.exposureTime.value = spccd.exposureTime.range[0] - detectors = [self.ccd, spccd] + detectors = [self.detector, spccd] # run alignment f = auto_align_grating_detector_offsets(spectrograph=self.spgr, detectors=detectors, selector=self.selector) @@ -337,17 +349,17 @@ def test_multi_detector_iteration(self): data = spccd.data.get(asap=False) # check data is not flat - self.assertNotEqual(data.max(), data.min) + self.assertNotEqual(data.max(), data.min()) def test_driver_raises_valueerror_on_hardware_limit(self): """ Verifies the driver successfully raises a ValueError when computing an offset outside of the hardware limits. """ - # A value guaranteed to exceed DET_OFFSET_MAX + GRAT_OFFSET_MAX + # out of bounds value out_of_bounds_target = 9999999.0 - # Test method directly to ensure the new bounds-checking logic works + # test method directly to ensure the new bounds-checking logic works with self.assertRaises(ValueError, msg="Driver should raise ValueError for out-of-bounds offset"): self.spgr._doSetGoffsetAbs(out_of_bounds_target) @@ -357,25 +369,25 @@ def test_hardware_limit_valueerror_handled(self): Verifies that when a goffset move exceeds hardware limits and throws a ValueError, the alignment algorithm catches it gracefully instead of crashing. """ - # Slightly misalign the peak so the algorithm wants to move + # slightly misalign the peak so the algorithm wants to move self.spgr.moveRelSync({"goffset": 50}) - # Temporarily bypass the software axis range clamp to ensure the raw command + # temporarily bypass the software axis range clamp to ensure the raw command # reaches the hardware driver to trigger the ValueError original_range = self.spgr.axes["goffset"].range try: self.spgr.axes["goffset"].range = (-9999999.0, 9999999.0) - # Execute auto-alignment procedure + # execute auto-alignment procedure f = sparc_auto_grating_offset(self.spgr, self.detector, max_it=10, gain=99999.0) result = f.result(timeout=100) - # Driver should raise ValueError for out-of-bounds move, which should be caught by the alignment code, + # driver should raise ValueError for out-of-bounds move, which should be caught by the alignment code self.assertFalse(result, "Alignment should gracefully fail (return False) when blocked by hardware limits.") finally: - # Restore original software limits and position + # restore original software limits and position self.spgr.axes["goffset"].range = original_range self.spgr.moveAbsSync({"goffset": self._original_position["goffset"]}) @@ -407,12 +419,16 @@ def test_find_peak_position_realistic_2d(self): image = np.array(image) - # run your actual function + # run function peak_pos = find_peak_position(image) self.assertAlmostEqual(peak_pos, true_center, places=1) def test_no_peak_present_image(self): + """ + Test peak detection on not-peak detector data. + """ + im_no_peak = hdf5.read_data(H5_FILE_2d_NO_PEAK + "no-peak-2d-1.h5") data = (im_no_peak[0].data) @@ -428,7 +444,7 @@ def test_peak_present_image(self): spectrum = np.squeeze(raw_data) if spectrum.ndim > 1: - # Average over the row dimension (256) to get the 1D spectrum (1024) + # convert 2D array into 1D by averaging spectrum = spectrum.mean(axis=0) clean_spec = spectrum - np.median(spectrum) @@ -437,7 +453,7 @@ def test_peak_present_image(self): snr = peak_val / (noise_std + 1e-6) peak_idx = np.argmax(clean_spec) - # Simple width estimation (standard deviation around peak) + # simple width estimation window = clean_spec[max(0, peak_idx - 2): min(len(clean_spec), peak_idx + 3)] x = np.arange(len(window)) w = window - window.min() @@ -447,7 +463,7 @@ def test_peak_present_image(self): var = np.sum(w * (x - mean) ** 2) / np.sum(w) width = np.sqrt(var) - # --- THE DEBUG LOGS --- + # debug logs logging.info("=" * 40) logging.info(f"IMAGE DEBUG SCORECARD") logging.info(f"Final Spectrum Shape: {spectrum.shape}") @@ -458,7 +474,7 @@ def test_peak_present_image(self): logging.info(f"Cleaned Spectrum Shape for Function: {spectrum.shape}") - # 3. Now pass the cleaned 1D array to your existing function + # pass the cleaned 1D array to the existing function present = peak_is_present(spectrum, snr_threshold=10.0, width_range=(0.5, 12.0)) self.assertTrue(present, "Peak should have been detected in the cleaned spectrum.") @@ -470,7 +486,7 @@ def test_cancel(self): Test cancelling alignment. """ f = sparc_auto_grating_offset(self.spgr, self.detector) - # Wait for the result or a timeout + # wait for the result or a timeout try: f.result(timeout=5) except: diff --git a/src/odemis/gui/cont/tabs/sparc2_align_tab.py b/src/odemis/gui/cont/tabs/sparc2_align_tab.py index 2b0d7e9f5c..3e3cdc1def 100644 --- a/src/odemis/gui/cont/tabs/sparc2_align_tab.py +++ b/src/odemis/gui/cont/tabs/sparc2_align_tab.py @@ -45,7 +45,7 @@ Sparc2AutoFocus, Sparc2ManualFocus, ) -from odemis.acq.align.goffset_ext import auto_align_grating_detector_offsets +from odemis.acq.align.goffset import auto_align_grating_detector_offsets from odemis.acq.stream_settings import StreamSettingsConfig from odemis.gui.comp import popup from odemis.gui.conf.data import get_hw_config, get_local_vas @@ -868,6 +868,13 @@ def add_axis(axisname, comp, label=None): # Auto-calibration button self.panel.btn_auto_calibrate.Bind(wx.EVT_BUTTON, self._on_btn_auto_calibrate) + # Auto-calibration state + self._auto_calibrate_future = model.InstantaneousFuture() + + # Create progress timer ONCE + self._progress_timer = wx.Timer(self.panel) + self.panel.Bind(wx.EVT_TIMER, self._on_progress_timer, self._progress_timer) + def _on_btn_auto_align(self, evt): """ Handle the "Auto alignment" button click. @@ -1549,6 +1556,9 @@ def _on_btn_auto_calibrate(self, evt): if hasattr(self, "_auto_calibrate_future") and not self._auto_calibrate_future.done(): self._auto_calibrate_future.cancel() + if self._progress_timer.IsRunning(): + self._progress_timer.Stop() + self.panel.btn_auto_calibrate.SetLabel("Auto calibrate") self.panel.gauge_auto_calibrate.SetValue(0) return @@ -1559,9 +1569,8 @@ def _on_btn_auto_calibrate(self, evt): wx.CallAfter(self._start_auto_calibration) - self._progress_timer = wx.Timer(self.panel) - self.panel.Bind(wx.EVT_TIMER, self._on_progress_timer, self._progress_timer) - self._progress_timer.Start(200) + if not self._progress_timer.IsRunning(): + self._progress_timer.Start(200) def _start_auto_calibration(self): main = self.tab_data_model.main @@ -1629,8 +1638,7 @@ def _on_auto_calibrate_done(self, f): finally: self._auto_calibrate_future = model.InstantaneousFuture() - if hasattr(self, "_progress_timer"): - self._progress_timer.Stop() + wx.CallAfter(self._progress_timer.Stop) wx.CallAfter(self.panel.btn_auto_calibrate.SetLabel, "Auto calibrate") wx.CallAfter(self.panel.gauge_auto_calibrate.SetValue, 0) From 450a0487a240b60c1875b7f1c00dd5b08e930aa8 Mon Sep 17 00:00:00 2001 From: Yu Xia de Jong Date: Tue, 12 May 2026 14:59:04 +0200 Subject: [PATCH 11/11] [feat] added turn on light + close slit for GUI --- src/odemis/acq/align/goffset.py | 55 ++++++++++++- src/odemis/acq/align/test/goffset_test.py | 86 +++++++++++++++++--- src/odemis/gui/cont/tabs/sparc2_align_tab.py | 49 ++++++++--- 3 files changed, 167 insertions(+), 23 deletions(-) diff --git a/src/odemis/acq/align/goffset.py b/src/odemis/acq/align/goffset.py index ba6d18df0e..97b1ab6cf2 100644 --- a/src/odemis/acq/align/goffset.py +++ b/src/odemis/acq/align/goffset.py @@ -6,6 +6,7 @@ from collections.abc import Iterable from concurrent.futures import CancelledError from concurrent.futures._base import CANCELLED, FINISHED, RUNNING +from odemis.acq.align import light from odemis import model from odemis.acq.align.autofocus import _mapDetectorToSelector from odemis.model import InstantaneousFuture @@ -565,6 +566,9 @@ def _total_alignment_time(n_gratings: int, def auto_align_grating_detector_offsets(spectrograph: model.Actuator, detectors: Union[model.Detector, List[model.Detector]], + opm: 'OpticalPathManager', + align_mode: str, + bl: model.Emitter, selector: Optional[model.Actuator] = None, streams: Optional[List['Stream']] = None) -> model.ProgressiveFuture: """ @@ -576,6 +580,9 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, :param spectrograph: spectrograph :param detectors: list of detectors + :param opm: OpticalPathManager + :param align_mode: optical alignment mode used to close slit/path + :param bl: brightlight emitter :param selector: optional selector to switch between detectors :param streams: optional list of streams to update with progress :return: ProgressiveFuture that will resolve to a dict mapping (grating, detector) @@ -596,7 +603,7 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, est_start = time.time() + 0.1 n_gratings = len(spectrograph.axes["grating"].choices) n_detectors = len(detectors) - a_time = _total_alignment_time(n_gratings, n_detectors) + a_time = (_total_alignment_time(n_gratings, n_detectors) + 10)# estimated time to turn on light and close slit f = model.ProgressiveFuture(start=est_start, end=est_start + a_time) f._progress = 0.0 f.task_canceller = _cancel_auto_align_grating_detector_offsets @@ -604,12 +611,15 @@ def auto_align_grating_detector_offsets(spectrograph: model.Actuator, f._task_lock = threading.Lock() f._task_state = RUNNING f._subfuture = InstantaneousFuture() - executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, selector, streams)) + executeAsyncTask(f, _do_auto_align_grating_detector_offsets, args=(f, spectrograph, detectors, opm, align_mode, bl, selector, streams)) return f def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, spectrograph: model.Actuator, detectors: List[model.Detector], + opm: 'OpticalPathManager', + align_mode: str, + bl: model.Emitter, selector: Optional[model.Actuator], streams: List['Stream'], stabilization_time: float = 10.0) -> Optional[Dict[Any, Any]]: @@ -619,10 +629,14 @@ def _do_auto_align_grating_detector_offsets(future: model.ProgressiveFuture, will be used for all subsequent gratings. - For multiple detectors, the grating alignment will only be adjusted for the first detector; subsequent detectors will be aligned by adjusting the detector offset with the grating alignment fixed. + - Before alignment, turn on the brightlight and closes the slit. :param future: ProgressiveFuture to update with progress and results :param spectrograph: spectrograph :param detectors: list of detectors + :param opm: OpticalPathManager + :param align_mode: optical alignment mode used to close slit/path + :param bl: brightlight emitter :param selector: optional selector to switch between detectors :param streams: optional list of streams to update with progress :param stabilization_time: time to wait after moving hardware before starting alignment (default: 15s) @@ -658,6 +672,37 @@ def is_current_detector(d): # start alignment for the first grating try: + _checkCancelled(future) + logging.info("Preparing optical path and turning on light") + + # make sure it's in 0th order (ie, show the image as-is), this also ensures the spectrograph + # is done with any previous actions. + spectrograph.moveAbsSync({"wavelength": 0}) + + # logging.info("Forcing brightlight ON") # bypassing sensor check for simulation + # bl.power.value = bl.power.range[1] + # time.sleep(1) + + logging.info("Turning on brightlight") + + future._subfuture = light.turnOnLight(bl, first_detector) + + try: + future._subfuture.result(timeout=60) + + except TimeoutError: + future._subfuture.cancel() + logging.warning("Brightlight did not confirm ON within 60s; continuing anyway") + + _checkCancelled(future) + + logging.info("Setting optical path to alignment mode: %s",align_mode) + future._subfuture = opm.setPath(align_mode,detector=first_detector) + future._subfuture.result() + + _checkCancelled(future) + future._subfuture = InstantaneousFuture() + g0 = gratings[0] logging.info("Starting alignment for initial grating: %s", g0) @@ -719,6 +764,12 @@ def is_current_detector(d): raise finally: + logging.info("Turning off brightlight") + try: + bl.power.value = (bl.power.range[0]) + except Exception: + logging.exception("Failed to turn off the light during alignment cleanup") + spectrograph.moveAbsSync(original_pos) if selector: selector.moveAbsSync(original_selector) diff --git a/src/odemis/acq/align/test/goffset_test.py b/src/odemis/acq/align/test/goffset_test.py index 1e2a694689..6f3b545feb 100644 --- a/src/odemis/acq/align/test/goffset_test.py +++ b/src/odemis/acq/align/test/goffset_test.py @@ -3,23 +3,23 @@ Test Sparc auto grating offset alignment """ -import os -import unittest import logging import numpy as np +import os +import threading +import unittest +from concurrent.futures import Future +from concurrent.futures._base import RUNNING from odemis import model from odemis.util import timeout -from odemis.acq.align.goffset_ext import( - find_peak_position, - peak_is_present, - estimate_goffset_scale, - sparc_auto_grating_offset, - auto_align_grating_detector_offsets, - log_detector_state -) - +from odemis.acq.align.goffset import(find_peak_position, peak_is_present, estimate_goffset_scale, + sparc_auto_grating_offset, auto_align_grating_detector_offsets, + _do_auto_align_grating_detector_offsets, log_detector_state) from odemis.dataio import hdf5 +from odemis.model import ProgressiveFuture +from unittest.mock import patch, MagicMock + logging.getLogger().setLevel(logging.DEBUG) HOME_PATH = os.path.expanduser("~") + "/" @@ -479,6 +479,70 @@ def test_peak_present_image(self): self.assertTrue(present, "Peak should have been detected in the cleaned spectrum.") + @patch('odemis.acq.align.goffset.sparc_auto_grating_offset') + @patch('odemis.acq.align.goffset.light.turnOnLight') + @patch('odemis.acq.align.goffset.time.sleep') + + def test_do_auto_align_mocked(self, mock_sleep, mock_turnOnLight, mock_sparc_offset): + """ + Test the full auto-calibration sequence using mocked hardware to bypass + timeouts, sleep delays, and physical hardware requirements. + """ + + # set up mocked return values using standard python futures + fake_light_future = Future() + fake_light_future.set_result(True) + mock_turnOnLight.return_value = fake_light_future + + fake_align_future = Future() + fake_align_future.set_result(True) + mock_sparc_offset.return_value = fake_align_future + + # set up fake hardware + fake_spectrograph = MagicMock() + + # Mock the initial position and the grating choices + fake_spectrograph.position.value = {"wavelength": 0, "grating": "grating1"} + fake_spectrograph.axes = {"grating": MagicMock(choices={"grating1": 0, "grating2": 1})} + + fake_detector = MagicMock() + fake_detector.name = "fake_ccd" + fake_detector.data.get.return_value = np.zeros((1024,)) + + fake_opm = MagicMock() + fake_path_future = Future() + fake_path_future.set_result(True) + fake_opm.setPath.return_value = fake_path_future + + fake_bl = MagicMock() + fake_bl.power.range = [0, 100] + + main_future = ProgressiveFuture(start=0, end=100) + main_future._task_lock = threading.Lock() + main_future._task_state = RUNNING + main_future._subfuture = Future() + + # run function + results = _do_auto_align_grating_detector_offsets(future=main_future, spectrograph=fake_spectrograph, + detectors=[fake_detector], opm=fake_opm, + align_mode="spec-focus", bl=fake_bl, selector=None, + streams=[]) + + # Assertions + # did it try to turn on the light? + mock_turnOnLight.assert_called_once_with(fake_bl, fake_detector) + + # did it set the optical path? + fake_opm.setPath.assert_called_once_with("spec-focus", detector=fake_detector) + + # did it try to align both gratings we mocked? + self.assertEqual(mock_sparc_offset.call_count, 2) + + # did it clean up and turn the light off in the finally block? + self.assertEqual(fake_bl.power.value, 0) + + # check the final returned results dict + self.assertEqual(results, {("grating1", "fake_ccd"): True, ("grating2", "fake_ccd"): True}) @timeout(100) def test_cancel(self): diff --git a/src/odemis/gui/cont/tabs/sparc2_align_tab.py b/src/odemis/gui/cont/tabs/sparc2_align_tab.py index 3e3cdc1def..95e8e87a6b 100644 --- a/src/odemis/gui/cont/tabs/sparc2_align_tab.py +++ b/src/odemis/gui/cont/tabs/sparc2_align_tab.py @@ -1574,28 +1574,57 @@ def _on_btn_auto_calibrate(self, evt): def _start_auto_calibration(self): main = self.tab_data_model.main - spectrograph = main.spectrograph + align_mode = self.tab_data_model.align_mode.value + opm = main.opm + + # Set the optical path according to the align mode + if align_mode == "streak-align": + if (main.streak_ccd + and main.spectrograph_ded + and main.streak_ccd.name in main.spectrograph_ded.affects.value): + opath = "streak-focus-ext" + else: + opath = "streak-focus" + elif align_mode == "tunnel-lens-align": + opath = "spec-focus-ext" + elif align_mode in ("lens-align", "lens2-align", "light-in-align"): + opath = "spec-focus" + else: + logging.warning("Auto calibration requested not compatible with requested alignment mode %s. Do nothing.", + align_mode) + return - # combine all detectors + # Pick the right hardware based on whether opath is external or internal + if opath in ("spec-focus-ext", "streak-focus-ext"): + bl = main.brightlight_ext + spectrograph = main.spectrograph_ded + selector = getattr(main, "spec_ded_det_selector", None) + else: + bl = main.brightlight + spectrograph = main.spectrograph + selector = getattr(main, "spec_det_selector", None) + + if not bl or not spectrograph: + logging.error("Missing brightlight or spectrograph. Cannot start auto-calibration.") + return + + # Combine all detectors detectors = getattr(main, "ccds", []) + getattr(main, "sp_ccds", []) if not detectors: - detectors = [main.ccd] - - selector = getattr(main, "detector_selector", None) + detectors = [main.ccd] if main.ccd else [] # If multiple detectors but no selector, only use the first detector if len(detectors) > 1 and selector is None: - logging.warning( - "Multiple detectors detected but no selector; aligning only the first detector") + logging.warning("Multiple detectors detected but no selector; aligning only the first detector") detectors = [detectors[0]] logging.info( - "Starting auto-calibration: detectors=%s selector=%s", - [d.name for d in detectors], selector.position.value if selector else None) + "Starting auto-calibration: detectors=%s selector=%s path=%s", + [d.name for d in detectors], selector.position.value if selector else None, opath) # Start alignment procedure self._auto_calibrate_future = auto_align_grating_detector_offsets( - spectrograph, detectors, selector=selector) + spectrograph, detectors, opm, opath, bl, selector=selector) # Bind progress & done callbacks self._auto_calibrate_future.add_done_callback(self._on_auto_calibrate_done)