Run multiple datasets and evaluate quality of derotation

.pre-commit-config.yaml

@@ -18,14 +18,14 @@ repos:
           - id: requirements-txt-fixer
           - id: trailing-whitespace
     - repo: https://github.com/astral-sh/ruff-pre-commit
-      rev: v0.8.6
+      rev: v0.9.10
       hooks:
         - id: ruff
           args: [ --config=pyproject.toml ]
         - id: ruff-format
           args: [ --config=pyproject.toml ]
     - repo: https://github.com/pre-commit/mirrors-mypy
-      rev: v1.13.0
+      rev: v1.15.0
       hooks:
           - id: mypy
             additional_dependencies:
@@ -39,8 +39,9 @@ repos:
             additional_dependencies: [setuptools-scm, wheel]
     - repo: https://github.com/codespell-project/codespell
       # Configuration for codespell is in pyproject.toml
-      rev: v2.3.0
+      rev: v2.4.1
       hooks:
       - id: codespell
         additional_dependencies:
         - tomli
+        args: ["--ignore-words-list=ptd"]

derotation/analysis/bayesian_optimization.py

@@ -0,0 +1,96 @@
+import logging
+from pathlib import Path
+from typing import Tuple
+
+import matplotlib.pyplot as plt
+import numpy as np
+from bayes_opt import BayesianOptimization
+
+from derotation.analysis.mean_images import calculate_mean_images
+from derotation.analysis.metrics import ptd_of_most_detected_blob
+from derotation.derotate_by_line import derotate_an_image_array_line_by_line
+
+
+class BO_for_derotation:
+    def __init__(
+        self,
+        movie: np.ndarray,
+        rot_deg_line: np.ndarray,
+        rot_deg_frame: np.ndarray,
+        blank_pixels_value: float,
+        center: Tuple[int, int],
+        delta: int,
+        init_points: int = 2,
+        n_iter: int = 10,
+        debug_plots_folder: Path = Path("./debug_plots"),
+    ):
+        self.movie = movie
+        self.rot_deg_line = rot_deg_line
+        self.rot_deg_frame = rot_deg_frame
+        self.blank_pixels_value = blank_pixels_value
+        x, y = center
+        self.pbounds = {
+            "x": (x - delta, x + delta),
+            "y": (y - delta, y + delta),
+        }
+        self.init_points = init_points
+        self.n_iter = n_iter
+        self.debug_plots_folder = debug_plots_folder
+
+        self.subfolder = self.debug_plots_folder / "bo"
+        #  remove previous dir
+        if self.subfolder.exists():
+            for file in self.subfolder.iterdir():
+                file.unlink()
+        else:
+            self.subfolder.mkdir(parents=True, exist_ok=True)
+
+    def optimize(self):
+        def derotate_and_get_metric(
+            x: float,
+            y: float,
+        ):
+            derotated_chunk = derotate_an_image_array_line_by_line(
+                image_stack=self.movie,
+                rot_deg_line=self.rot_deg_line,
+                blank_pixels_value=self.blank_pixels_value,
+                center=(int(x), int(y)),
+            )
+
+            mean_images = calculate_mean_images(
+                derotated_chunk, self.rot_deg_frame, round_decimals=0
+            )
+
+            plt.imshow(mean_images[0], cmap="gray")
+            plt.savefig(self.subfolder / f"mean_image_0_{x:.2f}_{y:.2f}.png")
+            plt.close()
+
+            ptd = ptd_of_most_detected_blob(
+                mean_images,
+                debug_plots_folder=self.subfolder,
+                image_names=[
+                    f"blobs_{x:.2f}_{y:.2f}.png",
+                    f"blob_centers_{x:.2f}_{y:.2f}.png",
+                ],
+            )
+
+            # we are maximizing the metric, so
+            # we need to return the negative of the metric
+            return -ptd
+
+        optimizer = BayesianOptimization(
+            f=derotate_and_get_metric,
+            pbounds=self.pbounds,
+            verbose=2,
+            random_state=1,
+        )
+
+        optimizer.maximize(
+            init_points=self.init_points,
+            n_iter=self.n_iter,
+        )
+
+        for i, res in enumerate(optimizer.res):
+            logging.info(f"Iteration {i}: {res}")
+
+        return optimizer.max

derotation/analysis/blob_detection.py

@@ -96,7 +96,7 @@
         fig, ax = plt.subplots(4, 3, figsize=(10, 10))
         for i, a in tqdm(enumerate(ax.flatten())):
             a.imshow(image_stack[i])
-            a.set_title(f"{i*5} degrees")
+            a.set_title(f"{i * 5} degrees")
             a.axis("off")
 
             for j, blob in enumerate(blobs[i][:4]):

derotation/analysis/fit_ellipse.py

@@ -10,6 +10,7 @@
 
 def fit_ellipse_to_points(
     centers: np.ndarray,
+    pixels_in_row: int = 256,
 ) -> Tuple[int, int, int, int, int]:
     """Fit an ellipse to the points using least squares optimization.
 
@@ -26,14 +27,19 @@
     """
     # Convert centers to numpy array
     centers = np.array(centers)
-    x = centers[:, 0]
-    y = centers[:, 1]
+    valid_points = centers[
+        ~np.isnan(centers).any(axis=1)
+    ]  # Remove rows with NaN
+    if len(valid_points) < 5:
+        raise ValueError("Not enough valid points to fit an ellipse.")
 
-    # Find the extreme points for the initial ellipse estimate
-    topmost = centers[np.argmin(y)]
-    rightmost = centers[np.argmax(x)]
-    bottommost = centers[np.argmax(y)]
-    leftmost = centers[np.argmin(x)]
+    x, y = valid_points[:, 0], valid_points[:, 1]
+
+    # Find extreme points for the initial ellipse estimate
+    topmost = valid_points[np.argmin(y)]
+    rightmost = valid_points[np.argmax(x)]
+    bottommost = valid_points[np.argmax(y)]
+    leftmost = valid_points[np.argmin(x)]
 
     # Initial parameters: (center_x, center_y, semi_major_axis,
     # semi_minor_axis, rotation_angle)
@@ -42,8 +48,12 @@
     )
     semi_major_axis = np.linalg.norm(rightmost - leftmost) / 2
     semi_minor_axis = np.linalg.norm(topmost - bottommost) / 2
-    rotation_angle = 0  # Start with no rotation
 
+    # Ensure axes are not zero
+    if semi_major_axis < 1e-3 or semi_minor_axis < 1e-3:
+        raise ValueError("Points are degenerate; cannot fit an ellipse.")
+
+    rotation_angle = 0  # Start with no rotation
     initial_params = [
         initial_center[0],
         initial_center[1],
@@ -53,6 +63,7 @@
     ]
 
     logging.info("Fitting ellipse to points...")
+    logging.info(f"Initial parameters: {initial_params}")
 
     # Objective function to minimize: sum of squared distances to ellipse
     def ellipse_residuals(params, x, y):
@@ -72,15 +83,26 @@
         ellipse_residuals,
         initial_params,
         args=(x, y),
-        loss="huber",  # minimize the influence of outliers
+        loss="huber",  # Minimize the influence of outliers
+        bounds=(
+            #  center_x, center_y, a, b, theta
+            [0, 0, 1e-3, 1e-3, -np.pi],
+            [
+                pixels_in_row,
+                pixels_in_row,
+                pixels_in_row,
+                pixels_in_row,
+                np.pi,
+            ],
+        ),
     )
 
+    if not result.success:
+        raise RuntimeError("Ellipse fitting did not converge.")
+
     # Extract optimized parameters
     center_x, center_y, a, b, theta = result.x
 
-    #  sometimes the fitted theta is a multiple of of 2pi
-    theta = theta % (2 * np.pi)
-
     return center_x, center_y, a, b, theta