[perf] Improve spot finding in estimate_grid_orientation_from_img

install/linux/usr/share/odemis/sim/fastem-sim-asm.odm.yaml

@@ -335,6 +335,7 @@ FASTEM-sim: {
             "percntl_im_rng": [1, 99],
             # List[int, int] min/max percentile of the dynamic range to which the histogram of the image will be stretched.
             "percntl_dyn_rng": [10, 90],
+            "pitch": 3.2e-6  # float [m], the expected pitch between two beamlets in meters
         },
     },
 }
@@ -366,7 +367,7 @@ FASTEM-sim: {
     class: static.OpticalLens,
     role: lens,
     init: {
-        mag: 60, # ratio, magnifying; higher magnification is a stronger simulated blur
+        mag: 40, # ratio, magnifying; higher magnification is a stronger simulated blur
         na: 0.95,  # numerical aperture
         ri: 1,  # refractive index
     },

scripts/spot-grid.py

@@ -20,6 +20,7 @@
 import wx.lib.wxcairo
 
 from odemis import dataio, model
+from odemis.acq.fastem import DEFAULT_PITCH
 from odemis.cli.video_displayer import VideoDisplayer
 from odemis.driver import ueye
 from odemis.gui.conf.file import AcquisitionConfig
@@ -35,6 +36,8 @@
 PIXEL_SIZE_SAMPLE_PLANE = 3.45e-6  # m
 DEFAULT_MAGNIFICATION = 40
 PIXEL_SIZE = PIXEL_SIZE_SAMPLE_PLANE / DEFAULT_MAGNIFICATION
+# 0.75 is a safety factor to allow for some variation in spot positions
+MIN_DIST_SPOTS = int(0.75 * DEFAULT_PITCH / PIXEL_SIZE)
 
 
 class VideoDisplayerGrid(VideoDisplayer):
@@ -43,16 +46,19 @@ class VideoDisplayerGrid(VideoDisplayer):
     It should be pretty much platform independent.
     """
 
-    def __init__(self, title="Live image", size=(640, 480), gridsize=None, pixel_size=PIXEL_SIZE):
+    def __init__(self, title="Live image", size=(640, 480), gridsize=None, pixel_size=PIXEL_SIZE,
+                 min_dist_spots=MIN_DIST_SPOTS):
         """
         Displays the window on the screen
         size (2-tuple int,int): X and Y size of the window at initialisation
         pixel_size (float): pixel size in m
+        min_dist_spots (int): minimum distance between spots in pixels
         Note that the size of the window automatically adapts afterwards to the
         coming pictures
         """
         self.app = ImageWindowApp(title, size, pixel_size)
         self.gridsize = (8, 8) if gridsize is None else gridsize
+        self.min_dist_spots = min_dist_spots
         self.acqui_conf = AcquisitionConfig()
 
     def new_image(self, data):
@@ -68,6 +74,7 @@ def new_image(self, data):
                 AffineTransform,
                 sigma=1.45,
                 threshold_rel=self.acqui_conf.spot_grid_threshold,
+                min_distance=self.min_dist_spots,
             )
             grid = unit_gridpoints(self.gridsize, mode="ji")
             self.app.spots = tform_ji.apply(grid)
@@ -244,17 +251,19 @@ def start_generate(self):
         self.notify(self._detector.array)
 
 
-def live_display(ccd, dataflow, pixel_size, kill_ccd=True, gridsize=None):
+def live_display(ccd, dataflow, pixel_size, kill_ccd=True, gridsize=None, min_dist_spots=MIN_DIST_SPOTS):
     """
     Acquire an image from one (or more) dataflow and display it with a spot grid overlay.
     ccd: a camera object
     dataflow: dataflow to access
     pixel_size (float): pixel size in m
     kill_ccd: True if it is required to terminate the ccd after closing the window
     gridsize: size of the grid of spots.
+    min_dist_spots: minimum distance between spots in pixels
     """
     # create a window
-    window = VideoDisplayerGrid("Live from %s.%s" % (ccd.role, "data"), ccd.resolution.value, gridsize, pixel_size)
+    window = VideoDisplayerGrid("Live from %s.%s" % (ccd.role, "data"), ccd.resolution.value, gridsize, pixel_size,
+                                min_dist_spots)
     im_passer = ImagePasser()
     t = threading.Thread(target=image_update, args=(im_passer, window))
     t.daemon = True
@@ -293,6 +302,8 @@ def main(args):
                         help="size of the grid of spots in x y, default 8 8")
     parser.add_argument("--magnification", dest="magnification", type=float,
                         help="magnification (typically 40 or 50)")
+    parser.add_argument("--pitch", dest="pitch", type=float, default=None,
+                        help=f"pitch in meters (defaults to  {DEFAULT_PITCH:0.1e})")
     parser.add_argument("--log-level", dest="loglev", metavar="<level>", type=int, choices=[0, 1, 2],
                         default=0, help="set verbosity level (0-2, default = 0)")
     options = parser.parse_args(args[1:])
@@ -327,21 +338,35 @@ def main(args):
         logging.warning("No magnification specified, falling back to %s.", magnification)
     pixel_size = PIXEL_SIZE_SAMPLE_PLANE / magnification
 
+    pitch = options.pitch
+    if not pitch:
+        try:
+            mppc = model.getComponent(role="mppc")
+            mppc_md = mppc.getMetadata()
+            pitch = mppc_md.get(model.MD_CALIB, {}).get("pitch", DEFAULT_PITCH)
+        except Exception as ex:
+            logging.debug("Failed to read pitch from mppc, ex: %s", ex)
+            pitch = DEFAULT_PITCH
+
+    # 0.75 is a safety factor to allow for some variation in spot positions
+    min_dist_spots = int(0.75 * pitch / pixel_size)
+
     if options.filename:
         logging.info("Will process image file %s" % options.filename)
         converter = dataio.find_fittest_converter(options.filename, default=None, mode=os.O_RDONLY)
         data = converter.read_data(options.filename)[0]
         fakeccd = StaticCCD(options.filename, "fakeccd", data)
-        live_display(fakeccd, fakeccd.data, pixel_size, gridsize=options.gridsize)
+        live_display(fakeccd, fakeccd.data, pixel_size, gridsize=options.gridsize, min_dist_spots=min_dist_spots)
     elif options.role:
         if get_backend_status() != BACKEND_RUNNING:
             raise ValueError("Backend is not running while role command is specified.")
         ccd = model.getComponent(role=options.role)
-        live_display(ccd, ccd.data, pixel_size, kill_ccd=False, gridsize=options.gridsize)
+        live_display(ccd, ccd.data, pixel_size, kill_ccd=False, gridsize=options.gridsize,
+                     min_dist_spots=min_dist_spots)
     else:
         ccd = ueye.Camera("camera", "ccd", device=None)
         ccd.SetFrameRate(2)
-        live_display(ccd, ccd.data, pixel_size, gridsize=options.gridsize)
+        live_display(ccd, ccd.data, pixel_size, gridsize=options.gridsize, min_dist_spots=min_dist_spots)
     return 0
 
 

src/odemis/acq/fastem.py

@@ -53,6 +53,8 @@
 # The executor is a single object, independent of how many times the module (fastem.py) is loaded.
 _executor = model.CancellableThreadPoolExecutor(max_workers=1)
 
+DEFAULT_PITCH = 3.2e-6  # distance between spots in m
+
 # TODO: Normally we do not use component names in code, only roles. Store in the roles in the SETTINGS_SELECTION,
 #  and at init lookup the role -> name conversion (using model.getComponent(role=role)).
 # Selection of components, VAs and values to save with the ROA acquisition, structured: {component: {VA: value}}
@@ -292,6 +294,16 @@ def __init__(self, scanner, multibeam, descanner, detector, stage, scan_stage, c
         # save the initial multibeam resolution, because the resolution will get updated if save_full_cells is True
         self._old_res = self._multibeam.resolution.value
 
+        # Calculate the expected minimum distance between spots in the grid on the diagnostic camera
+        detector_md = detector.getMetadata()
+        ccd_md = ccd.getMetadata()
+        self._exp_pitch_m = detector_md.get(model.MD_CALIB, {}).get("pitch", DEFAULT_PITCH)  # m
+        lens_mag = ccd_md.get(model.MD_LENS_MAG)
+        ccd_px_size = ccd_md.get(model.MD_SENSOR_PIXEL_SIZE)
+        exp_pitch_px = self._exp_pitch_m * lens_mag / ccd_px_size[0]
+        # 0.75 is a safety factor to allow for some variation in spot positions
+        self._min_dist_spots = int(0.75 * exp_pitch_px)
+
         beam_shift_path = fastem_util.create_image_dir("beam-shift-correction")
         # If there is a project name the path will be
         # [image-dir]/beam-shift-correction/[timestamp]/[project-name]/[roa-name]_[slice-idx]
@@ -457,8 +469,8 @@ def acquire_roa(self, dataflow):
                 logging.debug(f"Will run beam shift correction for field index {field_idx}")
                 try:
                     new_beam_shift = self.correct_beam_shift()
-                    # The difference in x or y should not be larger than 2 micrometers
-                    if any(map(lambda n, p: abs(n - p) > 2e-6, new_beam_shift, prev_beam_shift)):
+                    # The difference in x or y should not be larger than half a pitch
+                    if any(map(lambda n, p: abs(n - p) > 0.5 * self._exp_pitch_m, new_beam_shift, prev_beam_shift)):
                         raise ValueError(
                             f"Difference in beam shift is larger than 2 µm, therefore it most likely failed. "
                             f"Previous beam shift: {prev_beam_shift}, new beam shift: {new_beam_shift}"
@@ -663,7 +675,9 @@ def correct_beam_shift(self):
         # asap=False: wait until new image is acquired (don't read from buffer)
         ccd_image = self._ccd.data.get(asap=False)
         tform, error = estimate_grid_orientation_from_img(ccd_image, (8, 8), SimilarityTransform, sigma,
-                                                          threshold_rel=self._spot_grid_thresh)
+                                                          threshold_rel=self._spot_grid_thresh,
+                                                          min_distance=self._min_dist_spots,
+                                                          )
         logging.debug(f"Found center of grid at {tform.translation}, error: {error}.")
 
         # Determine the shift of the spots, by subtracting the good multiprobe position from the average (center)

src/odemis/acq/test/fastem_test.py

@@ -37,7 +37,7 @@
 from odemis.acq import fastem, stream
 from odemis.acq.acqmng import SettingsObserver
 from odemis.acq.align.fastem import Calibrations
-from odemis.acq.fastem import SETTINGS_SELECTION
+from odemis.acq.fastem import DEFAULT_PITCH, SETTINGS_SELECTION
 from odemis.gui.comp.fastem_roa import FastEMROA
 from odemis.gui.comp.overlay.shapes import EditableShape
 from odemis.gui.model.main_gui_data import FastEMMainGUIData
@@ -757,10 +757,10 @@ def test_get_abs_stage_movement(self):
             # Create an ROA with the coordinates of the field.
             roa_name = "test_megafield_id"
             roa = FastEMROA(shape=MockEditableShape(),
-                        main_data=self.main_data,
-                        overlap=0.0,
-                        name=roa_name,
-                        slice_index=0)
+                            main_data=self.main_data,
+                            overlap=0.0,
+                            name=roa_name,
+                            slice_index=0)
             roa.shape._points = points
             roa.shape.points.value = points
 
@@ -845,10 +845,10 @@ def test_get_abs_stage_movement_overlap(self):
             # Create an ROA with the coordinates of the field.
             roa_name = "test_megafield_id"
             roa = FastEMROA(shape=MockEditableShape(),
-                        main_data=self.main_data,
-                        overlap=overlap,
-                        name=roa_name,
-                        slice_index=0)
+                            main_data=self.main_data,
+                            overlap=overlap,
+                            name=roa_name,
+                            slice_index=0)
             roa.shape._points = points
             roa.shape.points.value = points
 
@@ -1130,10 +1130,10 @@ def test_get_pos_first_tile(self):
             # Create an ROA with the coordinates of the field.
             roa_name = "test_megafield_id"
             roa = FastEMROA(shape=MockEditableShape(),
-                        main_data=self.main_data,
-                        overlap=0.0,
-                        name=roa_name,
-                        slice_index=0)
+                            main_data=self.main_data,
+                            overlap=0.0,
+                            name=roa_name,
+                            slice_index=0)
             roa.shape._points = points
             roa.shape.points.value = points
 
@@ -1295,6 +1295,9 @@ def setUpClass(cls):
         cls.mppc.frameDuration.value = 0.1
         cls.mppc.cellCompleteResolution.value = (900, 900)
         cls.mppc.shape = (8, 8)
+        cls.mppc.configure_mock(
+            **{"getMetadata.return_value": {model.MD_CALIB: {"pitch": DEFAULT_PITCH},}}
+        )
 
         cls.multibeam = Mock()
         cls.multibeam.configure_mock(**{"getMetadata.return_value": {model.MD_SCAN_OFFSET_CALIB: [0.01, 0.01],
@@ -1328,7 +1331,11 @@ def setUpClass(cls):
         image[54:150:12, 54:150:12] = 1
         image = model.DataArray(input_array=image)
         cls.ccd.data.configure_mock(**{"get.return_value": image})
-        cls.ccd.configure_mock(**{"getMetadata.return_value": {model.MD_FAV_POS_ACTIVE: {"j": 100, "i": 100}}})
+        cls.ccd.configure_mock(**{"getMetadata.return_value": {
+            model.MD_FAV_POS_ACTIVE: {"j": 100, "i": 100},
+            model.MD_SENSOR_PIXEL_SIZE: (3.45e-6, 3.45e-6),
+            model.MD_LENS_MAG: 40,
+        }})
         cls.ccd.pointSpreadFunctionSize.value = 1
         cls.ccd.pixelSize.value = (1.0e-7, 1.0e-7)
 

src/odemis/util/peak_local_max.py

@@ -364,6 +364,7 @@ def peak_local_max(
     num_peaks: Optional[int] = None,
     footprint: Optional[np.ndarray] = None,
     p_norm: float = np.inf,
+    len_object: Optional[int] = None,
 ) -> np.ndarray:
     """
     Find peaks in an image as coordinate list.
@@ -378,7 +379,7 @@ def peak_local_max(
     image : ndarray
         Input image.
     min_distance : int, optional
-        The minimal allowed distance separating peaks. To find the
+        The minimal allowed distance in pixels separating peaks. To find the
         maximum number of peaks, use `min_distance=1`.
     threshold_abs : float, optional
         Minimum intensity of peaks. By default, the absolute threshold is
@@ -406,6 +407,12 @@ def peak_local_max(
         A finite large p may cause a ValueError if overflow can occur.
         ``inf`` corresponds to the Chebyshev distance and 2 to the
         Euclidean distance.
+    len_object : int, optional
+        The length of the object in pixels. This parameter is used to determine
+        the border width for peak exclusion when `exclude_border` is True. If
+        not provided, `min_distance` is used as the length of the object. This
+        can be useful when the object size is known and different from
+        `min_distance`.
 
     Returns
     -------
@@ -460,7 +467,10 @@ def peak_local_max(
             stacklevel=2,
         )
 
-    border_width = _get_excluded_border_width(image, min_distance, exclude_border)
+    if not len_object:
+        len_object = min_distance
+
+    border_width = _get_excluded_border_width(image, len_object, exclude_border)
 
     threshold = _get_threshold(image, threshold_abs, threshold_rel)
 

src/odemis/util/registration.py

@@ -532,6 +532,7 @@ def estimate_grid_orientation_from_img(
     threshold_abs: Optional[float] = None,
     threshold_rel: Optional[float] = None,
     num_spots: Optional[int] = None,
+    min_distance: Optional[int] = None,
 ) -> Tuple[T, float]:
     """
     Image based estimation of the orientation of a square grid of points.
@@ -559,6 +560,9 @@ def estimate_grid_orientation_from_img(
         `num_spots = shape[0] * shape[1]` as default when set to `None`. Set
         `num_spots = 0` to not impose a maximum. Note that this behavior is
         different from odemis.util.spot.find_spot_position().
+    min_distance : int, optional
+        The minimal allowed distance in pixels separating peaks. To find the
+        maximum number of peaks, use `min_distance=1`.
 
     Returns
     -------
@@ -578,5 +582,5 @@ def estimate_grid_orientation_from_img(
         num_spots = shape[0] * shape[1]
     elif num_spots == 0:
         num_spots = None
-    ji = find_spot_positions(image, sigma, threshold_abs, threshold_rel, num_spots)
+    ji = find_spot_positions(image, sigma, threshold_abs, threshold_rel, num_spots, min_distance)
     return estimate_grid_orientation(ji, shape, transform_type)

src/odemis/util/spot.py

@@ -369,6 +369,7 @@ def find_spot_positions(
     threshold_abs: Optional[float] = None,
     threshold_rel: Optional[float] = None,
     num_spots: Optional[int] = None,
+    min_distance: Optional[int] = None,
 ) -> numpy.ndarray:
     """
     Find the center coordinates of spots with the highest intensity in an
@@ -390,6 +391,9 @@ def find_spot_positions(
     num_spots : int, optional
         Maximum number of spots. When the number of spots exceeds `num_spots`,
         return `num_spots` peaks based on highest spot intensity.
+    min_distance : int, optional
+        The minimal allowed distance in pixels separating peaks. To find the
+        maximum number of peaks, use `min_distance=1`.
 
     Returns
     -------
@@ -399,15 +403,18 @@ def find_spot_positions(
 
     """
     size = int(round(3 * sigma))
-    min_distance = 2 * size
-    filtered = bandpass_filter(image, sigma, min_distance)
+    len_object = 2 * size  # typical length of a spot
+    min_distance = len_object if not min_distance else min_distance  # distance between spots
+    filtered = bandpass_filter(image, sigma, len_object)
     coordinates = peak_local_max(
         filtered,
         min_distance=min_distance,
         threshold_abs=threshold_abs,
         threshold_rel=threshold_rel,
+        exclude_border=False,
         num_peaks=num_spots,
         p_norm=2,
+        len_object=len_object,
     )
 
     # Improve coordinate estimate using radial symmetry center.

src/odemis/util/test/spot_test.py

@@ -339,6 +339,26 @@ def test_multiple(self):
         numpy.testing.assert_array_equal(sorted(indices), range(len(loc)))
         numpy.testing.assert_array_less(distances, 0.05)
 
+    def test_spots_close_to_edge(self):
+        """
+        `find_spot_positions` should find all spot positions in a generated
+        test image when the minimum distance between spots is larger than the
+        distance between the spots and the edge of the image.
+        """
+        # set a grid of 8 by 8 points to 1 at the top left of the image
+        image = numpy.zeros((256, 256))
+        image[4:100:12, 8:104:12] = 1
+        expected_ji = numpy.column_stack(numpy.where(image))
+
+        ji = spot.find_spot_positions(image, sigma=0.75, min_distance=12)
+
+        # Sort both arrays to ensure consistent ordering
+        expected_ji_sorted = expected_ji[numpy.lexsort((expected_ji[:, 1], expected_ji[:, 0]))]
+        ji_sorted = ji[numpy.lexsort((ji[:, 1], ji[:, 0]))]
+
+        # Check if the sorted arrays are equal
+        numpy.testing.assert_array_almost_equal(expected_ji_sorted, ji_sorted)
+
 
 if __name__ == "__main__":
     unittest.main()