liberate cellfinder from the concrete

Author: i-jey

ulc_mm_package/hardware/scope_routines.py

@@ -520,13 +520,6 @@ def find_cells_routine(
 
         img = yield
 
-        # Initial check for cells, return current motor position if cells found
-        cell_finder.add_image(mscope.motor.pos, img)
-        try:
-            return cell_finder.get_cells_found_position()
-        except NoCellsFound:
-            cell_finder.reset()
-
         # Defensive check, ensure the motor isn't moving (say for example,
         # if CellFinder was triggered by an OOF exception and SSAF just triggered a motor move)
         while mscope.motor.is_locked():
@@ -582,7 +575,7 @@ def find_cells_routine(
                 )
 
             self.logger.info("Looking for cells...")
-            # Perform a full focal stack and get the cross-correlation value for each image
+
             # If we're currently at the bottom, do the bottom-up sweep. Otherwise, do the top-down sweep.
             if mscope.motor.pos == 0:
                 for pos in range(0, mscope.motor.max_pos, steps_per_image):
@@ -603,18 +596,11 @@ def find_cells_routine(
                     except NoCellsFound:
                         pass
             else:
-                # Move from the current position to the bottom sweep as we're going down
-                for pos in range(mscope.motor.pos, 0, -steps_per_image):
-                    mscope.motor.move_abs(pos)
-                    img = yield
-                    cell_finder.add_image(mscope.motor.pos, img)
-                    try:
-                        return cell_finder.get_cells_found_position()
-                    except NoCellsFound:
-                        pass
-
-                # If cells not found on the way down, sweep all the way back up
-                for pos in range(0, mscope.motor.max_pos, steps_per_image):
+                # Sweep down from current position to bottom, then sweep from bottom to top
+                positions = list(range(mscope.motor.pos, 0, -steps_per_image)) + list(
+                    range(0, mscope.motor.max_pos, steps_per_image)
+                )
+                for pos in positions:
                     mscope.motor.move_abs(pos)
                     img = yield
                     cell_finder.add_image(mscope.motor.pos, img)

ulc_mm_package/image_processing/cell_finder.py

@@ -5,50 +5,6 @@
 
 
 from ulc_mm_package.image_processing.focus_metrics import downsample_image
-from ulc_mm_package.image_processing.processing_constants import (
-    RBC_THUMBNAIL_PATH,
-    CELLS_FOUND_THRESHOLD,
-    MIN_POINTS_ABOVE_THRESH,
-)
-
-RBC_THUMBNAIL = cv2.imread(RBC_THUMBNAIL_PATH, 0)
-
-
-def get_correlation_map(
-    img: np.ndarray,
-    template_img: np.ndarray = RBC_THUMBNAIL,
-    downsample_factor: int = 20,
-) -> np.ndarray:
-    """Downsamples the image and returns the 2D cross-correlation map.
-
-    Side note
-    ----------
-    The values generated by the cross-correlation map (using TM_CCOEFF)
-    are sensitive to the dimensions of the image and template. As such, the threshold
-    used is somewhat specific to the particular downsampling factor used.
-
-    TM_CCOEFF_NORMED (which normalizes the values to lie within 0-1) is not used
-    because we want to do a comparison of how well the template matches _between_ images.
-
-    This function has been tested for robustness with a downsample_factor=10
-    and the current template image (image_processing/thumbnail.png).
-
-    Parameters
-    ----------
-    img: np.ndarray
-        Base image
-    template_img: np.ndarray
-        Image to use as a template for the cross-correlation
-    downsample_factor: int=10
-
-    Returns
-    -------
-    np.ndarray:
-        Map of the 2D cross-correlation values
-    """
-
-    img_ds = downsample_image(img, downsample_factor)
-    return cv2.matchTemplate(img_ds, template_img, cv2.TM_CCOEFF)
 
 
 class LowDensity(Exception):
@@ -60,80 +16,47 @@ class NoCellsFound(Exception):
 
 
 class CellFinder:
-    def __init__(
-        self, template_path: str = RBC_THUMBNAIL_PATH, downsample_factor: int = 10
-    ):
-        self.thumbnail = downsample_image(
-            cv2.imread(template_path, 0), downsample_factor
-        )
+    def __init__(self, downsample_factor: int = 10, lookback_window: int = 3):
         self.downsample_factor = downsample_factor
+        self.lookback_window = lookback_window
         self.motor_pos: List[int] = []
-        self.confidences: List[float] = []
-        self.maps: List[np.ndarray] = []
+        self.sds: List[float] = []
 
     def add_image(self, motor_pos: int, img: np.ndarray) -> None:
-        """Check for cells for the given image, store the result + motor position the image was taken at."""
+        """Calculate the standard deviation of the given image (after downsmapling), store the result + motor position the image was taken at."""
+
+        img_ds = downsample_image(img, self.downsample_factor)
+        sd = np.std(img_ds)
 
         self.motor_pos.append(motor_pos)
-        xcorr_map = get_correlation_map(img, self.thumbnail, self.downsample_factor)
-        self.confidences.append(np.max(xcorr_map))
-        self.maps.append(xcorr_map)
+        self.sds.append(sd)
 
     def get_cells_found_position(self) -> Optional[int]:
-        """Check if the cross-correlation value exceeds the threshold for cell detection and there are
-        a sufficient number of points above the cells found threshold.
-
-        Returns
-        -------
-        int
-            Motor position if cells were found
-        Exceptions
-        ----------
-        NoCellsFound
-            Raised if the value of the maximum cross correlation value from the given images
-            does not exceed a threshold.
-        """
-
-        argmax = np.argmax(self.confidences)
-        max_val = self.confidences[argmax]
-
-        if max_val >= CELLS_FOUND_THRESHOLD:
-            return self.motor_pos[np.argmax(self.confidences)]
+        """Find the motor position whose image SD stands out compared to the rest.
 
-        raise NoCellsFound(
-            "None of the images at any of the motor positions had a maximum cross-correlation exceeding the CELLS_FOUND threshold"
-        )
+        Use the median-absolute-deviation and a monotonic decrease over a given lookback
+        period to assess whether the peak a) has objects in it and b) is prominent relative to
+        the surrounding measurements.
+        """
+        sds = np.asarray(self.sds, dtype=float)
+        if sds.size > self.lookback_window:
+            argmax = int(np.argmax(sds))
+            max_val = float(sds[argmax])
 
-    def sufficient_points_above_thresh(self, xcorr_map: np.ndarray) -> bool:
-        """Check if a sufficient number of points are above the cells found threshold
-        in the xcorr map.
+            median = np.median(sds)
+            mad = np.median(np.abs(sds - median))
 
-        Parameters
-        ----------
-        xcorr_map: np.ndarray
+            k = 1.4826  # Scale factor for normally distributed data
+            robust_z = (max_val - median) / (k * mad)
 
-        Returns
-        -------
-        bool
-        """
+            is_monotonically_decreasing = np.all(
+                np.sign(np.diff(sds[-self.lookback_window :])) == -1
+            )
+            if is_monotonically_decreasing and robust_z >= 3.0:
+                return self.motor_pos[argmax]
 
-        points = np.argwhere(xcorr_map >= CELLS_FOUND_THRESHOLD)
-        return len(points) > MIN_POINTS_ABOVE_THRESH
+        raise NoCellsFound("No image stands out as containing cells")
 
     def reset(self) -> None:
         self.motor_pos = []
-        self.confidences = []
-
-    def find_cells_cross_corr(self, img: np.ndarray) -> float:
-        """Returns the max value of the correlation between the RBC thumbnail and the given image (downsampled)
-
-        Returns
-        -------
-        float:
-            Max value of the 2D cross-correlation
-        """
-
-        cross_corr_map = get_correlation_map(
-            img, self.thumbnail, self.downsample_factor
-        )
-        return np.max(cross_corr_map)
+        self.sds = []