From 039bb923ca3dcc8cd5cd2ed12da9469fdbfac42c Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Fri, 29 Nov 2024 14:19:06 +0000
Subject: [PATCH 01/24] Playing around with dot products

---
 examples/rotating_3D.py | 142 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 142 insertions(+)
 create mode 100644 examples/rotating_3D.py

diff --git a/examples/rotating_3D.py b/examples/rotating_3D.py
new file mode 100644
index 0000000..4e22881
--- /dev/null
+++ b/examples/rotating_3D.py
@@ -0,0 +1,142 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+#  ----------------------------------------------------
+#  first start with two vectors - we are basically projecting a point
+
+#  vector v in the x-y plane laying on the x axis
+v = np.array([[10], [0]])
+
+# vector w in the x-y plane at 45 degrees from v
+w = np.array([[5], [5]])
+
+# project w onto v with the dot product
+w_dot_v = np.dot(w.T, v) / np.dot(v.T, v) * v
+
+# print the result
+print(f"Projection of w ({w}) onto v ({v}): {w_dot_v * v}")
+
+#  plot the vectors
+fig, ax = plt.subplots()
+ax.quiver(0, 0, *v, color="r", angles="xy", scale_units="xy", scale=1)
+ax.quiver(0, 0, *w, color="b", angles="xy", scale_units="xy", scale=1)
+ax.quiver(0, 0, *w_dot_v, color="g", angles="xy", scale_units="xy", scale=1)
+ax.set_xlim(-10, 10)
+ax.set_ylim(-10, 10)
+
+# plt.show()
+
+
+def single_vector_projection(v, w):
+    """
+    Project vector w onto vector v.
+    """
+    return np.dot(w.T, v) / np.dot(v.T, v) * v
+
+
+#  ----------------------------------------------------
+#  now let's project a line, an array of vectors, still in the x-y plane
+
+
+#  make a line of 10 points at 45 degrees from the x axis
+W = np.array([[i, i] for i in range(1, 10)])
+
+#  project the line_45 onto the line_on_x_axis
+W_dot_V = np.array([single_vector_projection(v, w) for w in W])
+
+print(f"Projection of line_45 (\n{W}) onto line_on_x_axis (\n{v}):")
+print(W_dot_V)
+
+#  plot the lines
+fig, ax = plt.subplots()
+for i in range(len(W)):
+    #  random rgb pattern
+    random_color = f"#{np.random.randint(0, 0xFFFFFF):06x}"
+    ax.quiver(
+        0, 0, *W[i], color=random_color, angles="xy", scale_units="xy", scale=1
+    )
+    ax.quiver(
+        0,
+        0,
+        *W_dot_V[i],
+        color=random_color,
+        angles="xy",
+        scale_units="xy",
+        scale=1,
+    )
+
+ax.set_xlim(-10, 10)
+ax.set_ylim(-10, 10)
+
+# plt.show()
+
+#  ----------------------------------------------------
+#  now in 3D space
+
+v_3D = np.array([[10], [0], [0]])
+
+W_3D = np.array([[i, 0, i] for i in range(1, 10)])
+
+W_dot_V_3D = np.array([single_vector_projection(v_3D, w) for w in W_3D])
+
+print(f"Projection of line_45 (\n{W_3D}) onto line_on_x_axis (\n{v_3D}):")
+print(W_dot_V_3D)
+
+#  plot the lines
+fig = plt.figure()
+ax = fig.add_subplot(111, projection="3d")
+for i in range(len(W_3D)):
+    #  random rgb pattern
+    random_color = f"#{np.random.randint(0, 0xFFFFFF):06x}"
+    ax.quiver(0, 0, 0, *W_3D[i], color=random_color)
+    ax.quiver(0, 0, 0, *W_dot_V_3D[i], color=random_color)
+
+ax.set_xlim(-10, 10)
+ax.set_ylim(-10, 10)
+ax.set_zlim(-10, 10)
+
+#  axis labels
+ax.set_xlabel("X")
+ax.set_ylabel("Y")
+ax.set_zlabel("Z")
+
+# plt.show()
+
+#  ----------------------------------------------------
+#  now in 3D space with a plane
+
+#  make a plane where x = 0
+y_plane = np.array([[0, i, j] for i in range(-10, 10) for j in range(-10, 10)])
+
+P_3D = np.array([[i, i, j] for i in range(-10, 10) for j in range(-10, 10)])
+
+#  multiply each dot on the P_3D plane by the corresponding vector in y_plane
+P_dot_V_3D = np.array(
+    [single_vector_projection(y_plane[i], p) for i, p in enumerate(P_3D)]
+)
+
+print(f"Projection of plane (\n{P_3D}) onto y_plane (\n{y_plane}):")
+print(P_dot_V_3D)
+
+#  plot the lines
+fig = plt.figure()
+ax = fig.add_subplot(111, projection="3d")
+for i in range(len(P_3D)):
+    #  random rgb pattern
+    random_color = f"#{np.random.randint(0, 0xFFFFFF):06x}"
+    #  plot them as dots now as they're too many
+    ax.scatter(*P_3D[i], color=random_color)
+    ax.scatter(*P_dot_V_3D[i], color=random_color)
+    #  show the plane also with arrows
+    ax.quiver(0, 0, 0, *y_plane[i], color=random_color, alpha=0.2)
+
+ax.set_xlim(-10, 10)
+ax.set_ylim(-10, 10)
+ax.set_zlim(-10, 10)
+
+#  axis labels
+ax.set_xlabel("X")
+ax.set_ylabel("Y")
+ax.set_zlabel("Z")
+
+plt.show()

From 3f9366eeb494872290e029c822ad4c7ce94f607b Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Fri, 29 Nov 2024 16:54:44 +0000
Subject: [PATCH 02/24] Testing homographies

---
 examples/homographies.py | 47 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 examples/homographies.py

diff --git a/examples/homographies.py b/examples/homographies.py
new file mode 100644
index 0000000..a03f952
--- /dev/null
+++ b/examples/homographies.py
@@ -0,0 +1,47 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import skimage as ski
+from skimage import transform
+
+#  load example  image
+img = ski.data.astronaut()
+
+#  add a stripe of 100 pixels around as padding
+n_pad = 0
+img = np.pad(img, ((n_pad, n_pad), (n_pad, n_pad), (0, 0)), mode="constant")
+
+#  image projected to a plane laying at 45 degrees from the x-y plane
+angle_deg = 45
+
+#  get the target dimensionality of the final projected
+#  image based on the angle
+angle_rad = np.radians(angle_deg)
+
+#  define the homography matrix using the cosine of the angle
+# in the x dimension
+homography_matrix = np.array([[1, 0, 0], [0, np.cos(angle_rad), 0], [0, 0, 1]])
+
+print(f"homography_matrix: {homography_matrix}")
+
+inverse_homography_matrix = np.linalg.inv(homography_matrix)
+print(f"inverse_homography_matrix: {inverse_homography_matrix}")
+
+#  apply the homography
+transformation_1 = transform.warp(img, inverse_homography_matrix)
+
+#  make also inverse homography
+transformation_2 = transform.warp(transformation_1, homography_matrix)
+
+#  plot the images
+fig, ax = plt.subplots(1, 3)
+ax[0].imshow(img)
+ax[0].set_title("Original image")
+ax[1].imshow(transformation_1)
+ax[1].set_title("Transformaion 1")
+ax[2].imshow(transformation_2)
+ax[2].set_title("Transformaion 2")
+
+for a in ax:
+    a.axis("off")
+
+plt.show()

From 4d97b2cf68026caa5e6a6beca77dade29bf155b8 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Fri, 29 Nov 2024 17:56:50 +0000
Subject: [PATCH 03/24] =?UTF-8?q?WIP:=20out=20of=20plane=20rotations=20?=
 =?UTF-8?q?=F0=9F=90=9B=F0=9F=90=9B=F0=9F=90=9B?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .../incremental_derotation_pipeline.py        | 14 +++-
 derotation/simulate/basic_rotator.py          | 53 +++++++++++-
 .../test_integration/rotation_out_of_plane.py | 84 +++++++++++++++++++
 ...er_of_rotation_by_joining_two_pipelines.py | 64 +++++++-------
 4 files changed, 178 insertions(+), 37 deletions(-)
 create mode 100644 tests/test_integration/rotation_out_of_plane.py

diff --git a/derotation/analysis/incremental_derotation_pipeline.py b/derotation/analysis/incremental_derotation_pipeline.py
index 25c1e3e..2723282 100644
--- a/derotation/analysis/incremental_derotation_pipeline.py
+++ b/derotation/analysis/incremental_derotation_pipeline.py
@@ -530,10 +530,16 @@ def get_coords_of_largest_blob(
         blobs = [
             blobs[i][blobs[i][:, 2].argsort()] for i in range(len(image_stack))
         ]
-
-        coord_first_blob_of_every_image = [
-            blobs[i][0][:2].astype(int) for i in range(len(image_stack))
-        ]
+        
+        try:
+            coord_first_blob_of_every_image = [
+                blobs[i][0][:2].astype(int) for i in range(len(image_stack))
+            ]
+        except IndexError:
+            raise ValueError(
+                "No blobs were found. Try changing the parameters of the "
+                + "blob detection algorithm."
+            )
 
         #  invert x, y order
         coord_first_blob_of_every_image = [
diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py
index 26ace59..2044d24 100644
--- a/derotation/simulate/basic_rotator.py
+++ b/derotation/simulate/basic_rotator.py
@@ -11,6 +11,7 @@ def __init__(
         angles: np.ndarray,
         image_stack: np.ndarray,
         center: Optional[Tuple[int, int]] = None,
+        rotation_plane_angle: Optional[float] = None,
     ) -> None:
         """Initializes the Rotator object.
         The Rotator aims to imitate a the scanning pattern of a multi-photon
@@ -60,6 +61,26 @@ def __init__(
             self.center = np.array(image_stack.shape[1:]) / 2
         else:
             self.center = np.array(center)
+        
+        if rotation_plane_angle is None:
+            self.rotation_plane_angle = 0
+        else:
+            self.rotation_plane_angle = rotation_plane_angle
+            self.create_homography_matrices()
+    
+    def create_homography_matrices(self) -> None:
+        #  expansion
+        self.homography_matrix = np.array(
+            [
+                [1, 0, 0],
+                [0, np.cos(np.radians(self.rotation_plane_angle)), 0],
+                [0, 0, 1],
+            ]
+        )
+
+        #  contraction
+        self.inverse_homography_matrix = np.linalg.inv(self.homography_matrix)
+        
 
     def rotate_by_line(self) -> np.ndarray:
         """Simulate the acquisition of a rotated image stack as if for each
@@ -90,6 +111,29 @@ def rotate_by_line(self) -> np.ndarray:
                         rotated_image_stack[i][j] = rotated_frame[j]
 
         return rotated_image_stack
+    
+    def apply_homography(self, image: np.ndarray, direction: str) -> np.ndarray:
+        if direction == "expand":
+            return affine_transform(
+                image,
+                self.homography_matrix,
+                offset=self.center - self.center,
+                output_shape=image.shape,
+                order=0,
+                mode="constant",
+                cval=self.get_blank_pixels_value(),
+            )
+        elif direction == "contract":
+            return affine_transform(
+                image,
+                self.inverse_homography_matrix,
+                offset=self.center - self.center,
+                output_shape=image.shape,
+                order=0,
+                mode="constant",
+                cval=self.get_blank_pixels_value(),
+            )
+
 
     def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
         """Rotate the entire image by a given angle. Uses affine transformation
@@ -110,6 +154,10 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
         np.ndarray
             The rotated image.
         """
+        if self.rotation_plane_angle != 0:
+            image = self.apply_homography(image, "expand")
+
+
         # Compute rotation in radians
         angle_rad = np.deg2rad(angle)
         cos, sin = np.cos(angle_rad), np.sin(angle_rad)
@@ -136,6 +184,9 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
             cval=self.get_blank_pixels_value(),
         )
 
+        if self.rotation_plane_angle != 0:
+            rotated_image = self.apply_homography(rotated_image, "contract")
+
         return rotated_image
 
     def get_blank_pixels_value(self) -> float:
@@ -150,4 +201,4 @@ def get_blank_pixels_value(self) -> float:
         float
             The minimum value of the image stack.
         """
-        return np.min(self.image_stack)
+        return 0 # np.min(self.image_stack)
diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/rotation_out_of_plane.py
new file mode 100644
index 0000000..6d281a4
--- /dev/null
+++ b/tests/test_integration/rotation_out_of_plane.py
@@ -0,0 +1,84 @@
+import numpy as np
+from derotation.simulate.basic_rotator import Rotator
+from matplotlib import pyplot as plt
+
+from test_finding_center_of_rotation_by_joining_two_pipelines import (
+    create_sample_image_with_two_cells,
+    create_image_stack, 
+    create_rotation_angles,
+    MockIncrementalPipeline
+)
+
+plane_angle = 5
+
+
+#  create a sample image with two cells
+cells = create_sample_image_with_two_cells()
+#  chnage all zeros 
+cells[cells == 0] = 80
+
+num_frames = 15
+image_stack = create_image_stack(cells, num_frames=num_frames)
+
+angles, _ = create_rotation_angles(image_stack.shape)
+
+#  rotate the image stack
+rotator = Rotator(angles, image_stack, rotation_plane_angle=plane_angle)
+rotated_image_stack = rotator.rotate_by_line()
+
+
+#  plot the original image, the rotated images 
+fig, ax = plt.subplots(2, num_frames // 2 + 2, figsize=(25, 10))
+
+print(ax.shape)
+ax[0, 0].imshow(image_stack[0], cmap="gray")
+ax[0, 0].set_title("Original image")
+ax[0, 0].axis("off")
+
+for n in range(len(rotated_image_stack) + 1):
+    row = n // 2 + 1
+    col = n % 2
+    ax[row, col].imshow(rotated_image_stack[n], cmap="gray")
+    ax[row, col].set_title(f"Rotated image {n}")
+    ax[row, col].axis("off")
+
+    # add information about the rotation angle
+    angles = rotator.angles[n]
+    angle_range = f"angles: {angles.min():.0f}-{angles.max():.0f}"
+    ax[row, col + 1].text(
+        0.5,
+        0.9,
+        angle_range,
+        horizontalalignment="center",
+        verticalalignment="center",
+        transform=ax[row, col + 1].transAxes,
+        color="white",
+    )
+
+plt.show()
+
+
+#  fit an ellipse to the centers of the largest blobs in each image
+DI = MockIncrementalPipeline(
+    rotated_image_stack, angles
+)
+
+coord_first_blob_of_every_image = DI.get_coords_of_largest_blob(
+        rotated_image_stack
+    )
+
+# Fit an ellipse to the largest blob centers and get its center
+center_x, center_y, a, b, theta = DI.fit_ellipse_to_points(
+    coord_first_blob_of_every_image
+)
+
+DI.plot_ellipse_fit_and_centers(
+    coord_first_blob_of_every_image,
+    center_x,
+    center_y,
+    a,
+    b,
+    theta,
+    rotated_image_stack,
+    num_frames,
+)
\ No newline at end of file
diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
index a5703f9..ff03992 100644
--- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
+++ b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
@@ -180,6 +180,37 @@ def create_rotation_angles(
 #  -----------------------------------------------------
 
 
+# Mock class to use the IncrementalPipeline
+class MockIncrementalPipeline(IncrementalPipeline):
+    def __init__(self, rotated_stack_incremental, incremental_angles):
+        # Overwrite the constructor and provide the mock data
+        self.image_stack = rotated_stack_incremental
+        self.rot_deg_frame = incremental_angles[
+            :: rotated_stack_incremental.shape[1]
+        ]
+        self.num_frames = rotated_stack_incremental.shape[0]
+
+        if __name__ == "__main__":
+            self.debugging_plots = True
+            self.debug_plots_folder = Path("debug/")
+        else:
+            self.debugging_plots = False
+
+    def calculate_mean_images(self, image_stack: np.ndarray) -> list:
+        #  Overwrite original method as it is too bound
+        #  to signal coming from a real motor
+        angles_subset = copy.deepcopy(self.rot_deg_frame)
+        rounded_angles = np.round(angles_subset)
+
+        mean_images = []
+        for i in np.arange(10, 360, 10):
+            images = image_stack[rounded_angles == i]
+            mean_image = np.mean(images, axis=0)
+
+            mean_images.append(mean_image)
+
+        return mean_images
+
 def get_center_of_rotation(
     rotated_stack_incremental: np.ndarray, incremental_angles: np.ndarray
 ) -> Tuple[int, int]:
@@ -206,39 +237,8 @@ def get_center_of_rotation(
         The center of rotation
     """
 
-    # Mock class to use the IncrementalPipeline
-    class MockIncrementalPipeline(IncrementalPipeline):
-        def __init__(self):
-            # Overwrite the constructor and provide the mock data
-            self.image_stack = rotated_stack_incremental
-            self.rot_deg_frame = incremental_angles[
-                :: rotated_stack_incremental.shape[1]
-            ]
-            self.num_frames = rotated_stack_incremental.shape[0]
-
-            if __name__ == "__main__":
-                self.debugging_plots = True
-                self.debug_plots_folder = Path("debug/")
-            else:
-                self.debugging_plots = False
-
-        def calculate_mean_images(self, image_stack: np.ndarray) -> list:
-            #  Overwrite original method as it is too bound
-            #  to signal coming from a real motor
-            angles_subset = copy.deepcopy(self.rot_deg_frame)
-            rounded_angles = np.round(angles_subset)
-
-            mean_images = []
-            for i in np.arange(10, 360, 10):
-                images = image_stack[rounded_angles == i]
-                mean_image = np.mean(images, axis=0)
-
-                mean_images.append(mean_image)
-
-            return mean_images
-
     # Use the mock class to find the center of rotation
-    pipeline = MockIncrementalPipeline()
+    pipeline = MockIncrementalPipeline(rotated_stack_incremental, incremental_angles)
     center_of_rotation = pipeline.find_center_of_rotation()
 
     return center_of_rotation

From c8a36c0616802c6382ec2b7e23d0f38f894a862f Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 2 Dec 2024 16:14:17 +0000
Subject: [PATCH 04/24] Add working simulation of rotation out of plane

---
 derotation/simulate/basic_rotator.py          | 103 ++++++++++++++----
 .../test_integration/rotation_out_of_plane.py |  89 +++++++--------
 ...er_of_rotation_by_joining_two_pipelines.py |   8 +-
 3 files changed, 127 insertions(+), 73 deletions(-)

diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py
index 2044d24..9def7e5 100644
--- a/derotation/simulate/basic_rotator.py
+++ b/derotation/simulate/basic_rotator.py
@@ -1,7 +1,7 @@
-import copy
 from typing import Optional, Tuple
 
 import numpy as np
+import tqdm
 from scipy.ndimage import affine_transform
 
 
@@ -50,10 +50,6 @@ def __init__(
             + f"({image_stack.shape[0] * image_stack.shape[1]})"
         )
 
-        #  reshape the angles to the shape of the image stack to ease indexing
-        self.angles = angles.reshape(
-            image_stack.shape[0], image_stack.shape[1]
-        )
         self.image_stack = image_stack
         self.num_lines_per_frame = image_stack.shape[1]
 
@@ -61,15 +57,28 @@ def __init__(
             self.center = np.array(image_stack.shape[1:]) / 2
         else:
             self.center = np.array(center)
-        
+
         if rotation_plane_angle is None:
-            self.rotation_plane_angle = 0
+            self.rotation_plane_angle: float = 0
+            #  reshape the angles to the shape of the image
+            #  stack to ease indexing
+            self.angles = angles.reshape(
+                image_stack.shape[0], image_stack.shape[1]
+            )
         else:
             self.rotation_plane_angle = rotation_plane_angle
             self.create_homography_matrices()
-    
+            print(f"Pixel shift: {self.ps}")
+            print(f"New image size: {self.image_size}")
+            #  reshape angles depending on new image size
+            #  angles might be cropped
+            self.angles = angles[
+                : image_stack.shape[0] * self.image_size
+            ].reshape(image_stack.shape[0], self.image_size)
+            print("New angles shape:", self.angles.shape)
+
     def create_homography_matrices(self) -> None:
-        #  expansion
+        #  from the scanning plane to the rotation plane
         self.homography_matrix = np.array(
             [
                 [1, 0, 0],
@@ -78,9 +87,35 @@ def create_homography_matrices(self) -> None:
             ]
         )
 
-        #  contraction
+        #  from the rotation plane to the scanning plane
         self.inverse_homography_matrix = np.linalg.inv(self.homography_matrix)
-        
+
+        #  store pixels shift based on inverse homography
+        line_length = self.image_stack.shape[2]
+        self.ps = (
+            line_length
+            - np.round(
+                np.abs(
+                    line_length * np.cos(np.radians(self.rotation_plane_angle))
+                )
+            ).astype(int)
+            + 1
+        )
+
+        #  round to the nearest even number
+        if self.ps % 2 != 0:
+            self.ps -= 1
+
+        #  final image size depends on the pixel shift
+        self.image_size = line_length - self.ps
+
+    def crop_image(self, image: np.ndarray) -> np.ndarray:
+        return image[
+            # centered in the rows
+            self.ps // 2 : -self.ps // 2,
+            # take the left side of the image
+            : self.image_size,
+        ]
 
     def rotate_by_line(self) -> np.ndarray:
         """Simulate the acquisition of a rotated image stack as if for each
@@ -95,9 +130,14 @@ def rotate_by_line(self) -> np.ndarray:
             The rotated image stack of the same shape as the input image stack,
             i.e. (num_frames, num_lines_per_frame, num_pixels_per_line).
         """
-        rotated_image_stack = copy.deepcopy(self.image_stack)
+        rotated_image_stack = np.empty(
+            (self.image_stack.shape[0], self.image_size, self.image_size),
+            dtype=np.float64,
+        )
 
-        for i, image in enumerate(self.image_stack):
+        for i, image in tqdm.tqdm(
+            enumerate(self.image_stack), total=self.image_stack.shape[0]
+        ):
             is_this_frame_rotating = not np.all(
                 # don't bother if rotation is less than 0.01 degrees
                 np.isclose(self.angles[i], 0, atol=1e-2)
@@ -105,15 +145,21 @@ def rotate_by_line(self) -> np.ndarray:
             if is_this_frame_rotating:
                 for j, angle in enumerate(self.angles[i]):
                     if angle == 0:
-                        continue
+                        rotated_image_stack[i][j] = self.crop_image(image)[j]
                     else:
                         rotated_frame = self.rotate(image, angle)
                         rotated_image_stack[i][j] = rotated_frame[j]
+                        # plt.imshow(rotated_image_stack[i])
+                        # plt.savefig(f"debug/projections/image_{i}_{j}.png")
+                        # plt.close()
 
         return rotated_image_stack
-    
-    def apply_homography(self, image: np.ndarray, direction: str) -> np.ndarray:
-        if direction == "expand":
+
+    def apply_homography(
+        self, image: np.ndarray, direction: str
+    ) -> np.ndarray:
+        if direction == "scanning_to_rotation_plane":
+            # backward transformation
             return affine_transform(
                 image,
                 self.homography_matrix,
@@ -123,7 +169,8 @@ def apply_homography(self, image: np.ndarray, direction: str) -> np.ndarray:
                 mode="constant",
                 cval=self.get_blank_pixels_value(),
             )
-        elif direction == "contract":
+        elif direction == "rotation_to_scanning_plane":
+            # forward transformation
             return affine_transform(
                 image,
                 self.inverse_homography_matrix,
@@ -134,7 +181,6 @@ def apply_homography(self, image: np.ndarray, direction: str) -> np.ndarray:
                 cval=self.get_blank_pixels_value(),
             )
 
-
     def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
         """Rotate the entire image by a given angle. Uses affine transformation
         with no interpolation.
@@ -154,9 +200,6 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
         np.ndarray
             The rotated image.
         """
-        if self.rotation_plane_angle != 0:
-            image = self.apply_homography(image, "expand")
-
 
         # Compute rotation in radians
         angle_rad = np.deg2rad(angle)
@@ -185,7 +228,19 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
         )
 
         if self.rotation_plane_angle != 0:
-            rotated_image = self.apply_homography(rotated_image, "contract")
+            rotated_image = self.apply_homography(
+                rotated_image, "rotation_to_scanning_plane"
+            )
+
+        # plt.imshow(rotated_image, vmin=0, vmax=255)
+        # plt.savefig(f"debug/projections/rotated_image_projected{angle:.3f}.png")
+        # plt.close()
+
+        rotated_image = self.crop_image(rotated_image)
+
+        # plt.imshow(rotated_image, vmin=0, vmax=255)
+        # plt.savefig(f"debug/projections/rotated_image_projected_cropped{angle:.3f}.png")
+        # plt.close()
 
         return rotated_image
 
@@ -201,4 +256,4 @@ def get_blank_pixels_value(self) -> float:
         float
             The minimum value of the image stack.
         """
-        return 0 # np.min(self.image_stack)
+        return 0  # np.min(self.image_stack)
diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/rotation_out_of_plane.py
index 6d281a4..d69c610 100644
--- a/tests/test_integration/rotation_out_of_plane.py
+++ b/tests/test_integration/rotation_out_of_plane.py
@@ -1,84 +1,79 @@
 import numpy as np
-from derotation.simulate.basic_rotator import Rotator
 from matplotlib import pyplot as plt
-
 from test_finding_center_of_rotation_by_joining_two_pipelines import (
-    create_sample_image_with_two_cells,
-    create_image_stack, 
+    create_image_stack,
     create_rotation_angles,
-    MockIncrementalPipeline
+    create_sample_image_with_two_cells,
 )
 
-plane_angle = 5
+from derotation.simulate.basic_rotator import Rotator
 
+plane_angle = 25
+num_frames = 50
+pad = 20
 
 #  create a sample image with two cells
-cells = create_sample_image_with_two_cells()
-#  chnage all zeros 
+cells = create_sample_image_with_two_cells(lines_per_frame=100)
+cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
 cells[cells == 0] = 80
 
-num_frames = 15
+
 image_stack = create_image_stack(cells, num_frames=num_frames)
+print("Image stack shape:", image_stack.shape)
 
-angles, _ = create_rotation_angles(image_stack.shape)
+_, angles = create_rotation_angles(image_stack.shape)
+print("Angles shape:", angles.shape)
 
 #  rotate the image stack
 rotator = Rotator(angles, image_stack, rotation_plane_angle=plane_angle)
 rotated_image_stack = rotator.rotate_by_line()
 
+print("Rotated image stack shape:", rotated_image_stack.shape)
 
-#  plot the original image, the rotated images 
-fig, ax = plt.subplots(2, num_frames // 2 + 2, figsize=(25, 10))
+#  plot the original image, the rotated images
+row_n = 5
+fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 5))
 
 print(ax.shape)
-ax[0, 0].imshow(image_stack[0], cmap="gray")
+ax[0, 0].imshow(image_stack[0], cmap="gray", vmin=0, vmax=255)
 ax[0, 0].set_title("Original image")
 ax[0, 0].axis("off")
 
-for n in range(len(rotated_image_stack) + 1):
-    row = n // 2 + 1
-    col = n % 2
-    ax[row, col].imshow(rotated_image_stack[n], cmap="gray")
-    ax[row, col].set_title(f"Rotated image {n}")
-    ax[row, col].axis("off")
+for n in range(1, len(rotated_image_stack) + 1):
+    row = n % row_n
+    col = n // row_n
+    ax[row, col].imshow(
+        rotated_image_stack[n - 1], cmap="gray", vmin=0, vmax=255
+    )
+    ax[row, col].set_title(n)
 
     # add information about the rotation angle
-    angles = rotator.angles[n]
-    angle_range = f"angles: {angles.min():.0f}-{angles.max():.0f}"
-    ax[row, col + 1].text(
+    angles = rotator.angles[n - 1]
+    angle_range = f"{angles.min():.0f}-{angles.max():.0f}"
+    ax[row, col].text(
         0.5,
         0.9,
         angle_range,
         horizontalalignment="center",
         verticalalignment="center",
-        transform=ax[row, col + 1].transAxes,
+        transform=ax[row, col].transAxes,
         color="white",
     )
 
-plt.show()
-
-
-#  fit an ellipse to the centers of the largest blobs in each image
-DI = MockIncrementalPipeline(
-    rotated_image_stack, angles
+#  last one, max projection
+ax[row_n - 1, num_frames // row_n].imshow(
+    rotated_image_stack.max(axis=0), cmap="gray"
 )
-
-coord_first_blob_of_every_image = DI.get_coords_of_largest_blob(
-        rotated_image_stack
-    )
-
-# Fit an ellipse to the largest blob centers and get its center
-center_x, center_y, a, b, theta = DI.fit_ellipse_to_points(
-    coord_first_blob_of_every_image
+#  plot also an x in the center of the image in red
+ax[row_n - 1, num_frames // row_n].plot(
+    rotated_image_stack.shape[2] / 2,
+    rotated_image_stack.shape[1] / 2,
+    "rx",
+    markersize=10,
 )
+ax[row_n - 1, num_frames // row_n].set_title("Max projection")
+
+for a in ax.ravel():
+    a.axis("off")
 
-DI.plot_ellipse_fit_and_centers(
-    coord_first_blob_of_every_image,
-    center_x,
-    center_y,
-    a,
-    b,
-    theta,
-    rotated_image_stack,
-    num_frames,
-)
\ No newline at end of file
+plt.show()
diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
index ff03992..c3c2d6b 100644
--- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
+++ b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
@@ -63,6 +63,7 @@
 def create_sample_image_with_two_cells(
     center_of_bright_cell: Tuple[int, int] = (50, 10),
     center_of_dimmer_cell: Tuple[int, int] = (60, 60),
+    lines_per_frame: int = 100,
 ) -> np.ndarray:
     """Create a 2D image with two circles, one bright and one dim
     by default in the top center and bottom right, respectively.
@@ -84,7 +85,7 @@ def create_sample_image_with_two_cells(
     """
 
     # Initialize a black image of size 100x100
-    image = np.zeros((100, 100), dtype=np.uint8)
+    image = np.zeros((lines_per_frame, lines_per_frame), dtype=np.uint8)
 
     # Define the circle's parameters
     radius = 5  # radius of the circle
@@ -211,6 +212,7 @@ def calculate_mean_images(self, image_stack: np.ndarray) -> list:
 
         return mean_images
 
+
 def get_center_of_rotation(
     rotated_stack_incremental: np.ndarray, incremental_angles: np.ndarray
 ) -> Tuple[int, int]:
@@ -238,7 +240,9 @@ def get_center_of_rotation(
     """
 
     # Use the mock class to find the center of rotation
-    pipeline = MockIncrementalPipeline(rotated_stack_incremental, incremental_angles)
+    pipeline = MockIncrementalPipeline(
+        rotated_stack_incremental, incremental_angles
+    )
     center_of_rotation = pipeline.find_center_of_rotation()
 
     return center_of_rotation

From 5f6dd711f9525473c45dab8e607a50b9b8f38ed4 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 2 Dec 2024 16:50:16 +0000
Subject: [PATCH 05/24] Fix some bugs found via pytest

---
 derotation/simulate/basic_rotator.py | 39 ++++++++++++++--------------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py
index 9def7e5..08ee0d3 100644
--- a/derotation/simulate/basic_rotator.py
+++ b/derotation/simulate/basic_rotator.py
@@ -60,22 +60,21 @@ def __init__(
 
         if rotation_plane_angle is None:
             self.rotation_plane_angle: float = 0
-            #  reshape the angles to the shape of the image
-            #  stack to ease indexing
-            self.angles = angles.reshape(
-                image_stack.shape[0], image_stack.shape[1]
-            )
+            self.ps: int = 0
+            self.image_size = image_stack.shape[1]
         else:
             self.rotation_plane_angle = rotation_plane_angle
+
             self.create_homography_matrices()
             print(f"Pixel shift: {self.ps}")
             print(f"New image size: {self.image_size}")
-            #  reshape angles depending on new image size
-            #  angles might be cropped
-            self.angles = angles[
-                : image_stack.shape[0] * self.image_size
-            ].reshape(image_stack.shape[0], self.image_size)
-            print("New angles shape:", self.angles.shape)
+
+        #  reshape the angles to the shape of the image
+        #  stack to ease indexing
+        self.angles = angles[: image_stack.shape[0] * self.image_size].reshape(
+            image_stack.shape[0], self.image_size
+        )
+        print("New angles shape:", self.angles.shape)
 
     def create_homography_matrices(self) -> None:
         #  from the scanning plane to the rotation plane
@@ -110,12 +109,15 @@ def create_homography_matrices(self) -> None:
         self.image_size = line_length - self.ps
 
     def crop_image(self, image: np.ndarray) -> np.ndarray:
-        return image[
-            # centered in the rows
-            self.ps // 2 : -self.ps // 2,
-            # take the left side of the image
-            : self.image_size,
-        ]
+        if self.ps == 0:
+            return image
+        else:
+            return image[
+                # centered in the rows
+                self.ps // 2 : -self.ps // 2,
+                # take the left side of the image
+                : self.image_size,
+            ]
 
     def rotate_by_line(self) -> np.ndarray:
         """Simulate the acquisition of a rotated image stack as if for each
@@ -231,13 +233,12 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
             rotated_image = self.apply_homography(
                 rotated_image, "rotation_to_scanning_plane"
             )
+            rotated_image = self.crop_image(rotated_image)
 
         # plt.imshow(rotated_image, vmin=0, vmax=255)
         # plt.savefig(f"debug/projections/rotated_image_projected{angle:.3f}.png")
         # plt.close()
 
-        rotated_image = self.crop_image(rotated_image)
-
         # plt.imshow(rotated_image, vmin=0, vmax=255)
         # plt.savefig(f"debug/projections/rotated_image_projected_cropped{angle:.3f}.png")
         # plt.close()

From f159f612cf286233e33596dd50865025c90a56ee Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 2 Dec 2024 16:50:41 +0000
Subject: [PATCH 06/24] Remove exploratory scripts

---
 examples/homographies.py |  47 -------------
 examples/rotating_3D.py  | 142 ---------------------------------------
 2 files changed, 189 deletions(-)
 delete mode 100644 examples/homographies.py
 delete mode 100644 examples/rotating_3D.py

diff --git a/examples/homographies.py b/examples/homographies.py
deleted file mode 100644
index a03f952..0000000
--- a/examples/homographies.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-import skimage as ski
-from skimage import transform
-
-#  load example  image
-img = ski.data.astronaut()
-
-#  add a stripe of 100 pixels around as padding
-n_pad = 0
-img = np.pad(img, ((n_pad, n_pad), (n_pad, n_pad), (0, 0)), mode="constant")
-
-#  image projected to a plane laying at 45 degrees from the x-y plane
-angle_deg = 45
-
-#  get the target dimensionality of the final projected
-#  image based on the angle
-angle_rad = np.radians(angle_deg)
-
-#  define the homography matrix using the cosine of the angle
-# in the x dimension
-homography_matrix = np.array([[1, 0, 0], [0, np.cos(angle_rad), 0], [0, 0, 1]])
-
-print(f"homography_matrix: {homography_matrix}")
-
-inverse_homography_matrix = np.linalg.inv(homography_matrix)
-print(f"inverse_homography_matrix: {inverse_homography_matrix}")
-
-#  apply the homography
-transformation_1 = transform.warp(img, inverse_homography_matrix)
-
-#  make also inverse homography
-transformation_2 = transform.warp(transformation_1, homography_matrix)
-
-#  plot the images
-fig, ax = plt.subplots(1, 3)
-ax[0].imshow(img)
-ax[0].set_title("Original image")
-ax[1].imshow(transformation_1)
-ax[1].set_title("Transformaion 1")
-ax[2].imshow(transformation_2)
-ax[2].set_title("Transformaion 2")
-
-for a in ax:
-    a.axis("off")
-
-plt.show()
diff --git a/examples/rotating_3D.py b/examples/rotating_3D.py
deleted file mode 100644
index 4e22881..0000000
--- a/examples/rotating_3D.py
+++ /dev/null
@@ -1,142 +0,0 @@
-import matplotlib.pyplot as plt
-import numpy as np
-
-#  ----------------------------------------------------
-#  first start with two vectors - we are basically projecting a point
-
-#  vector v in the x-y plane laying on the x axis
-v = np.array([[10], [0]])
-
-# vector w in the x-y plane at 45 degrees from v
-w = np.array([[5], [5]])
-
-# project w onto v with the dot product
-w_dot_v = np.dot(w.T, v) / np.dot(v.T, v) * v
-
-# print the result
-print(f"Projection of w ({w}) onto v ({v}): {w_dot_v * v}")
-
-#  plot the vectors
-fig, ax = plt.subplots()
-ax.quiver(0, 0, *v, color="r", angles="xy", scale_units="xy", scale=1)
-ax.quiver(0, 0, *w, color="b", angles="xy", scale_units="xy", scale=1)
-ax.quiver(0, 0, *w_dot_v, color="g", angles="xy", scale_units="xy", scale=1)
-ax.set_xlim(-10, 10)
-ax.set_ylim(-10, 10)
-
-# plt.show()
-
-
-def single_vector_projection(v, w):
-    """
-    Project vector w onto vector v.
-    """
-    return np.dot(w.T, v) / np.dot(v.T, v) * v
-
-
-#  ----------------------------------------------------
-#  now let's project a line, an array of vectors, still in the x-y plane
-
-
-#  make a line of 10 points at 45 degrees from the x axis
-W = np.array([[i, i] for i in range(1, 10)])
-
-#  project the line_45 onto the line_on_x_axis
-W_dot_V = np.array([single_vector_projection(v, w) for w in W])
-
-print(f"Projection of line_45 (\n{W}) onto line_on_x_axis (\n{v}):")
-print(W_dot_V)
-
-#  plot the lines
-fig, ax = plt.subplots()
-for i in range(len(W)):
-    #  random rgb pattern
-    random_color = f"#{np.random.randint(0, 0xFFFFFF):06x}"
-    ax.quiver(
-        0, 0, *W[i], color=random_color, angles="xy", scale_units="xy", scale=1
-    )
-    ax.quiver(
-        0,
-        0,
-        *W_dot_V[i],
-        color=random_color,
-        angles="xy",
-        scale_units="xy",
-        scale=1,
-    )
-
-ax.set_xlim(-10, 10)
-ax.set_ylim(-10, 10)
-
-# plt.show()
-
-#  ----------------------------------------------------
-#  now in 3D space
-
-v_3D = np.array([[10], [0], [0]])
-
-W_3D = np.array([[i, 0, i] for i in range(1, 10)])
-
-W_dot_V_3D = np.array([single_vector_projection(v_3D, w) for w in W_3D])
-
-print(f"Projection of line_45 (\n{W_3D}) onto line_on_x_axis (\n{v_3D}):")
-print(W_dot_V_3D)
-
-#  plot the lines
-fig = plt.figure()
-ax = fig.add_subplot(111, projection="3d")
-for i in range(len(W_3D)):
-    #  random rgb pattern
-    random_color = f"#{np.random.randint(0, 0xFFFFFF):06x}"
-    ax.quiver(0, 0, 0, *W_3D[i], color=random_color)
-    ax.quiver(0, 0, 0, *W_dot_V_3D[i], color=random_color)
-
-ax.set_xlim(-10, 10)
-ax.set_ylim(-10, 10)
-ax.set_zlim(-10, 10)
-
-#  axis labels
-ax.set_xlabel("X")
-ax.set_ylabel("Y")
-ax.set_zlabel("Z")
-
-# plt.show()
-
-#  ----------------------------------------------------
-#  now in 3D space with a plane
-
-#  make a plane where x = 0
-y_plane = np.array([[0, i, j] for i in range(-10, 10) for j in range(-10, 10)])
-
-P_3D = np.array([[i, i, j] for i in range(-10, 10) for j in range(-10, 10)])
-
-#  multiply each dot on the P_3D plane by the corresponding vector in y_plane
-P_dot_V_3D = np.array(
-    [single_vector_projection(y_plane[i], p) for i, p in enumerate(P_3D)]
-)
-
-print(f"Projection of plane (\n{P_3D}) onto y_plane (\n{y_plane}):")
-print(P_dot_V_3D)
-
-#  plot the lines
-fig = plt.figure()
-ax = fig.add_subplot(111, projection="3d")
-for i in range(len(P_3D)):
-    #  random rgb pattern
-    random_color = f"#{np.random.randint(0, 0xFFFFFF):06x}"
-    #  plot them as dots now as they're too many
-    ax.scatter(*P_3D[i], color=random_color)
-    ax.scatter(*P_dot_V_3D[i], color=random_color)
-    #  show the plane also with arrows
-    ax.quiver(0, 0, 0, *y_plane[i], color=random_color, alpha=0.2)
-
-ax.set_xlim(-10, 10)
-ax.set_ylim(-10, 10)
-ax.set_zlim(-10, 10)
-
-#  axis labels
-ax.set_xlabel("X")
-ax.set_ylabel("Y")
-ax.set_zlabel("Z")
-
-plt.show()

From b68494c077def9ba56b1d5433b49db1115d5c563 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 12:29:09 +0000
Subject: [PATCH 07/24] Remove comments from basic rotator and fix bug

---
 derotation/simulate/basic_rotator.py | 15 +++------------
 1 file changed, 3 insertions(+), 12 deletions(-)

diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py
index 08ee0d3..596c1ff 100644
--- a/derotation/simulate/basic_rotator.py
+++ b/derotation/simulate/basic_rotator.py
@@ -151,9 +151,8 @@ def rotate_by_line(self) -> np.ndarray:
                     else:
                         rotated_frame = self.rotate(image, angle)
                         rotated_image_stack[i][j] = rotated_frame[j]
-                        # plt.imshow(rotated_image_stack[i])
-                        # plt.savefig(f"debug/projections/image_{i}_{j}.png")
-                        # plt.close()
+            else:
+                rotated_image_stack[i] = self.crop_image(image)
 
         return rotated_image_stack
 
@@ -235,14 +234,6 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
             )
             rotated_image = self.crop_image(rotated_image)
 
-        # plt.imshow(rotated_image, vmin=0, vmax=255)
-        # plt.savefig(f"debug/projections/rotated_image_projected{angle:.3f}.png")
-        # plt.close()
-
-        # plt.imshow(rotated_image, vmin=0, vmax=255)
-        # plt.savefig(f"debug/projections/rotated_image_projected_cropped{angle:.3f}.png")
-        # plt.close()
-
         return rotated_image
 
     def get_blank_pixels_value(self) -> float:
@@ -257,4 +248,4 @@ def get_blank_pixels_value(self) -> float:
         float
             The minimum value of the image stack.
         """
-        return 0  # np.min(self.image_stack)
+        return np.min(self.image_stack)

From 0601588c6d199966afaac6cccfd58b0a1c9aa9c8 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 12:30:16 +0000
Subject: [PATCH 08/24] Give upper and lower bounds to blob log

---
 ...est_finding_center_of_rotation_by_joining_two_pipelines.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
index c3c2d6b..7bf2603 100644
--- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
+++ b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
@@ -475,7 +475,9 @@ def test_blob_detection_on_derotated_stack(
     # If not, the derotation was not successful and the test fails
 
     # Detect the blobs in the derotated stack
-    blobs = [blob_log(img) for img in derotated_sinusoidal]
+    blobs = [
+        blob_log(img, min_sigma=3, max_sigma=5) for img in derotated_sinusoidal
+    ]
 
     # Get the center of the blobs
     # for every frame, place first the blob with the smallest x value

From 31a63a1c3121e56f1f6be56fadb28ca12d844d15 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 12:32:58 +0000
Subject: [PATCH 09/24] Bring mock class back in the method where it was

---
 ...er_of_rotation_by_joining_two_pipelines.py | 67 +++++++++----------
 1 file changed, 32 insertions(+), 35 deletions(-)

diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
index 7bf2603..0aa3f69 100644
--- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
+++ b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
@@ -181,38 +181,6 @@ def create_rotation_angles(
 #  -----------------------------------------------------
 
 
-# Mock class to use the IncrementalPipeline
-class MockIncrementalPipeline(IncrementalPipeline):
-    def __init__(self, rotated_stack_incremental, incremental_angles):
-        # Overwrite the constructor and provide the mock data
-        self.image_stack = rotated_stack_incremental
-        self.rot_deg_frame = incremental_angles[
-            :: rotated_stack_incremental.shape[1]
-        ]
-        self.num_frames = rotated_stack_incremental.shape[0]
-
-        if __name__ == "__main__":
-            self.debugging_plots = True
-            self.debug_plots_folder = Path("debug/")
-        else:
-            self.debugging_plots = False
-
-    def calculate_mean_images(self, image_stack: np.ndarray) -> list:
-        #  Overwrite original method as it is too bound
-        #  to signal coming from a real motor
-        angles_subset = copy.deepcopy(self.rot_deg_frame)
-        rounded_angles = np.round(angles_subset)
-
-        mean_images = []
-        for i in np.arange(10, 360, 10):
-            images = image_stack[rounded_angles == i]
-            mean_image = np.mean(images, axis=0)
-
-            mean_images.append(mean_image)
-
-        return mean_images
-
-
 def get_center_of_rotation(
     rotated_stack_incremental: np.ndarray, incremental_angles: np.ndarray
 ) -> Tuple[int, int]:
@@ -239,10 +207,39 @@ def get_center_of_rotation(
         The center of rotation
     """
 
+    # Mock class to use the IncrementalPipeline
+    class MockIncrementalPipeline(IncrementalPipeline):
+        def __init__(self):
+            # Overwrite the constructor and provide the mock data
+            self.image_stack = rotated_stack_incremental
+            self.rot_deg_frame = incremental_angles[
+                :: rotated_stack_incremental.shape[1]
+            ]
+            self.num_frames = rotated_stack_incremental.shape[0]
+
+            if __name__ == "__main__":
+                self.debugging_plots = True
+                self.debug_plots_folder = Path("debug/")
+            else:
+                self.debugging_plots = False
+
+        def calculate_mean_images(self, image_stack: np.ndarray) -> list:
+            #  Overwrite original method as it is too bound
+            #  to signal coming from a real motor
+            angles_subset = copy.deepcopy(self.rot_deg_frame)
+            rounded_angles = np.round(angles_subset)
+
+            mean_images = []
+            for i in np.arange(10, 360, 10):
+                images = image_stack[rounded_angles == i]
+                mean_image = np.mean(images, axis=0)
+
+                mean_images.append(mean_image)
+
+            return mean_images
+
     # Use the mock class to find the center of rotation
-    pipeline = MockIncrementalPipeline(
-        rotated_stack_incremental, incremental_angles
-    )
+    pipeline = MockIncrementalPipeline()
     center_of_rotation = pipeline.find_center_of_rotation()
 
     return center_of_rotation

From 7d1a07293556708646d258f1ddcb283bbc700ad9 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 14:19:55 +0000
Subject: [PATCH 10/24] Small refactoring

---
 derotation/simulate/basic_rotator.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py
index 596c1ff..fe324e0 100644
--- a/derotation/simulate/basic_rotator.py
+++ b/derotation/simulate/basic_rotator.py
@@ -12,6 +12,7 @@ def __init__(
         image_stack: np.ndarray,
         center: Optional[Tuple[int, int]] = None,
         rotation_plane_angle: Optional[float] = None,
+        blank_pixel_val: Optional[float] = None,
     ) -> None:
         """Initializes the Rotator object.
         The Rotator aims to imitate a the scanning pattern of a multi-photon
@@ -76,6 +77,9 @@ def __init__(
         )
         print("New angles shape:", self.angles.shape)
 
+        if blank_pixel_val is None:
+            self.blank_pixel_val = self.get_blank_pixels_value()
+
     def create_homography_matrices(self) -> None:
         #  from the scanning plane to the rotation plane
         self.homography_matrix = np.array(
@@ -149,7 +153,7 @@ def rotate_by_line(self) -> np.ndarray:
                     if angle == 0:
                         rotated_image_stack[i][j] = self.crop_image(image)[j]
                     else:
-                        rotated_frame = self.rotate(image, angle)
+                        rotated_frame = self.rotate_sample(image, angle)
                         rotated_image_stack[i][j] = rotated_frame[j]
             else:
                 rotated_image_stack[i] = self.crop_image(image)
@@ -182,7 +186,7 @@ def apply_homography(
                 cval=self.get_blank_pixels_value(),
             )
 
-    def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
+    def rotate_sample(self, image: np.ndarray, angle: float) -> np.ndarray:
         """Rotate the entire image by a given angle. Uses affine transformation
         with no interpolation.
 
@@ -225,7 +229,7 @@ def rotate(self, image: np.ndarray, angle: float) -> np.ndarray:
             output_shape=image.shape,  # Keep original shape
             order=0,
             mode="constant",  # NO interpolation
-            cval=self.get_blank_pixels_value(),
+            cval=self.blank_pixel_val,
         )
 
         if self.rotation_plane_angle != 0:

From 6a55da2c13aa4649fd12878cda05cb95004d988b Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 14:21:35 +0000
Subject: [PATCH 11/24] Rename module

---
 .../simulate/{basic_rotator.py => line_scanning_microscope.py}  | 0
 examples/rotate_and_derotate_a_square.py                        | 2 +-
 .../test_finding_center_of_rotation_by_joining_two_pipelines.py | 2 +-
 tests/test_regression/recreate_target/shared.py                 | 2 +-
 .../test_basic_rotator_with_different_center_of_rotation.py     | 2 +-
 tests/test_unit/test_rotator.py                                 | 2 +-
 6 files changed, 5 insertions(+), 5 deletions(-)
 rename derotation/simulate/{basic_rotator.py => line_scanning_microscope.py} (100%)

diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/line_scanning_microscope.py
similarity index 100%
rename from derotation/simulate/basic_rotator.py
rename to derotation/simulate/line_scanning_microscope.py
diff --git a/examples/rotate_and_derotate_a_square.py b/examples/rotate_and_derotate_a_square.py
index dc7ec52..36834b5 100644
--- a/examples/rotate_and_derotate_a_square.py
+++ b/examples/rotate_and_derotate_a_square.py
@@ -2,7 +2,7 @@
 import numpy as np
 
 from derotation.derotate_by_line import derotate_an_image_array_line_by_line
-from derotation.simulate.basic_rotator import Rotator
+from derotation.simulate.line_scanning_microscope import Rotator
 
 #  make a simple image, a square in a black background
 image = np.empty((100, 100))
diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
index 0aa3f69..8a4e796 100644
--- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
+++ b/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
@@ -53,7 +53,7 @@
     IncrementalPipeline,
 )
 from derotation.derotate_by_line import derotate_an_image_array_line_by_line
-from derotation.simulate.basic_rotator import Rotator
+from derotation.simulate.line_scanning_microscope import Rotator
 
 #  -----------------------------------------------------
 #  Prepare the 3D image stack and the rotation angles
diff --git a/tests/test_regression/recreate_target/shared.py b/tests/test_regression/recreate_target/shared.py
index e7ef854..d4430a8 100644
--- a/tests/test_regression/recreate_target/shared.py
+++ b/tests/test_regression/recreate_target/shared.py
@@ -5,7 +5,7 @@
 from PIL import Image
 
 from derotation.derotate_by_line import derotate_an_image_array_line_by_line
-from derotation.simulate.basic_rotator import Rotator
+from derotation.simulate.line_scanning_microscope import Rotator
 
 NUMBER_OF_FRAMES = 3
 ROTATED_IMAGES_PATH = Path("tests/test_regression/images/rotator")
diff --git a/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py b/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py
index 2b7f5f6..e8ca73a 100644
--- a/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py
+++ b/tests/test_regression/test_basic_rotator_with_different_center_of_rotation.py
@@ -3,7 +3,7 @@
 from assertions import compare_images
 from PIL import Image
 
-from derotation.simulate.basic_rotator import Rotator
+from derotation.simulate.line_scanning_microscope import Rotator
 from tests.test_regression.recreate_target.shared import (
     ROTATED_IMAGES_PATH,
     center_formatting,
diff --git a/tests/test_unit/test_rotator.py b/tests/test_unit/test_rotator.py
index 4879e97..6f5f822 100644
--- a/tests/test_unit/test_rotator.py
+++ b/tests/test_unit/test_rotator.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from derotation.simulate.basic_rotator import Rotator
+from derotation.simulate.line_scanning_microscope import Rotator
 
 
 def test_Rotator_constructor():

From 6521f4ee21a33c00446e3c987508f3e91273c36c Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 14:22:29 +0000
Subject: [PATCH 12/24] Making the script almost a test (not complete)

---
 .../test_integration/rotation_out_of_plane.py | 162 ++++++++++++------
 1 file changed, 105 insertions(+), 57 deletions(-)

diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/rotation_out_of_plane.py
index d69c610..e12ce40 100644
--- a/tests/test_integration/rotation_out_of_plane.py
+++ b/tests/test_integration/rotation_out_of_plane.py
@@ -6,74 +6,122 @@
     create_sample_image_with_two_cells,
 )
 
-from derotation.simulate.basic_rotator import Rotator
+from derotation.simulate.line_scanning_microscope import Rotator
 
-plane_angle = 25
-num_frames = 50
-pad = 20
 
-#  create a sample image with two cells
-cells = create_sample_image_with_two_cells(lines_per_frame=100)
-cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
-cells[cells == 0] = 80
+def setup():
+    plane_angle = 25
+    num_frames = 50
+    pad = 20
 
+    #  create a sample image with two cells
+    cells = create_sample_image_with_two_cells(lines_per_frame=100)
+    cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
+    cells[cells == 0] = 80
 
-image_stack = create_image_stack(cells, num_frames=num_frames)
-print("Image stack shape:", image_stack.shape)
+    image_stack = create_image_stack(cells, num_frames=num_frames)
+    print("Image stack shape:", image_stack.shape)
 
-_, angles = create_rotation_angles(image_stack.shape)
-print("Angles shape:", angles.shape)
+    _, angles = create_rotation_angles(image_stack.shape)
+    print("Angles shape:", angles.shape)
 
-#  rotate the image stack
-rotator = Rotator(angles, image_stack, rotation_plane_angle=plane_angle)
-rotated_image_stack = rotator.rotate_by_line()
+    #  rotate the image stack
+    rotator = Rotator(angles, image_stack, rotation_plane_angle=plane_angle)
+    rotated_image_stack = rotator.rotate_by_line()
+
+    print("Rotated image stack shape:", rotated_image_stack.shape)
+
+    return image_stack, rotated_image_stack, rotator, num_frames
 
-print("Rotated image stack shape:", rotated_image_stack.shape)
 
 #  plot the original image, the rotated images
-row_n = 5
-fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 5))
-
-print(ax.shape)
-ax[0, 0].imshow(image_stack[0], cmap="gray", vmin=0, vmax=255)
-ax[0, 0].set_title("Original image")
-ax[0, 0].axis("off")
-
-for n in range(1, len(rotated_image_stack) + 1):
-    row = n % row_n
-    col = n // row_n
-    ax[row, col].imshow(
-        rotated_image_stack[n - 1], cmap="gray", vmin=0, vmax=255
+def make_plot(image_stack, rotated_image_stack, rotator, num_frames):
+    row_n = 5
+    fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 25))
+
+    print(ax.shape)
+    ax[0, 0].imshow(image_stack[0], cmap="gray", vmin=0, vmax=255)
+    ax[0, 0].set_title("Original image")
+    ax[0, 0].axis("off")
+
+    for n in range(1, len(rotated_image_stack) + 1):
+        row = n % row_n
+        col = n // row_n
+        ax[row, col].imshow(
+            rotated_image_stack[n - 1], cmap="gray", vmin=0, vmax=255
+        )
+        ax[row, col].set_title(n)
+
+        # add information about the rotation angle
+        angles = rotator.angles[n - 1]
+        angle_range = f"{angles.min():.0f}-{angles.max():.0f}"
+        ax[row, col].text(
+            0.5,
+            0.9,
+            angle_range,
+            horizontalalignment="center",
+            verticalalignment="center",
+            transform=ax[row, col].transAxes,
+            color="white",
+        )
+
+    #  last one, max projection
+    ax[row_n - 1, num_frames // row_n].imshow(
+        rotated_image_stack.max(axis=0), cmap="gray"
     )
-    ax[row, col].set_title(n)
-
-    # add information about the rotation angle
-    angles = rotator.angles[n - 1]
-    angle_range = f"{angles.min():.0f}-{angles.max():.0f}"
-    ax[row, col].text(
-        0.5,
-        0.9,
-        angle_range,
-        horizontalalignment="center",
-        verticalalignment="center",
-        transform=ax[row, col].transAxes,
-        color="white",
+    #  plot also an x in the center of the image in red
+    ax[row_n - 1, num_frames // row_n].plot(
+        rotated_image_stack.shape[2] / 2,
+        rotated_image_stack.shape[1] / 2,
+        "rx",
+        markersize=10,
     )
+    ax[row_n - 1, num_frames // row_n].set_title("Max projection")
 
-#  last one, max projection
-ax[row_n - 1, num_frames // row_n].imshow(
-    rotated_image_stack.max(axis=0), cmap="gray"
-)
-#  plot also an x in the center of the image in red
-ax[row_n - 1, num_frames // row_n].plot(
-    rotated_image_stack.shape[2] / 2,
-    rotated_image_stack.shape[1] / 2,
-    "rx",
-    markersize=10,
-)
-ax[row_n - 1, num_frames // row_n].set_title("Max projection")
+    for a in ax.ravel():
+        a.axis("off")
+
+    #  save in debug folder
+    plt.savefig("debug/rotation_out_of_plane.png")
+
+
+# def test_max_projection_with_rotation_out_of_plane():
+#     _, rotated_image_stack, *_ = setup()
+
+#     max_projection = rotated_image_stack.max(axis=0)
+#     edges = canny(max_projection, sigma=7)
+
+#     fig, ax = plt.subplots(1, 3, figsize=(10, 5))
+#     ax[0].imshow(max_projection, cmap="gray")
+#     ax[0].set_title("Max projection")
+#     ax[0].axis("off")
+
+#     #  plot edges
+#     ax[1].imshow(edges, cmap="gray")
+#     ax[1].set_title("Edges")
+#     ax[1].axis("off")
+
+#     result = hough_ellipse(edges, accuracy=100, threshold=1,
+#                            min_size=0, max_size=200)
+#     result.sort(order="accumulator")
+
+#     #  get the best ellipse
+#     best = result[-1]
+#     yc, xc, a, b = [int(x) for x in best[1:5]]
+
+#     #  plot the ellipse
+#     cy, cx = ellipse_perimeter(yc, xc, a, b)
+#     max_projection[cy, cx] = 255
+
+#     #  show ellipse
+#     ax[2].imshow(max_projection, cmap="gray")
+#     ax[2].set_title("Ellipse")
+#     ax[2].axis("off")
+
+#     plt.savefig("debug/rotation_out_of_plane_max_projection.png")
 
-for a in ax.ravel():
-    a.axis("off")
 
-plt.show()
+if __name__ == "__main__":
+    image_stack, rotated_image_stack, rotator, num_frames = setup()
+    make_plot(image_stack, rotated_image_stack, rotator, num_frames)
+    # test_max_projection_with_rotation_out_of_plane()

From 8398f3bff0aa8173cad8cde96c77c513f2d19051 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 14:26:15 +0000
Subject: [PATCH 13/24] Bug fix

---
 derotation/simulate/line_scanning_microscope.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py
index fe324e0..ee80784 100644
--- a/derotation/simulate/line_scanning_microscope.py
+++ b/derotation/simulate/line_scanning_microscope.py
@@ -79,6 +79,8 @@ def __init__(
 
         if blank_pixel_val is None:
             self.blank_pixel_val = self.get_blank_pixels_value()
+        else:
+            self.blank_pixel_val = blank_pixel_val
 
     def create_homography_matrices(self) -> None:
         #  from the scanning plane to the rotation plane

From c6adc72ff7d4389ce10aeb674ee3312d920627c7 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 15:10:40 +0000
Subject: [PATCH 14/24] Add feature to handle rotation out of plane not
 vertically aligned

---
 .../simulate/line_scanning_microscope.py      | 37 ++++++++++++++-----
 1 file changed, 27 insertions(+), 10 deletions(-)

diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py
index ee80784..99e7c8c 100644
--- a/derotation/simulate/line_scanning_microscope.py
+++ b/derotation/simulate/line_scanning_microscope.py
@@ -12,6 +12,7 @@ def __init__(
         image_stack: np.ndarray,
         center: Optional[Tuple[int, int]] = None,
         rotation_plane_angle: Optional[float] = None,
+        rotation_plane_orientation: Optional[float] = None,
         blank_pixel_val: Optional[float] = None,
     ) -> None:
         """Initializes the Rotator object.
@@ -65,6 +66,10 @@ def __init__(
             self.image_size = image_stack.shape[1]
         else:
             self.rotation_plane_angle = rotation_plane_angle
+            if rotation_plane_orientation is not None:
+                self.rotation_plane_orientation = rotation_plane_orientation
+            else:
+                self.rotation_plane_orientation = 0
 
             self.create_homography_matrices()
             print(f"Pixel shift: {self.ps}")
@@ -155,8 +160,19 @@ def rotate_by_line(self) -> np.ndarray:
                     if angle == 0:
                         rotated_image_stack[i][j] = self.crop_image(image)[j]
                     else:
-                        rotated_frame = self.rotate_sample(image, angle)
-                        rotated_image_stack[i][j] = rotated_frame[j]
+                        # rotate the whole image by the angle
+                        rotated_image = self.rotate_sample(image, angle)
+
+                        # if the rotation plane angle is not 0,
+                        # apply the homography
+                        if self.rotation_plane_angle != 0:
+                            rotated_image = self.apply_homography(
+                                rotated_image, "rotation_to_scanning_plane"
+                            )
+                            rotated_image = self.crop_image(rotated_image)
+
+                        # store the rotated image line
+                        rotated_image_stack[i][j] = rotated_image[j]
             else:
                 rotated_image_stack[i] = self.crop_image(image)
 
@@ -178,7 +194,7 @@ def apply_homography(
             )
         elif direction == "rotation_to_scanning_plane":
             # forward transformation
-            return affine_transform(
+            image = affine_transform(
                 image,
                 self.inverse_homography_matrix,
                 offset=self.center - self.center,
@@ -188,6 +204,14 @@ def apply_homography(
                 cval=self.get_blank_pixels_value(),
             )
 
+            #  rotate the image back to the scanning plane angle
+            if self.rotation_plane_orientation != 0:
+                image = self.rotate_sample(
+                    self.rotation_plane_orientation, image
+                )
+
+            return image
+
     def rotate_sample(self, image: np.ndarray, angle: float) -> np.ndarray:
         """Rotate the entire image by a given angle. Uses affine transformation
         with no interpolation.
@@ -207,7 +231,6 @@ def rotate_sample(self, image: np.ndarray, angle: float) -> np.ndarray:
         np.ndarray
             The rotated image.
         """
-
         # Compute rotation in radians
         angle_rad = np.deg2rad(angle)
         cos, sin = np.cos(angle_rad), np.sin(angle_rad)
@@ -234,12 +257,6 @@ def rotate_sample(self, image: np.ndarray, angle: float) -> np.ndarray:
             cval=self.blank_pixel_val,
         )
 
-        if self.rotation_plane_angle != 0:
-            rotated_image = self.apply_homography(
-                rotated_image, "rotation_to_scanning_plane"
-            )
-            rotated_image = self.crop_image(rotated_image)
-
         return rotated_image
 
     def get_blank_pixels_value(self) -> float:

From 4ec1af0a85ed28b45cea0bf20eef27c6eebc3e68 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 15:53:10 +0000
Subject: [PATCH 15/24] Add working tests

---
 .../test_integration/rotation_out_of_plane.py | 160 ++++++++++++++----
 1 file changed, 127 insertions(+), 33 deletions(-)

diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/rotation_out_of_plane.py
index e12ce40..9bbbd7b 100644
--- a/tests/test_integration/rotation_out_of_plane.py
+++ b/tests/test_integration/rotation_out_of_plane.py
@@ -1,5 +1,8 @@
 import numpy as np
 from matplotlib import pyplot as plt
+from skimage.draw import ellipse_perimeter
+from skimage.feature import canny
+from skimage.transform import hough_ellipse
 from test_finding_center_of_rotation_by_joining_two_pipelines import (
     create_image_stack,
     create_rotation_angles,
@@ -9,7 +12,7 @@
 from derotation.simulate.line_scanning_microscope import Rotator
 
 
-def setup():
+def rotate_with_homography():
     plane_angle = 25
     num_frames = 50
     pad = 20
@@ -26,7 +29,44 @@ def setup():
     print("Angles shape:", angles.shape)
 
     #  rotate the image stack
-    rotator = Rotator(angles, image_stack, rotation_plane_angle=plane_angle)
+    rotator = Rotator(
+        angles,
+        image_stack,
+        rotation_plane_angle=plane_angle,
+        blank_pixel_val=0,
+    )
+    rotated_image_stack = rotator.rotate_by_line()
+
+    print("Rotated image stack shape:", rotated_image_stack.shape)
+
+    return image_stack, rotated_image_stack, rotator, num_frames
+
+
+def rotate_with_homography_and_orientation():
+    plane_angle = 25  # in degrees
+    num_frames = 50
+    pad = 20
+    orientation = 45  # in degrees
+
+    #  create a sample image with two cells
+    cells = create_sample_image_with_two_cells(lines_per_frame=100)
+    cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
+    cells[cells == 0] = 80
+
+    image_stack = create_image_stack(cells, num_frames=num_frames)
+    print("Image stack shape:", image_stack.shape)
+
+    _, angles = create_rotation_angles(image_stack.shape)
+    print("Angles shape:", angles.shape)
+
+    #  rotate the image stack
+    rotator = Rotator(
+        angles,
+        image_stack,
+        rotation_plane_angle=plane_angle,
+        blank_pixel_val=0,
+        rotation_plane_orientation=orientation,
+    )
     rotated_image_stack = rotator.rotate_by_line()
 
     print("Rotated image stack shape:", rotated_image_stack.shape)
@@ -35,7 +75,7 @@ def setup():
 
 
 #  plot the original image, the rotated images
-def make_plot(image_stack, rotated_image_stack, rotator, num_frames):
+def make_plot(image_stack, rotated_image_stack, rotator, num_frames, title=""):
     row_n = 5
     fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 25))
 
@@ -82,46 +122,100 @@ def make_plot(image_stack, rotated_image_stack, rotator, num_frames):
         a.axis("off")
 
     #  save in debug folder
-    plt.savefig("debug/rotation_out_of_plane.png")
+    plt.savefig(f"debug/{title}.png")
+
 
+def test_max_projection_with_rotation_out_of_plane():
+    _, rotated_image_stack, *_ = rotate_with_homography()
 
-# def test_max_projection_with_rotation_out_of_plane():
-#     _, rotated_image_stack, *_ = setup()
+    max_projection = rotated_image_stack.max(axis=0)
+    edges = canny(max_projection, sigma=7, low_threshold=35, high_threshold=50)
 
-#     max_projection = rotated_image_stack.max(axis=0)
-#     edges = canny(max_projection, sigma=7)
+    result = hough_ellipse(
+        edges, accuracy=10, threshold=100, min_size=30, max_size=100
+    )
+    result.sort(order="accumulator")
+
+    #  get the best ellipse
+    best = result[-1]
+    yc, xc, a, b, o = [int(best[i]) for i in range(1, len(best))]
+
+    assert o < 5, "Orientation should be close to 0"
 
-#     fig, ax = plt.subplots(1, 3, figsize=(10, 5))
-#     ax[0].imshow(max_projection, cmap="gray")
-#     ax[0].set_title("Max projection")
-#     ax[0].axis("off")
 
-#     #  plot edges
-#     ax[1].imshow(edges, cmap="gray")
-#     ax[1].set_title("Edges")
-#     ax[1].axis("off")
+def test_max_projection_with_rotation_out_of_plane_plus_orientation():
+    _, rotated_image_stack, *_ = rotate_with_homography_and_orientation()
 
-#     result = hough_ellipse(edges, accuracy=100, threshold=1,
-#                            min_size=0, max_size=200)
-#     result.sort(order="accumulator")
+    max_projection = rotated_image_stack.max(axis=0)
+    edges = canny(max_projection, sigma=7, low_threshold=35, high_threshold=50)
+
+    result = hough_ellipse(
+        edges, accuracy=10, threshold=100, min_size=30, max_size=100
+    )
+    result.sort(order="accumulator")
 
-#     #  get the best ellipse
-#     best = result[-1]
-#     yc, xc, a, b = [int(x) for x in best[1:5]]
+    #  get the best ellipse
+    best = result[-1]
+    yc, xc, a, b, o = [int(best[i]) for i in range(1, len(best))]
 
-#     #  plot the ellipse
-#     cy, cx = ellipse_perimeter(yc, xc, a, b)
-#     max_projection[cy, cx] = 255
+    assert np.allclose(o, 45, atol=5), "Orientation should be close to 45"
 
-#     #  show ellipse
-#     ax[2].imshow(max_projection, cmap="gray")
-#     ax[2].set_title("Ellipse")
-#     ax[2].axis("off")
 
-#     plt.savefig("debug/rotation_out_of_plane_max_projection.png")
+def plot_max_projection_with_rotation_out_of_plane(
+    max_projection, edges, yc, xc, a, b, o
+):
+    #  plot for debugging purposes
+
+    fig, ax = plt.subplots(1, 3, figsize=(10, 5))
+    ax[0].imshow(max_projection, cmap="gray")
+    ax[0].set_title("Max projection")
+    ax[0].axis("off")
+
+    #  plot edges
+    ax[1].imshow(edges, cmap="gray")
+    ax[1].set_title("Edges")
+    ax[1].axis("off")
+
+    print(f"Ellipse center: ({yc}, {xc}), a: {a}, b: {b}, orientation: {o}")
+
+    #  plot the ellipse
+    cy, cx = ellipse_perimeter(yc, xc, a, b, orientation=o)
+    empty_image = np.zeros_like(max_projection)
+    empty_image[cy, cx] = 255
+
+    #  show ellipse
+    ax[2].imshow(empty_image, cmap="gray")
+    ax[2].set_title("Ellipse")
+    ax[2].axis("off")
+
+    plt.savefig("debug/rotation_out_of_plane_max_projection.png")
 
 
 if __name__ == "__main__":
-    image_stack, rotated_image_stack, rotator, num_frames = setup()
-    make_plot(image_stack, rotated_image_stack, rotator, num_frames)
-    # test_max_projection_with_rotation_out_of_plane()
+    (
+        image_stack,
+        rotated_image_stack,
+        rotator,
+        num_frames,
+    ) = rotate_with_homography()
+    make_plot(
+        image_stack,
+        rotated_image_stack,
+        rotator,
+        num_frames,
+        title="rotation_out_of_plane",
+    )
+
+    (
+        image_stack,
+        rotated_image_stack,
+        rotator,
+        num_frames,
+    ) = rotate_with_homography_and_orientation()
+    make_plot(
+        image_stack,
+        rotated_image_stack,
+        rotator,
+        num_frames,
+        title="rotation_out_of_plane_plus_orientation",
+    )

From 53140ed7cddfac89a94e3e09f8843efca0167d8a Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 15:53:28 +0000
Subject: [PATCH 16/24] Update gitignore

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index f43115e..e1ba14f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -83,3 +83,4 @@ venv/
 **/_version.py
 
 .conda/
+debug/*

From 4bb3779ff01af6d886959f0719355e8931add2b7 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 15:56:29 +0000
Subject: [PATCH 17/24] Refactor tests

---
 .../test_integration/rotation_out_of_plane.py | 199 +++++++-----------
 1 file changed, 81 insertions(+), 118 deletions(-)

diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/rotation_out_of_plane.py
index 9bbbd7b..d9e68b0 100644
--- a/tests/test_integration/rotation_out_of_plane.py
+++ b/tests/test_integration/rotation_out_of_plane.py
@@ -1,6 +1,7 @@
+from typing import Optional, Tuple
+
 import numpy as np
 from matplotlib import pyplot as plt
-from skimage.draw import ellipse_perimeter
 from skimage.feature import canny
 from skimage.transform import hough_ellipse
 from test_finding_center_of_rotation_by_joining_two_pipelines import (
@@ -12,54 +13,40 @@
 from derotation.simulate.line_scanning_microscope import Rotator
 
 
-def rotate_with_homography():
-    plane_angle = 25
-    num_frames = 50
-    pad = 20
-
-    #  create a sample image with two cells
-    cells = create_sample_image_with_two_cells(lines_per_frame=100)
-    cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
-    cells[cells == 0] = 80
-
-    image_stack = create_image_stack(cells, num_frames=num_frames)
-    print("Image stack shape:", image_stack.shape)
-
-    _, angles = create_rotation_angles(image_stack.shape)
-    print("Angles shape:", angles.shape)
-
-    #  rotate the image stack
-    rotator = Rotator(
-        angles,
-        image_stack,
-        rotation_plane_angle=plane_angle,
-        blank_pixel_val=0,
-    )
-    rotated_image_stack = rotator.rotate_by_line()
-
-    print("Rotated image stack shape:", rotated_image_stack.shape)
-
-    return image_stack, rotated_image_stack, rotator, num_frames
-
-
-def rotate_with_homography_and_orientation():
-    plane_angle = 25  # in degrees
-    num_frames = 50
-    pad = 20
-    orientation = 45  # in degrees
-
-    #  create a sample image with two cells
+def rotate_image_stack(
+    plane_angle: int,
+    num_frames: int,
+    pad: int,
+    orientation: Optional[int] = None,
+) -> Tuple[np.ndarray, np.ndarray, Rotator, int]:
+    """
+    Rotates an image stack with a specified plane angle and
+    optional orientation.
+
+    Parameters
+    ----------
+    plane_angle : int
+        The angle of the rotation plane in degrees.
+    num_frames : int
+        Number of frames in the image stack.
+    pad : int
+        Padding size to apply to the sample image.
+    orientation : int, optional
+        Orientation of the rotation plane in degrees, by default None.
+
+    Returns
+    -------
+    Tuple[np.ndarray, np.ndarray, Rotator, int]
+        Original image stack, rotated image stack, Rotator instance,
+        and number of frames.
+    """
     cells = create_sample_image_with_two_cells(lines_per_frame=100)
     cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
     cells[cells == 0] = 80
 
     image_stack = create_image_stack(cells, num_frames=num_frames)
-    print("Image stack shape:", image_stack.shape)
-
     _, angles = create_rotation_angles(image_stack.shape)
-    print("Angles shape:", angles.shape)
 
-    #  rotate the image stack
     rotator = Rotator(
         angles,
         image_stack,
@@ -68,18 +55,35 @@ def rotate_with_homography_and_orientation():
         rotation_plane_orientation=orientation,
     )
     rotated_image_stack = rotator.rotate_by_line()
-
-    print("Rotated image stack shape:", rotated_image_stack.shape)
-
     return image_stack, rotated_image_stack, rotator, num_frames
 
 
-#  plot the original image, the rotated images
-def make_plot(image_stack, rotated_image_stack, rotator, num_frames, title=""):
+def make_plot(
+    image_stack: np.ndarray,
+    rotated_image_stack: np.ndarray,
+    rotator: Rotator,
+    num_frames: int,
+    title: str = "",
+) -> None:
+    """
+    Plots the original and rotated image stacks.
+
+    Parameters
+    ----------
+    image_stack : np.ndarray
+        Original image stack.
+    rotated_image_stack : np.ndarray
+        Rotated image stack.
+    rotator : Rotator
+        Rotator instance containing rotation metadata.
+    num_frames : int
+        Number of frames in the image stack.
+    title : str, optional
+        Title for the saved plot, by default "".
+    """
     row_n = 5
     fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 25))
 
-    print(ax.shape)
     ax[0, 0].imshow(image_stack[0], cmap="gray", vmin=0, vmax=255)
     ax[0, 0].set_title("Original image")
     ax[0, 0].axis("off")
@@ -92,7 +96,6 @@ def make_plot(image_stack, rotated_image_stack, rotator, num_frames, title=""):
         )
         ax[row, col].set_title(n)
 
-        # add information about the rotation angle
         angles = rotator.angles[n - 1]
         angle_range = f"{angles.min():.0f}-{angles.max():.0f}"
         ax[row, col].text(
@@ -105,11 +108,9 @@ def make_plot(image_stack, rotated_image_stack, rotator, num_frames, title=""):
             color="white",
         )
 
-    #  last one, max projection
     ax[row_n - 1, num_frames // row_n].imshow(
         rotated_image_stack.max(axis=0), cmap="gray"
     )
-    #  plot also an x in the center of the image in red
     ax[row_n - 1, num_frames // row_n].plot(
         rotated_image_stack.shape[2] / 2,
         rotated_image_stack.shape[1] / 2,
@@ -121,83 +122,48 @@ def make_plot(image_stack, rotated_image_stack, rotator, num_frames, title=""):
     for a in ax.ravel():
         a.axis("off")
 
-    #  save in debug folder
     plt.savefig(f"debug/{title}.png")
 
 
-def test_max_projection_with_rotation_out_of_plane():
-    _, rotated_image_stack, *_ = rotate_with_homography()
-
-    max_projection = rotated_image_stack.max(axis=0)
-    edges = canny(max_projection, sigma=7, low_threshold=35, high_threshold=50)
-
-    result = hough_ellipse(
-        edges, accuracy=10, threshold=100, min_size=30, max_size=100
-    )
-    result.sort(order="accumulator")
-
-    #  get the best ellipse
-    best = result[-1]
-    yc, xc, a, b, o = [int(best[i]) for i in range(1, len(best))]
-
-    assert o < 5, "Orientation should be close to 0"
-
-
-def test_max_projection_with_rotation_out_of_plane_plus_orientation():
-    _, rotated_image_stack, *_ = rotate_with_homography_and_orientation()
-
+def test_max_projection(
+    rotated_image_stack: np.ndarray,
+    expected_orientation: Optional[int] = None,
+    atol: int = 5,
+) -> None:
+    """
+    Tests the max projection for an image stack to validate rotation.
+
+    Parameters
+    ----------
+    rotated_image_stack : np.ndarray
+        Rotated image stack.
+    expected_orientation : int, optional
+        Expected orientation of the rotation plane, by default None.
+    atol : int, optional
+        Allowed tolerance for orientation, by default 5.
+    """
     max_projection = rotated_image_stack.max(axis=0)
     edges = canny(max_projection, sigma=7, low_threshold=35, high_threshold=50)
-
     result = hough_ellipse(
         edges, accuracy=10, threshold=100, min_size=30, max_size=100
     )
     result.sort(order="accumulator")
 
-    #  get the best ellipse
     best = result[-1]
     yc, xc, a, b, o = [int(best[i]) for i in range(1, len(best))]
 
-    assert np.allclose(o, 45, atol=5), "Orientation should be close to 45"
-
-
-def plot_max_projection_with_rotation_out_of_plane(
-    max_projection, edges, yc, xc, a, b, o
-):
-    #  plot for debugging purposes
-
-    fig, ax = plt.subplots(1, 3, figsize=(10, 5))
-    ax[0].imshow(max_projection, cmap="gray")
-    ax[0].set_title("Max projection")
-    ax[0].axis("off")
-
-    #  plot edges
-    ax[1].imshow(edges, cmap="gray")
-    ax[1].set_title("Edges")
-    ax[1].axis("off")
-
-    print(f"Ellipse center: ({yc}, {xc}), a: {a}, b: {b}, orientation: {o}")
-
-    #  plot the ellipse
-    cy, cx = ellipse_perimeter(yc, xc, a, b, orientation=o)
-    empty_image = np.zeros_like(max_projection)
-    empty_image[cy, cx] = 255
-
-    #  show ellipse
-    ax[2].imshow(empty_image, cmap="gray")
-    ax[2].set_title("Ellipse")
-    ax[2].axis("off")
-
-    plt.savefig("debug/rotation_out_of_plane_max_projection.png")
+    if expected_orientation is not None:
+        assert np.allclose(
+            o, expected_orientation, atol=atol
+        ), f"Orientation should be close to {expected_orientation}"
+    else:
+        assert o < 5, "Orientation should be close to 0"
 
 
 if __name__ == "__main__":
-    (
-        image_stack,
-        rotated_image_stack,
-        rotator,
-        num_frames,
-    ) = rotate_with_homography()
+    image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
+        plane_angle=25, num_frames=50, pad=20
+    )
     make_plot(
         image_stack,
         rotated_image_stack,
@@ -206,12 +172,9 @@ def plot_max_projection_with_rotation_out_of_plane(
         title="rotation_out_of_plane",
     )
 
-    (
-        image_stack,
-        rotated_image_stack,
-        rotator,
-        num_frames,
-    ) = rotate_with_homography_and_orientation()
+    image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
+        plane_angle=25, num_frames=50, pad=20, orientation=45
+    )
     make_plot(
         image_stack,
         rotated_image_stack,

From fd19ca11e7f1d5c8435152ea042599f3dde224c1 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 17:56:13 +0000
Subject: [PATCH 18/24] Bug fix

---
 derotation/simulate/line_scanning_microscope.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py
index 99e7c8c..00476d3 100644
--- a/derotation/simulate/line_scanning_microscope.py
+++ b/derotation/simulate/line_scanning_microscope.py
@@ -52,6 +52,12 @@ def __init__(
             + f"({image_stack.shape[0] * image_stack.shape[1]})"
         )
 
+        if rotation_plane_orientation is not None:
+            assert rotation_plane_angle is not None, (
+                "If rotation_plane_orientation is provided, "
+                + "rotation_plane_angle should be provided as well."
+            )
+
         self.image_stack = image_stack
         self.num_lines_per_frame = image_stack.shape[1]
 
@@ -207,7 +213,7 @@ def apply_homography(
             #  rotate the image back to the scanning plane angle
             if self.rotation_plane_orientation != 0:
                 image = self.rotate_sample(
-                    self.rotation_plane_orientation, image
+                    image, self.rotation_plane_orientation
                 )
 
             return image

From f301d96a1da91a3f7d859ae0a4bea8dfc97a7b12 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Thu, 5 Dec 2024 17:57:16 +0000
Subject: [PATCH 19/24] =?UTF-8?q?WIP=20tests=20-=20breaking=20=F0=9F=90=9B?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 ...plane.py => test_rotation_out_of_plane.py} | 147 +++++++++++++++---
 1 file changed, 127 insertions(+), 20 deletions(-)
 rename tests/test_integration/{rotation_out_of_plane.py => test_rotation_out_of_plane.py} (51%)

diff --git a/tests/test_integration/rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py
similarity index 51%
rename from tests/test_integration/rotation_out_of_plane.py
rename to tests/test_integration/test_rotation_out_of_plane.py
index d9e68b0..e951cf5 100644
--- a/tests/test_integration/rotation_out_of_plane.py
+++ b/tests/test_integration/test_rotation_out_of_plane.py
@@ -1,22 +1,24 @@
 from typing import Optional, Tuple
 
 import numpy as np
+import pytest
 from matplotlib import pyplot as plt
+from skimage.draw import ellipse_perimeter
 from skimage.feature import canny
 from skimage.transform import hough_ellipse
-from test_finding_center_of_rotation_by_joining_two_pipelines import (
+
+from derotation.simulate.line_scanning_microscope import Rotator
+from tests.test_integration.test_finding_center_of_rotation_by_joining_two_pipelines import (
     create_image_stack,
     create_rotation_angles,
     create_sample_image_with_two_cells,
 )
 
-from derotation.simulate.line_scanning_microscope import Rotator
-
 
 def rotate_image_stack(
-    plane_angle: int,
-    num_frames: int,
-    pad: int,
+    plane_angle: int = 0,
+    num_frames: int = 50,
+    pad: int = 20,
     orientation: Optional[int] = None,
 ) -> Tuple[np.ndarray, np.ndarray, Rotator, int]:
     """
@@ -125,42 +127,147 @@ def make_plot(
     plt.savefig(f"debug/{title}.png")
 
 
+@pytest.mark.parametrize(
+    "plane_angle,expected_orientation",
+    [
+        (0, None),
+        (15, None),
+        (25, None),
+        (25, 45),
+        (25, 90),
+        (25, 135),
+        (40, None),
+        (40, 45),
+        (40, 90),
+        (40, 135),
+    ],
+)
 def test_max_projection(
-    rotated_image_stack: np.ndarray,
-    expected_orientation: Optional[int] = None,
-    atol: int = 5,
+    plane_angle: int,
+    expected_orientation: Optional[int],
+    atol: int = 15,
 ) -> None:
     """
-    Tests the max projection for an image stack to validate rotation.
+    Tests if the max projection of an image stack fits an ellipse,
+    and verify orientation and major/minor axes are close to expected values:
+    - if plane_angle is 0, the major and minor axes should be close to each other;
+    - if expected_orientation is not None, the orientation should be close to it.
+    - if expected_orientation is None, the orientation should be close to 0.
+
+    Giving by default an allowed tolerance of 15 degrees as I care about
+    the test tracking the general shape of the ellipse, not the exact values.
 
     Parameters
     ----------
-    rotated_image_stack : np.ndarray
-        Rotated image stack.
+    plane_angle : int
+        The angle of the rotation plane in degrees.
     expected_orientation : int, optional
         Expected orientation of the rotation plane, by default None.
     atol : int, optional
         Allowed tolerance for orientation, by default 5.
     """
+    _, rotated_image_stack, *_ = rotate_image_stack(
+        plane_angle=plane_angle, orientation=expected_orientation
+    )
+
     max_projection = rotated_image_stack.max(axis=0)
-    edges = canny(max_projection, sigma=7, low_threshold=35, high_threshold=50)
+    edges = canny(max_projection, sigma=7, low_threshold=30, high_threshold=50)
     result = hough_ellipse(
-        edges, accuracy=10, threshold=100, min_size=30, max_size=100
+        edges, accuracy=9, threshold=110, min_size=35, max_size=100
     )
     result.sort(order="accumulator")
 
     best = result[-1]
-    yc, xc, a, b, o = [int(best[i]) for i in range(1, len(best))]
+    yc, xc, a, b, orientation = [int(best[i]) for i in range(1, len(best))]
 
-    if expected_orientation is not None:
-        assert np.allclose(
-            o, expected_orientation, atol=atol
-        ), f"Orientation should be close to {expected_orientation}"
+    #  orientation: Major axis orientation in clockwise direction as radians.
+    #  expected orientation is calculated frm the minor axis, so convert
+    orientation = np.abs(np.rad2deg(orientation))
+
+    if expected_orientation is None and plane_angle == 0:
+        #  lower tolerance for the major and minor axes
+        if not np.allclose(a, b, atol=atol - 5):
+            plot_max_projection_with_rotation_out_of_plane(
+                max_projection,
+                edges,
+                yc,
+                xc,
+                a,
+                b,
+                orientation,
+                f"debug/error_{plane_angle}_{expected_orientation}.png",
+            )
+            assert (
+                False
+            ), f"Major and minor axes should be close to each other, instead got {a} and {b}"
+    elif expected_orientation is not None:
+        if not np.allclose(orientation, expected_orientation, atol=atol):
+            plot_max_projection_with_rotation_out_of_plane(
+                max_projection,
+                edges,
+                yc,
+                xc,
+                a,
+                b,
+                orientation,
+                f"debug/error_{plane_angle}_{expected_orientation}.png",
+            )
+            assert (
+                False
+            ), f"Orientation should be close to {expected_orientation}, instead got {orientation}"
     else:
-        assert o < 5, "Orientation should be close to 0"
+        # check that the major axis (a) differs from the minor axis (b) as calculated with cosine
+        decrease = np.cos(np.deg2rad(plane_angle))
+        expected_b = a - a * decrease
+        #  lower tolerance for the major and minor axes
+        if np.allclose(b, expected_b, atol=atol - 5):
+            plot_max_projection_with_rotation_out_of_plane(
+                max_projection,
+                edges,
+                yc,
+                xc,
+                a,
+                b,
+                orientation,
+                f"debug/error_{plane_angle}_{expected_orientation}.png",
+            )
+            assert (
+                False
+            ), f"Difference between major and minor axes should be close to {expected_b}, instead got {b}"
+
+
+def plot_max_projection_with_rotation_out_of_plane(
+    max_projection, edges, yc, xc, a, b, o, title
+):
+    #  plot for debugging purposes
+
+    fig, ax = plt.subplots(1, 3, figsize=(10, 5))
+    ax[0].imshow(max_projection, cmap="gray")
+    ax[0].set_title("Max projection")
+    ax[0].axis("off")
+
+    #  plot edges
+    ax[1].imshow(edges, cmap="gray")
+    ax[1].set_title("Edges")
+    ax[1].axis("off")
+
+    print(f"Ellipse center: ({yc}, {xc}), a: {a}, b: {b}, orientation: {o}")
+
+    #  plot the ellipse
+    cy, cx = ellipse_perimeter(yc, xc, a, b, orientation=o)
+    empty_image = np.zeros_like(max_projection)
+    empty_image[cy, cx] = 255
+
+    #  show ellipse
+    ax[2].imshow(empty_image, cmap="gray")
+    ax[2].set_title("Ellipse")
+    ax[2].axis("off")
+
+    plt.savefig(title)
 
 
 if __name__ == "__main__":
+    #  To be used for debugging purposes
     image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
         plane_angle=25, num_frames=50, pad=20
     )

From 1f33526b908b1d53ed4834aecf9b54a15f6a33a2 Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 9 Dec 2024 18:19:12 +0000
Subject: [PATCH 20/24] =?UTF-8?q?Add=20working=20tests=20=E2=9C=A8=20reusi?=
 =?UTF-8?q?ng=20ellipse=20fit=20function=20created=20before=20-=20required?=
 =?UTF-8?q?=20refactoring?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 derotation/analysis/fit_ellipse.py            | 159 ++++++++++++++++
 .../incremental_derotation_pipeline.py        | 170 ++----------------
 ....py => test_finding_center_of_rotation.py} |  20 +--
 .../test_rotation_out_of_plane.py             | 147 ++++++++-------
 4 files changed, 267 insertions(+), 229 deletions(-)
 create mode 100644 derotation/analysis/fit_ellipse.py
 rename tests/test_integration/{test_finding_center_of_rotation_by_joining_two_pipelines.py => test_finding_center_of_rotation.py} (98%)

diff --git a/derotation/analysis/fit_ellipse.py b/derotation/analysis/fit_ellipse.py
new file mode 100644
index 0000000..46ec9a9
--- /dev/null
+++ b/derotation/analysis/fit_ellipse.py
@@ -0,0 +1,159 @@
+import logging
+from pathlib import Path
+from typing import Tuple
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.patches import Ellipse
+from scipy.optimize import OptimizeResult, least_squares
+
+
+def fit_ellipse_to_points(
+    centers: np.ndarray,
+) -> Tuple[int, int, int, int, int]:
+    """Fit an ellipse to the points using least squares optimization.
+
+    Parameters
+    ----------
+    centers : np.ndarray
+        The centers of the largest blob in each image.
+
+    Returns
+    -------
+    Tuple[int, int, int, int, int]
+        The center of the ellipse (center_x, center_y), the semi-major
+        axis (a), the semi-minor axis (b), and the rotation angle (theta).
+    """
+    # Convert centers to numpy array
+    centers = np.array(centers)
+    x = centers[:, 0]
+    y = centers[:, 1]
+
+    # 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)]
+
+    # Initial parameters: (center_x, center_y, semi_major_axis,
+    # semi_minor_axis, rotation_angle)
+    initial_center = np.mean(
+        [topmost, bottommost, leftmost, rightmost], axis=0
+    )
+    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
+
+    initial_params = [
+        initial_center[0],
+        initial_center[1],
+        semi_major_axis,
+        semi_minor_axis,
+        rotation_angle,
+    ]
+
+    logging.info("Fitting ellipse to points...")
+
+    # Objective function to minimize: sum of squared distances to ellipse
+    def ellipse_residuals(params, x, y):
+        center_x, center_y, a, b, theta = params  # theta is in radians
+        cos_angle = np.cos(theta)
+        sin_angle = np.sin(theta)
+
+        # Rotate the points to align with the ellipse axes
+        x_rot = cos_angle * (x - center_x) + sin_angle * (y - center_y)
+        y_rot = -sin_angle * (x - center_x) + cos_angle * (y - center_y)
+
+        # Ellipse equation: (x_rot^2 / a^2) + (y_rot^2 / b^2) = 1
+        return (x_rot / a) ** 2 + (y_rot / b) ** 2 - 1
+
+    # Use least squares optimization to fit the ellipse to the points
+    result: OptimizeResult = least_squares(
+        ellipse_residuals,
+        initial_params,
+        args=(x, y),
+        loss="huber",  # minimize the influence of outliers
+    )
+
+    # Extract optimized parameters
+    center_x, center_y, a, b, theta = result.x
+
+    return center_x, center_y, a, b, theta
+
+
+def plot_ellipse_fit_and_centers(
+    centers: np.ndarray,
+    center_x: int,
+    center_y: int,
+    a: int,
+    b: int,
+    theta: int,
+    image_stack: np.ndarray,
+    debug_plots_folder: Path,
+    saving_name: str = "ellipse_fit.png",
+):
+    """Plot the fitted ellipse on the largest blob centers.
+
+    Parameters
+    ----------
+    centers : np.ndarray
+        The centers of the largest blob in each image.
+    center_x : int
+        The x-coordinate of the center of the ellipse
+    center_y : int
+        The y-coordinate of the center of the ellipse
+    a : int
+        The semi-major axis of the ellipse
+    b : int
+        The semi-minor axis of the ellipse
+    theta : int
+        The rotation angle of the ellipse
+    """
+    # Convert centers to numpy array
+    centers = np.array(centers)
+    x = centers[:, 0]
+    y = centers[:, 1]
+
+    # Create the plot
+    fig, ax = plt.subplots(figsize=(8, 8))
+
+    max_projection = image_stack.max(axis=0)
+    #  plot behind a frame of the original image
+    ax.imshow(max_projection, cmap="gray")
+    ax.scatter(x, y, label="Largest Blob Centers", color="red")
+
+    # Plot fitted ellipse
+    ellipse = Ellipse(
+        (center_x, center_y),
+        width=2 * a,
+        height=2 * b,
+        angle=np.degrees(theta),
+        edgecolor="blue",
+        facecolor="none",
+        label="Fitted Ellipse",
+    )
+    ax.add_patch(ellipse)
+
+    # Plot center of fitted ellipse
+    ax.scatter(
+        center_x,
+        center_y,
+        color="green",
+        marker="x",
+        s=100,
+        label="Ellipse Center",
+    )
+
+    # Add some plot formatting
+    ax.set_xlim(0, image_stack.shape[1])
+    ax.set_ylim(
+        image_stack.shape[1], 0
+    )  # Invert y-axis to match image coordinate system
+    ax.set_aspect("equal")
+    ax.legend()
+    ax.grid(True)
+    ax.set_title("Fitted Ellipse on largest blob centers")
+    ax.axis("off")
+
+    #  save
+    plt.savefig(debug_plots_folder / saving_name)
diff --git a/derotation/analysis/incremental_derotation_pipeline.py b/derotation/analysis/incremental_derotation_pipeline.py
index 2723282..f0f8873 100644
--- a/derotation/analysis/incremental_derotation_pipeline.py
+++ b/derotation/analysis/incremental_derotation_pipeline.py
@@ -6,12 +6,14 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import pandas as pd
-from matplotlib.patches import Ellipse
 from scipy.ndimage import rotate
-from scipy.optimize import OptimizeResult, least_squares
 from skimage.feature import blob_log
 from tqdm import tqdm
 
+from derotation.analysis.fit_ellipse import (
+    fit_ellipse_to_points,
+    plot_ellipse_fit_and_centers,
+)
 from derotation.analysis.full_derotation_pipeline import FullPipeline
 
 
@@ -489,10 +491,23 @@ def find_center_of_rotation(self) -> Tuple[int, int]:
         )
 
         # Fit an ellipse to the largest blob centers and get its center
-        center_x, center_y, a, b, theta = self.fit_ellipse_to_points(
+        center_x, center_y, a, b, theta = fit_ellipse_to_points(
             coord_first_blob_of_every_image
         )
 
+        if self.debugging_plots:
+            plot_ellipse_fit_and_centers(
+                coord_first_blob_of_every_image,
+                center_x,
+                center_y,
+                a,
+                b,
+                theta,
+                image_stack=mean_images,
+                debug_plots_folder=self.debug_plots_folder,
+                saving_name="ellipse_fit.png",
+            )
+
         logging.info(
             f"Center of ellipse: ({center_x:.2f}, {center_y:.2f}), "
             f"semi-major axis: {a:.2f}, semi-minor axis: {b:.2f}"
@@ -530,7 +545,7 @@ def get_coords_of_largest_blob(
         blobs = [
             blobs[i][blobs[i][:, 2].argsort()] for i in range(len(image_stack))
         ]
-        
+
         try:
             coord_first_blob_of_every_image = [
                 blobs[i][0][:2].astype(int) for i in range(len(image_stack))
@@ -579,150 +594,3 @@ def plot_blob_detection(self, blobs: list, image_stack: np.ndarray):
 
         # save the plot
         plt.savefig(self.debug_plots_folder / "blobs.png")
-
-    def fit_ellipse_to_points(
-        self, centers: np.ndarray
-    ) -> Tuple[int, int, int, int, int]:
-        """Fit an ellipse to the points using least squares optimization.
-
-        Parameters
-        ----------
-        centers : np.ndarray
-            The centers of the largest blob in each image.
-
-        Returns
-        -------
-        Tuple[int, int, int, int, int]
-            The center of the ellipse (center_x, center_y), the semi-major
-            axis (a), the semi-minor axis (b), and the rotation angle (theta).
-        """
-        # Convert centers to numpy array
-        centers = np.array(centers)
-        x = centers[:, 0]
-        y = centers[:, 1]
-
-        # 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)]
-
-        # Initial parameters: (center_x, center_y, semi_major_axis,
-        # semi_minor_axis, rotation_angle)
-        initial_center = np.mean(
-            [topmost, bottommost, leftmost, rightmost], axis=0
-        )
-        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
-
-        initial_params = [
-            initial_center[0],
-            initial_center[1],
-            semi_major_axis,
-            semi_minor_axis,
-            rotation_angle,
-        ]
-
-        logging.info("Fitting ellipse to points...")
-
-        # Objective function to minimize: sum of squared distances to ellipse
-        def ellipse_residuals(params, x, y):
-            center_x, center_y, a, b, theta = params
-            cos_angle = np.cos(theta)
-            sin_angle = np.sin(theta)
-
-            # Rotate the points to align with the ellipse axes
-            x_rot = cos_angle * (x - center_x) + sin_angle * (y - center_y)
-            y_rot = -sin_angle * (x - center_x) + cos_angle * (y - center_y)
-
-            # Ellipse equation: (x_rot^2 / a^2) + (y_rot^2 / b^2) = 1
-            return (x_rot / a) ** 2 + (y_rot / b) ** 2 - 1
-
-        # Use least squares optimization to fit the ellipse to the points
-        result: OptimizeResult = least_squares(
-            ellipse_residuals, initial_params, args=(x, y)
-        )
-
-        # Extract optimized parameters
-        center_x, center_y, a, b, theta = result.x
-
-        if self.debugging_plots:
-            self.plot_ellipse_fit_and_centers(
-                centers, center_x, center_y, a, b, theta
-            )
-
-        return center_x, center_y, a, b, theta
-
-    def plot_ellipse_fit_and_centers(
-        self,
-        centers: np.ndarray,
-        center_x: int,
-        center_y: int,
-        a: int,
-        b: int,
-        theta: int,
-    ):
-        """Plot the fitted ellipse on the largest blob centers.
-
-        Parameters
-        ----------
-        centers : np.ndarray
-            The centers of the largest blob in each image.
-        center_x : int
-            The x-coordinate of the center of the ellipse
-        center_y : int
-            The y-coordinate of the center of the ellipse
-        a : int
-            The semi-major axis of the ellipse
-        b : int
-            The semi-minor axis of the ellipse
-        theta : int
-            The rotation angle of the ellipse
-        """
-        # Convert centers to numpy array
-        centers = np.array(centers)
-        x = centers[:, 0]
-        y = centers[:, 1]
-
-        # Create the plot
-        fig, ax = plt.subplots(figsize=(8, 8))
-        #  plot behind a frame of the original image
-        ax.imshow(self.image_stack[0], cmap="gray")
-        ax.scatter(x, y, label="Largest Blob Centers", color="red")
-
-        # Plot fitted ellipse
-        ellipse = Ellipse(
-            (center_x, center_y),
-            width=2 * a,
-            height=2 * b,
-            angle=np.degrees(theta),
-            edgecolor="blue",
-            facecolor="none",
-            label="Fitted Ellipse",
-        )
-        ax.add_patch(ellipse)
-
-        # Plot center of fitted ellipse
-        ax.scatter(
-            center_x,
-            center_y,
-            color="green",
-            marker="x",
-            s=100,
-            label="Ellipse Center",
-        )
-
-        # Add some plot formatting
-        ax.set_xlim(0, self.image_stack.shape[1])
-        ax.set_ylim(
-            self.image_stack.shape[1], 0
-        )  # Invert y-axis to match image coordinate system
-        ax.set_aspect("equal")
-        ax.legend()
-        ax.grid(True)
-        ax.set_title("Fitted Ellipse on largest blob centers")
-        ax.axis("off")
-
-        #  save
-        plt.savefig(self.debug_plots_folder / "ellipse_fit.png")
diff --git a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py b/tests/test_integration/test_finding_center_of_rotation.py
similarity index 98%
rename from tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
rename to tests/test_integration/test_finding_center_of_rotation.py
index 8a4e796..46c38d1 100644
--- a/tests/test_integration/test_finding_center_of_rotation_by_joining_two_pipelines.py
+++ b/tests/test_integration/test_finding_center_of_rotation.py
@@ -64,6 +64,8 @@ def create_sample_image_with_two_cells(
     center_of_bright_cell: Tuple[int, int] = (50, 10),
     center_of_dimmer_cell: Tuple[int, int] = (60, 60),
     lines_per_frame: int = 100,
+    second_cell=True,
+    radius: int = 5,
 ) -> np.ndarray:
     """Create a 2D image with two circles, one bright and one dim
     by default in the top center and bottom right, respectively.
@@ -88,13 +90,8 @@ def create_sample_image_with_two_cells(
     image = np.zeros((lines_per_frame, lines_per_frame), dtype=np.uint8)
 
     # Define the circle's parameters
-    radius = 5  # radius of the circle
     white_value = 255  # white color for the circle
 
-    #  add an extra gray circle at the bottom right
-    radius2 = 5
-    gray_value = 128
-
     # Draw a white circle in the top center
     y, x = np.ogrid[: image.shape[0], : image.shape[1]]
     mask = (x - center_of_bright_cell[0]) ** 2 + (
@@ -102,11 +99,14 @@ def create_sample_image_with_two_cells(
     ) ** 2 <= radius**2
     image[mask] = white_value
 
-    # Draw a gray circle in the bottom right
-    mask2 = (x - center_of_dimmer_cell[0]) ** 2 + (
-        y - center_of_dimmer_cell[1]
-    ) ** 2 <= radius2**2
-    image[mask2] = gray_value
+    if second_cell:
+        #  add an extra gray circle at the bottom right
+        gray_value = 128
+        # Draw a gray circle in the bottom right
+        mask2 = (x - center_of_dimmer_cell[0]) ** 2 + (
+            y - center_of_dimmer_cell[1]
+        ) ** 2 <= radius**2
+        image[mask2] = gray_value
 
     return image
 
diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py
index e951cf5..78c7d35 100644
--- a/tests/test_integration/test_rotation_out_of_plane.py
+++ b/tests/test_integration/test_rotation_out_of_plane.py
@@ -1,14 +1,17 @@
+from pathlib import Path
 from typing import Optional, Tuple
 
 import numpy as np
 import pytest
 from matplotlib import pyplot as plt
 from skimage.draw import ellipse_perimeter
-from skimage.feature import canny
-from skimage.transform import hough_ellipse
 
+from derotation.analysis.fit_ellipse import (
+    fit_ellipse_to_points,
+    plot_ellipse_fit_and_centers,
+)
 from derotation.simulate.line_scanning_microscope import Rotator
-from tests.test_integration.test_finding_center_of_rotation_by_joining_two_pipelines import (
+from tests.test_integration.test_finding_center_of_rotation import (
     create_image_stack,
     create_rotation_angles,
     create_sample_image_with_two_cells,
@@ -42,7 +45,9 @@ def rotate_image_stack(
         Original image stack, rotated image stack, Rotator instance,
         and number of frames.
     """
-    cells = create_sample_image_with_two_cells(lines_per_frame=100)
+    cells = create_sample_image_with_two_cells(
+        lines_per_frame=100, second_cell=False, radius=1
+    )
     cells = np.pad(cells, ((pad, pad), (pad, pad)), mode="constant")
     cells[cells == 0] = 80
 
@@ -135,23 +140,23 @@ def make_plot(
         (25, None),
         (25, 45),
         (25, 90),
-        (25, 135),
         (40, None),
         (40, 45),
         (40, 90),
-        (40, 135),
     ],
 )
 def test_max_projection(
     plane_angle: int,
-    expected_orientation: Optional[int],
+    exp_orientation: Optional[int],
     atol: int = 15,
 ) -> None:
     """
     Tests if the max projection of an image stack fits an ellipse,
     and verify orientation and major/minor axes are close to expected values:
-    - if plane_angle is 0, the major and minor axes should be close to each other;
-    - if expected_orientation is not None, the orientation should be close to it.
+    - if plane_angle is 0, the major and minor axes should be close to each
+      other;
+    - if expected_orientation is not None, the orientation should be close to
+      it.
     - if expected_orientation is None, the orientation should be close to 0.
 
     Giving by default an allowed tolerance of 15 degrees as I care about
@@ -167,89 +172,95 @@ def test_max_projection(
         Allowed tolerance for orientation, by default 5.
     """
     _, rotated_image_stack, *_ = rotate_image_stack(
-        plane_angle=plane_angle, orientation=expected_orientation
+        plane_angle=plane_angle, orientation=exp_orientation
     )
 
-    max_projection = rotated_image_stack.max(axis=0)
-    edges = canny(max_projection, sigma=7, low_threshold=30, high_threshold=50)
-    result = hough_ellipse(
-        edges, accuracy=9, threshold=110, min_size=35, max_size=100
+    # get the location of the brightest pixel in any frame as a center
+    centers = np.array(
+        [
+            np.unravel_index(np.argmax(frame), frame.shape)
+            for frame in rotated_image_stack
+        ]
     )
-    result.sort(order="accumulator")
-
-    best = result[-1]
-    yc, xc, a, b, orientation = [int(best[i]) for i in range(1, len(best))]
+    #  invert centers to (x, y) format
+    centers = np.array([(x, y) for y, x in centers])
 
-    #  orientation: Major axis orientation in clockwise direction as radians.
-    #  expected orientation is calculated frm the minor axis, so convert
-    orientation = np.abs(np.rad2deg(orientation))
+    xc, yc, a, b, orientation = fit_ellipse_to_points(centers)
 
-    if expected_orientation is None and plane_angle == 0:
+    if exp_orientation is None and plane_angle == 0:
         #  lower tolerance for the major and minor axes
         if not np.allclose(a, b, atol=atol - 5):
-            plot_max_projection_with_rotation_out_of_plane(
-                max_projection,
-                edges,
-                yc,
-                xc,
-                a,
-                b,
-                orientation,
-                f"debug/error_{plane_angle}_{expected_orientation}.png",
+            plot_ellipse_fit_and_centers(
+                image_stack=rotated_image_stack,
+                centers=centers,
+                center_x=xc,
+                center_y=yc,
+                a=a,
+                b=b,
+                theta=orientation,
+                debug_plots_folder=Path("debug/"),
+                saving_name=f"ellipse_fit_{plane_angle}_{exp_orientation}.png",
             )
             assert (
                 False
-            ), f"Major and minor axes should be close to each other, instead got {a} and {b}"
-    elif expected_orientation is not None:
-        if not np.allclose(orientation, expected_orientation, atol=atol):
-            plot_max_projection_with_rotation_out_of_plane(
-                max_projection,
-                edges,
-                yc,
-                xc,
-                a,
-                b,
-                orientation,
-                f"debug/error_{plane_angle}_{expected_orientation}.png",
+            ), f"Major and minor axes should be close, instead got {a} and {b}"
+    elif exp_orientation is not None:
+        #  Major axis orientation in clockwise direction as radians.
+        exp_orientation = np.deg2rad(exp_orientation)
+        if not np.allclose(orientation, exp_orientation, atol=atol):
+            plot_ellipse_fit_and_centers(
+                image_stack=rotated_image_stack,
+                centers=centers,
+                center_x=xc,
+                center_y=yc,
+                a=a,
+                b=b,
+                theta=orientation,
+                debug_plots_folder=Path("debug/"),
+                saving_name=f"ellipse_fit_{plane_angle}_{exp_orientation}.png",
+            )
+            assert False, (
+                f"Orientation should be close to {exp_orientation}, "
+                f"instead got {orientation}"
             )
-            assert (
-                False
-            ), f"Orientation should be close to {expected_orientation}, instead got {orientation}"
     else:
-        # check that the major axis (a) differs from the minor axis (b) as calculated with cosine
+        # check that the major axis (a) differs from the minor axis (b)
+        # as calculated with cosine
         decrease = np.cos(np.deg2rad(plane_angle))
         expected_b = a - a * decrease
         #  lower tolerance for the major and minor axes
         if np.allclose(b, expected_b, atol=atol - 5):
-            plot_max_projection_with_rotation_out_of_plane(
-                max_projection,
-                edges,
-                yc,
-                xc,
-                a,
-                b,
-                orientation,
-                f"debug/error_{plane_angle}_{expected_orientation}.png",
+            plot_ellipse_fit_and_centers(
+                image_stack=rotated_image_stack,
+                centers=centers,
+                center_x=xc,
+                center_y=yc,
+                a=a,
+                b=b,
+                theta=orientation,
+                debug_plots_folder=Path("debug/"),
+                saving_name=f"ellipse_fit_{plane_angle}_{exp_orientation}.png",
+            )
+            assert False, (
+                f"Difference between major and minor axes should be close to "
+                f"{expected_b}, instead got {b}"
             )
-            assert (
-                False
-            ), f"Difference between major and minor axes should be close to {expected_b}, instead got {b}"
 
 
 def plot_max_projection_with_rotation_out_of_plane(
-    max_projection, edges, yc, xc, a, b, o, title
+    max_projection, yc, xc, a, b, o, title
 ):
     #  plot for debugging purposes
 
-    fig, ax = plt.subplots(1, 3, figsize=(10, 5))
+    fig, ax = plt.subplots(1, 2, figsize=(10, 5))
     ax[0].imshow(max_projection, cmap="gray")
     ax[0].set_title("Max projection")
     ax[0].axis("off")
 
-    #  plot edges
-    ax[1].imshow(edges, cmap="gray")
-    ax[1].set_title("Edges")
-    ax[1].axis("off")
+    # #  plot edges
+    # ax[1].imshow(edges, cmap="gray")
+    # ax[1].set_title("Edges")
+    # ax[1].axis("off")
 
     print(f"Ellipse center: ({yc}, {xc}), a: {a}, b: {b}, orientation: {o}")
 
@@ -259,9 +270,9 @@ def plot_max_projection_with_rotation_out_of_plane(
     empty_image[cy, cx] = 255
 
     #  show ellipse
-    ax[2].imshow(empty_image, cmap="gray")
-    ax[2].set_title("Ellipse")
-    ax[2].axis("off")
+    ax[1].imshow(empty_image, cmap="gray")
+    ax[1].set_title("Fitted ellipse")
+    ax[1].axis("off")
 
     plt.savefig(title)
 

From 5b32532be68910e5bbc488f0ec3692e8c9ae0bce Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 9 Dec 2024 18:47:42 +0000
Subject: [PATCH 21/24] Add mode doscrings and some refactoring

---
 derotation/analysis/fit_ellipse.py            |   6 +
 .../simulate/line_scanning_microscope.py      |  68 +++++++++-
 examples/elliptical_rotations.py              | 100 ++++++++++++++
 .../test_finding_center_of_rotation.py        |   8 +-
 .../test_rotation_out_of_plane.py             | 124 ------------------
 5 files changed, 178 insertions(+), 128 deletions(-)
 create mode 100644 examples/elliptical_rotations.py

diff --git a/derotation/analysis/fit_ellipse.py b/derotation/analysis/fit_ellipse.py
index 46ec9a9..f92a887 100644
--- a/derotation/analysis/fit_ellipse.py
+++ b/derotation/analysis/fit_ellipse.py
@@ -108,6 +108,12 @@ def plot_ellipse_fit_and_centers(
         The semi-minor axis of the ellipse
     theta : int
         The rotation angle of the ellipse
+    image_stack : np.ndarray
+        The image stack to plot the ellipse on.
+    debug_plots_folder : Path
+        The folder to save the debug plot to.
+    saving_name : str, optional
+        The name of the file to save the plot to, by default "ellipse_fit.png".
     """
     # Convert centers to numpy array
     centers = np.array(centers)
diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py
index 00476d3..e94ce3d 100644
--- a/derotation/simulate/line_scanning_microscope.py
+++ b/derotation/simulate/line_scanning_microscope.py
@@ -21,6 +21,11 @@ def __init__(
         the acquisition of a given line as if it was instantaneous, happening
         while the sample was rotated at a given angle.
 
+        It is also possible to simulate the acquisition of a movie from a
+        rotation plane that differs from the scanning plane. To achieve this,
+        provide the rotation_plane_angle and if you want the orientation as
+        well.
+
         The purpouse of the Rotator object is to imitate the acquisition of
         rotated samples in order to validate the derotation algorithms and in
         the future, build a forward model of the transformation.
@@ -38,12 +43,25 @@ def __init__(
             num_lines_per_frame, num_pixels_per_line). In case you want to
             rotate a single frame, provide an (1, num_lines_per_frame,
             num_pixels_per_line) image stack.
-
+        center : Tuple[int, int], optional
+            The center of rotation. If None, the center is going to be the
+            center of the image, by default None
+        rotation_plane_angle : Optional[float], optional
+            The z angle of the rotation plane in degrees in relation to the
+            scanning plane. If None or 0, the rotation plane is the same as
+            the  scanning plane, by default None.
+        rotation_plane_orientation : Optional[float], optional
+            The angle of the rotation plane in the x-y plane in degrees,
+            transposed into the rotation plane. If None or 0, the rotation
+            plane is the same as the scanning plane, by default None.
         Raises
         ------
-        AssertionError
+        AssertionError (1)
             If the number of angles is not equal to the number of lines
             per frame multiplied by the number of frames
+        AssertionError (2)
+            If rotation_plane_orientation is provided, but rotation_plane_angle
+            is not provided.
         """
         #  there should be one angle per line pe frame
         assert len(angles) == image_stack.shape[0] * image_stack.shape[1], (
@@ -53,6 +71,9 @@ def __init__(
         )
 
         if rotation_plane_orientation is not None:
+            # The rotation plane angle makes sense only if the orientation is
+            # provided, as without it the rotation is going to be circular
+            # and not elliptical.
             assert rotation_plane_angle is not None, (
                 "If rotation_plane_orientation is provided, "
                 + "rotation_plane_angle should be provided as well."
@@ -78,6 +99,7 @@ def __init__(
                 self.rotation_plane_orientation = 0
 
             self.create_homography_matrices()
+            self.calculate_pixel_shift()
             print(f"Pixel shift: {self.ps}")
             print(f"New image size: {self.image_size}")
 
@@ -94,6 +116,22 @@ def __init__(
             self.blank_pixel_val = blank_pixel_val
 
     def create_homography_matrices(self) -> None:
+        """
+        Create the homography matrices to simulate the acquisition of a
+        rotated sample from a different plane than the scanning plane.
+
+        The homography matrix is used to transform the image from the scanning
+        plane to the rotation plane and vice-versa.
+
+        The homography matrix is defined as:
+        H = [[1, 0, 0],
+             [0, cos(theta), 0],
+             [0, 0, 1]]
+        where theta is the rotation plane angle in degrees.
+
+        Currently, we are using only the inverse homography matrix to transform
+        the image from the rotation plane to the scanning plane.
+        """
         #  from the scanning plane to the rotation plane
         self.homography_matrix = np.array(
             [
@@ -106,6 +144,11 @@ def create_homography_matrices(self) -> None:
         #  from the rotation plane to the scanning plane
         self.inverse_homography_matrix = np.linalg.inv(self.homography_matrix)
 
+    def calculate_pixel_shift(self) -> None:
+        """
+        Calculate the pixel shift and the new image size based on the rotation
+        plane angle.
+        """
         #  store pixels shift based on inverse homography
         line_length = self.image_stack.shape[2]
         self.ps = (
@@ -126,6 +169,7 @@ def create_homography_matrices(self) -> None:
         self.image_size = line_length - self.ps
 
     def crop_image(self, image: np.ndarray) -> np.ndarray:
+        """Crop the image to the new size based on the pixel shift."""
         if self.ps == 0:
             return image
         else:
@@ -138,7 +182,8 @@ def crop_image(self, image: np.ndarray) -> np.ndarray:
 
     def rotate_by_line(self) -> np.ndarray:
         """Simulate the acquisition of a rotated image stack as if for each
-        line acquired, the sample was rotated at a given angle.
+        line acquired, the sample was rotated at a given angle in a given
+        center and plane of rotation.
 
         Each frame is rotated n_lines_per_frame times, where n_lines_per_frame
         is the number of lines per frame in the image stack.
@@ -187,6 +232,23 @@ def rotate_by_line(self) -> np.ndarray:
     def apply_homography(
         self, image: np.ndarray, direction: str
     ) -> np.ndarray:
+        """Apply the homography matrix to the image to simulate the acquisition
+        of a rotated sample from a different plane than the scanning plane.
+
+        Parameters
+        ----------
+        image : np.ndarray
+            The image to apply the homography.
+        direction : str
+            The direction of the transformation. It can be either
+            "scanning_to_rotation_plane" or "rotation_to_scanning_plane".
+
+        Returns
+        -------
+        np.ndarray
+            The transformed image.
+        """
+
         if direction == "scanning_to_rotation_plane":
             # backward transformation
             return affine_transform(
diff --git a/examples/elliptical_rotations.py b/examples/elliptical_rotations.py
new file mode 100644
index 0000000..975de05
--- /dev/null
+++ b/examples/elliptical_rotations.py
@@ -0,0 +1,100 @@
+#  Visualize the acquisition of a movie where the rotation axis is not
+#  aligned with the image plane.
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from derotation.simulate.line_scanning_microscope import Rotator
+from tests.test_integration.test_rotation_out_of_plane import (
+    rotate_image_stack,
+)
+
+
+def make_plot(
+    image_stack: np.ndarray,
+    rotated_image_stack: np.ndarray,
+    rotator: Rotator,
+    num_frames: int,
+    title: str = "",
+) -> None:
+    """
+    Make figure with rotated image stacks frame by frame.
+
+    Parameters
+    ----------
+    image_stack : np.ndarray
+        Original image stack.
+    rotated_image_stack : np.ndarray
+        Rotated image stack.
+    rotator : Rotator
+        Rotator instance containing rotation metadata.
+    num_frames : int
+        Number of frames in the image stack.
+    title : str, optional
+        Title for the saved plot, by default "".
+    """
+    row_n = 5
+    fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 25))
+
+    ax[0, 0].imshow(image_stack[0], cmap="gray", vmin=0, vmax=255)
+    ax[0, 0].set_title("Original image")
+    ax[0, 0].axis("off")
+
+    for n in range(1, len(rotated_image_stack) + 1):
+        row = n % row_n
+        col = n // row_n
+        ax[row, col].imshow(
+            rotated_image_stack[n - 1], cmap="gray", vmin=0, vmax=255
+        )
+        ax[row, col].set_title(n)
+
+        angles = rotator.angles[n - 1]
+        angle_range = f"{angles.min():.0f}-{angles.max():.0f}"
+        ax[row, col].text(
+            0.5,
+            0.9,
+            angle_range,
+            horizontalalignment="center",
+            verticalalignment="center",
+            transform=ax[row, col].transAxes,
+            color="white",
+        )
+
+    ax[row_n - 1, num_frames // row_n].imshow(
+        rotated_image_stack.max(axis=0), cmap="gray"
+    )
+    ax[row_n - 1, num_frames // row_n].plot(
+        rotated_image_stack.shape[2] / 2,
+        rotated_image_stack.shape[1] / 2,
+        "rx",
+        markersize=10,
+    )
+    ax[row_n - 1, num_frames // row_n].set_title("Max projection")
+
+    for a in ax.ravel():
+        a.axis("off")
+
+    plt.savefig(f"debug/{title}.png")
+
+
+image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
+    plane_angle=25, num_frames=50, pad=20
+)
+make_plot(
+    image_stack,
+    rotated_image_stack,
+    rotator,
+    num_frames,
+    title="rotation_out_of_plane",
+)
+
+image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
+    plane_angle=25, num_frames=50, pad=20, orientation=45
+)
+make_plot(
+    image_stack,
+    rotated_image_stack,
+    rotator,
+    num_frames,
+    title="rotation_out_of_plane_plus_orientation",
+)
diff --git a/tests/test_integration/test_finding_center_of_rotation.py b/tests/test_integration/test_finding_center_of_rotation.py
index 46c38d1..4d229ef 100644
--- a/tests/test_integration/test_finding_center_of_rotation.py
+++ b/tests/test_integration/test_finding_center_of_rotation.py
@@ -67,7 +67,7 @@ def create_sample_image_with_two_cells(
     second_cell=True,
     radius: int = 5,
 ) -> np.ndarray:
-    """Create a 2D image with two circles, one bright and one dim
+    """Create a 2D image with two circles, one bright and one dim (optional)
     by default in the top center and bottom right, respectively.
 
     Location of the circles can be changed by providing the
@@ -79,6 +79,12 @@ def create_sample_image_with_two_cells(
         Location of brightest cell, by default (50, 10)
     center_of_dimmer_cell : Tuple[int, int], optional
         Location of dimmer cell, by default (60, 60)
+    lines_per_frame : int, optional
+        Number of lines per frame, by default 100
+    second_cell : bool, optional
+        Add an extra dimmer cell, by default True
+    radius : int, optional
+        Radius of the circles, by default 5
 
     Returns
     -------
diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py
index 78c7d35..ec40248 100644
--- a/tests/test_integration/test_rotation_out_of_plane.py
+++ b/tests/test_integration/test_rotation_out_of_plane.py
@@ -3,8 +3,6 @@
 
 import numpy as np
 import pytest
-from matplotlib import pyplot as plt
-from skimage.draw import ellipse_perimeter
 
 from derotation.analysis.fit_ellipse import (
     fit_ellipse_to_points,
@@ -65,73 +63,6 @@ def rotate_image_stack(
     return image_stack, rotated_image_stack, rotator, num_frames
 
 
-def make_plot(
-    image_stack: np.ndarray,
-    rotated_image_stack: np.ndarray,
-    rotator: Rotator,
-    num_frames: int,
-    title: str = "",
-) -> None:
-    """
-    Plots the original and rotated image stacks.
-
-    Parameters
-    ----------
-    image_stack : np.ndarray
-        Original image stack.
-    rotated_image_stack : np.ndarray
-        Rotated image stack.
-    rotator : Rotator
-        Rotator instance containing rotation metadata.
-    num_frames : int
-        Number of frames in the image stack.
-    title : str, optional
-        Title for the saved plot, by default "".
-    """
-    row_n = 5
-    fig, ax = plt.subplots(row_n, num_frames // row_n + 1, figsize=(40, 25))
-
-    ax[0, 0].imshow(image_stack[0], cmap="gray", vmin=0, vmax=255)
-    ax[0, 0].set_title("Original image")
-    ax[0, 0].axis("off")
-
-    for n in range(1, len(rotated_image_stack) + 1):
-        row = n % row_n
-        col = n // row_n
-        ax[row, col].imshow(
-            rotated_image_stack[n - 1], cmap="gray", vmin=0, vmax=255
-        )
-        ax[row, col].set_title(n)
-
-        angles = rotator.angles[n - 1]
-        angle_range = f"{angles.min():.0f}-{angles.max():.0f}"
-        ax[row, col].text(
-            0.5,
-            0.9,
-            angle_range,
-            horizontalalignment="center",
-            verticalalignment="center",
-            transform=ax[row, col].transAxes,
-            color="white",
-        )
-
-    ax[row_n - 1, num_frames // row_n].imshow(
-        rotated_image_stack.max(axis=0), cmap="gray"
-    )
-    ax[row_n - 1, num_frames // row_n].plot(
-        rotated_image_stack.shape[2] / 2,
-        rotated_image_stack.shape[1] / 2,
-        "rx",
-        markersize=10,
-    )
-    ax[row_n - 1, num_frames // row_n].set_title("Max projection")
-
-    for a in ax.ravel():
-        a.axis("off")
-
-    plt.savefig(f"debug/{title}.png")
-
-
 @pytest.mark.parametrize(
     "plane_angle,expected_orientation",
     [
@@ -245,58 +176,3 @@ def test_max_projection(
                 f"Difference between major and minor axes should be close to "
                 f"{expected_b}, instead got {b}"
             )
-
-
-def plot_max_projection_with_rotation_out_of_plane(
-    max_projection, yc, xc, a, b, o, title
-):
-    #  plot for debugging purposes
-
-    fig, ax = plt.subplots(1, 2, figsize=(10, 5))
-    ax[0].imshow(max_projection, cmap="gray")
-    ax[0].set_title("Max projection")
-    ax[0].axis("off")
-
-    # #  plot edges
-    # ax[1].imshow(edges, cmap="gray")
-    # ax[1].set_title("Edges")
-    # ax[1].axis("off")
-
-    print(f"Ellipse center: ({yc}, {xc}), a: {a}, b: {b}, orientation: {o}")
-
-    #  plot the ellipse
-    cy, cx = ellipse_perimeter(yc, xc, a, b, orientation=o)
-    empty_image = np.zeros_like(max_projection)
-    empty_image[cy, cx] = 255
-
-    #  show ellipse
-    ax[1].imshow(empty_image, cmap="gray")
-    ax[1].set_title("Fitted ellipse")
-    ax[1].axis("off")
-
-    plt.savefig(title)
-
-
-if __name__ == "__main__":
-    #  To be used for debugging purposes
-    image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
-        plane_angle=25, num_frames=50, pad=20
-    )
-    make_plot(
-        image_stack,
-        rotated_image_stack,
-        rotator,
-        num_frames,
-        title="rotation_out_of_plane",
-    )
-
-    image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
-        plane_angle=25, num_frames=50, pad=20, orientation=45
-    )
-    make_plot(
-        image_stack,
-        rotated_image_stack,
-        rotator,
-        num_frames,
-        title="rotation_out_of_plane_plus_orientation",
-    )

From 7097c744023379eede8abbd8d5463f33e4bbea7f Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 9 Dec 2024 19:07:20 +0000
Subject: [PATCH 22/24] Fix docstring

---
 tests/test_integration/test_rotation_out_of_plane.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py
index ec40248..9700a73 100644
--- a/tests/test_integration/test_rotation_out_of_plane.py
+++ b/tests/test_integration/test_rotation_out_of_plane.py
@@ -97,7 +97,7 @@ def test_max_projection(
     ----------
     plane_angle : int
         The angle of the rotation plane in degrees.
-    expected_orientation : int, optional
+    exp_orientation : int, optional
         Expected orientation of the rotation plane, by default None.
     atol : int, optional
         Allowed tolerance for orientation, by default 5.

From 4e855e8b04f4ee56dc9badebc259176b99da086a Mon Sep 17 00:00:00 2001
From: lauraporta <29216006+lauraporta@users.noreply.github.com>
Date: Mon, 9 Dec 2024 19:09:50 +0000
Subject: [PATCH 23/24] Fix tests

---
 tests/test_integration/test_rotation_out_of_plane.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py
index 9700a73..9df92a7 100644
--- a/tests/test_integration/test_rotation_out_of_plane.py
+++ b/tests/test_integration/test_rotation_out_of_plane.py
@@ -64,7 +64,7 @@ def rotate_image_stack(
 
 
 @pytest.mark.parametrize(
-    "plane_angle,expected_orientation",
+    "plane_angle,exp_orientation",
     [
         (0, None),
         (15, None),
@@ -86,9 +86,9 @@ def test_max_projection(
     and verify orientation and major/minor axes are close to expected values:
     - if plane_angle is 0, the major and minor axes should be close to each
       other;
-    - if expected_orientation is not None, the orientation should be close to
+    - if exp_orientation is not None, the orientation should be close to
       it.
-    - if expected_orientation is None, the orientation should be close to 0.
+    - if exp_orientation is None, the orientation should be close to 0.
 
     Giving by default an allowed tolerance of 15 degrees as I care about
     the test tracking the general shape of the ellipse, not the exact values.

From 4129dacdbdff4f72b9fd6ebe997e0573d7565932 Mon Sep 17 00:00:00 2001
From: lauraporta <ucqflpo@ucl.ac.uk>
Date: Tue, 17 Dec 2024 13:19:52 +0000
Subject: [PATCH 24/24] Apply suggestions from Igor

---
 .../simulate/line_scanning_microscope.py      | 142 ++++++++----------
 examples/elliptical_rotations.py              |   4 +
 .../test_rotation_out_of_plane.py             |  28 ++--
 3 files changed, 77 insertions(+), 97 deletions(-)

diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py
index e94ce3d..3298391 100644
--- a/derotation/simulate/line_scanning_microscope.py
+++ b/derotation/simulate/line_scanning_microscope.py
@@ -11,12 +11,12 @@ def __init__(
         angles: np.ndarray,
         image_stack: np.ndarray,
         center: Optional[Tuple[int, int]] = None,
-        rotation_plane_angle: Optional[float] = None,
-        rotation_plane_orientation: Optional[float] = None,
+        rotation_plane_angle: float = 0,
+        rotation_plane_orientation: float = 0,
         blank_pixel_val: Optional[float] = None,
     ) -> None:
         """Initializes the Rotator object.
-        The Rotator aims to imitate a the scanning pattern of a multi-photon
+        The Rotator aims to imitate the scanning pattern of a multi-photon
         microscope while the speciment is rotating. Currently, it approximates
         the acquisition of a given line as if it was instantaneous, happening
         while the sample was rotated at a given angle.
@@ -26,7 +26,7 @@ def __init__(
         provide the rotation_plane_angle and if you want the orientation as
         well.
 
-        The purpouse of the Rotator object is to imitate the acquisition of
+        The purpose of the Rotator object is to imitate the acquisition of
         rotated samples in order to validate the derotation algorithms and in
         the future, build a forward model of the transformation.
 
@@ -46,14 +46,17 @@ def __init__(
         center : Tuple[int, int], optional
             The center of rotation. If None, the center is going to be the
             center of the image, by default None
-        rotation_plane_angle : Optional[float], optional
+        rotation_plane_angle : float, optional
             The z angle of the rotation plane in degrees in relation to the
-            scanning plane. If None or 0, the rotation plane is the same as
+            scanning plane. If 0, the rotation plane is the same as
             the  scanning plane, by default None.
-        rotation_plane_orientation : Optional[float], optional
+        rotation_plane_orientation : float, optional
             The angle of the rotation plane in the x-y plane in degrees,
-            transposed into the rotation plane. If None or 0, the rotation
+            transposed into the rotation plane. If 0, the rotation
             plane is the same as the scanning plane, by default None.
+        blank_pixel_val : Optional[float], optional
+            The value to fill the blank pixels with. If None, it is going to be
+            the minimum value of the image stack, by default None.
         Raises
         ------
         AssertionError (1)
@@ -87,17 +90,12 @@ def __init__(
         else:
             self.center = np.array(center)
 
-        if rotation_plane_angle is None:
-            self.rotation_plane_angle: float = 0
+        self.rotation_plane_angle = np.deg2rad(rotation_plane_angle)
+        self.rotation_plane_orientation = rotation_plane_orientation
+        if rotation_plane_angle == 0:
             self.ps: int = 0
             self.image_size = image_stack.shape[1]
         else:
-            self.rotation_plane_angle = rotation_plane_angle
-            if rotation_plane_orientation is not None:
-                self.rotation_plane_orientation = rotation_plane_orientation
-            else:
-                self.rotation_plane_orientation = 0
-
             self.create_homography_matrices()
             self.calculate_pixel_shift()
             print(f"Pixel shift: {self.ps}")
@@ -136,7 +134,7 @@ def create_homography_matrices(self) -> None:
         self.homography_matrix = np.array(
             [
                 [1, 0, 0],
-                [0, np.cos(np.radians(self.rotation_plane_angle)), 0],
+                [0, np.cos(self.rotation_plane_angle), 0],
                 [0, 0, 1],
             ]
         )
@@ -151,19 +149,12 @@ def calculate_pixel_shift(self) -> None:
         """
         #  store pixels shift based on inverse homography
         line_length = self.image_stack.shape[2]
-        self.ps = (
-            line_length
-            - np.round(
-                np.abs(
-                    line_length * np.cos(np.radians(self.rotation_plane_angle))
-                )
-            ).astype(int)
-            + 1
-        )
+        self.ps = line_length - np.round(
+            np.abs(line_length * np.cos(self.rotation_plane_angle))
+        ).astype(int)
 
         #  round to the nearest even number
-        if self.ps % 2 != 0:
-            self.ps -= 1
+        self.ps += self.ps % 2
 
         #  final image size depends on the pixel shift
         self.image_size = line_length - self.ps
@@ -202,46 +193,45 @@ def rotate_by_line(self) -> np.ndarray:
         for i, image in tqdm.tqdm(
             enumerate(self.image_stack), total=self.image_stack.shape[0]
         ):
-            is_this_frame_rotating = not np.all(
+            if np.all(
                 # don't bother if rotation is less than 0.01 degrees
                 np.isclose(self.angles[i], 0, atol=1e-2)
-            )
-            if is_this_frame_rotating:
-                for j, angle in enumerate(self.angles[i]):
-                    if angle == 0:
-                        rotated_image_stack[i][j] = self.crop_image(image)[j]
-                    else:
-                        # rotate the whole image by the angle
-                        rotated_image = self.rotate_sample(image, angle)
-
-                        # if the rotation plane angle is not 0,
-                        # apply the homography
-                        if self.rotation_plane_angle != 0:
-                            rotated_image = self.apply_homography(
-                                rotated_image, "rotation_to_scanning_plane"
+            ):
+                rotated_image_stack[i] = self.crop_image(image)
+                continue
+
+            for j, angle in enumerate(self.angles[i]):
+                if angle == 0:
+                    rotated_image_stack[i][j] = self.crop_image(image)[j]
+                else:
+                    # rotate the whole image by the angle
+                    rotated_image = self.rotate_sample(image, angle)
+
+                    # if the rotation plane angle is not 0,
+                    # apply the homography
+                    if self.rotation_plane_angle != 0:
+                        rotated_image = (
+                            self.homography_rotation_to_scanning_plane(
+                                rotated_image
                             )
-                            rotated_image = self.crop_image(rotated_image)
+                        )
+                        rotated_image = self.crop_image(rotated_image)
 
-                        # store the rotated image line
-                        rotated_image_stack[i][j] = rotated_image[j]
-            else:
-                rotated_image_stack[i] = self.crop_image(image)
+                    # store the rotated image line
+                    rotated_image_stack[i][j] = rotated_image[j]
 
         return rotated_image_stack
 
-    def apply_homography(
-        self, image: np.ndarray, direction: str
+    def homography_rotation_to_scanning_plane(
+        self, image: np.ndarray
     ) -> np.ndarray:
-        """Apply the homography matrix to the image to simulate the acquisition
-        of a rotated sample from a different plane than the scanning plane.
+        """Apply the homography to the image to simulate the acquisition of a
+        rotated sample from a different plane than the scanning plane.
 
         Parameters
         ----------
         image : np.ndarray
             The image to apply the homography.
-        direction : str
-            The direction of the transformation. It can be either
-            "scanning_to_rotation_plane" or "rotation_to_scanning_plane".
 
         Returns
         -------
@@ -249,36 +239,22 @@ def apply_homography(
             The transformed image.
         """
 
-        if direction == "scanning_to_rotation_plane":
-            # backward transformation
-            return affine_transform(
-                image,
-                self.homography_matrix,
-                offset=self.center - self.center,
-                output_shape=image.shape,
-                order=0,
-                mode="constant",
-                cval=self.get_blank_pixels_value(),
-            )
-        elif direction == "rotation_to_scanning_plane":
-            # forward transformation
-            image = affine_transform(
-                image,
-                self.inverse_homography_matrix,
-                offset=self.center - self.center,
-                output_shape=image.shape,
-                order=0,
-                mode="constant",
-                cval=self.get_blank_pixels_value(),
-            )
+        # forward transformation
+        image = affine_transform(
+            image,
+            self.inverse_homography_matrix,
+            offset=self.center,
+            output_shape=image.shape,
+            order=0,
+            mode="constant",
+            cval=self.blank_pixel_val,
+        )
 
-            #  rotate the image back to the scanning plane angle
-            if self.rotation_plane_orientation != 0:
-                image = self.rotate_sample(
-                    image, self.rotation_plane_orientation
-                )
+        #  rotate the image back to the scanning plane angle
+        if self.rotation_plane_orientation != 0:
+            image = self.rotate_sample(image, self.rotation_plane_orientation)
 
-            return image
+        return image
 
     def rotate_sample(self, image: np.ndarray, angle: float) -> np.ndarray:
         """Rotate the entire image by a given angle. Uses affine transformation
diff --git a/examples/elliptical_rotations.py b/examples/elliptical_rotations.py
index 975de05..69ab8ae 100644
--- a/examples/elliptical_rotations.py
+++ b/examples/elliptical_rotations.py
@@ -1,6 +1,8 @@
 #  Visualize the acquisition of a movie where the rotation axis is not
 #  aligned with the image plane.
 
+from pathlib import Path
+
 import matplotlib.pyplot as plt
 import numpy as np
 
@@ -77,6 +79,8 @@ def make_plot(
     plt.savefig(f"debug/{title}.png")
 
 
+Path("debug").mkdir(exist_ok=True)
+
 image_stack, rotated_image_stack, rotator, num_frames = rotate_image_stack(
     plane_angle=25, num_frames=50, pad=20
 )
diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py
index 9df92a7..2a84189 100644
--- a/tests/test_integration/test_rotation_out_of_plane.py
+++ b/tests/test_integration/test_rotation_out_of_plane.py
@@ -1,5 +1,5 @@
 from pathlib import Path
-from typing import Optional, Tuple
+from typing import Tuple
 
 import numpy as np
 import pytest
@@ -17,10 +17,10 @@
 
 
 def rotate_image_stack(
-    plane_angle: int = 0,
+    plane_angle: float = 0,
     num_frames: int = 50,
     pad: int = 20,
-    orientation: Optional[int] = None,
+    orientation: float = 0,
 ) -> Tuple[np.ndarray, np.ndarray, Rotator, int]:
     """
     Rotates an image stack with a specified plane angle and
@@ -28,13 +28,13 @@ def rotate_image_stack(
 
     Parameters
     ----------
-    plane_angle : int
+    plane_angle : float, optional
         The angle of the rotation plane in degrees.
     num_frames : int
         Number of frames in the image stack.
     pad : int
         Padding size to apply to the sample image.
-    orientation : int, optional
+    orientation : float, optional
         Orientation of the rotation plane in degrees, by default None.
 
     Returns
@@ -66,19 +66,19 @@ def rotate_image_stack(
 @pytest.mark.parametrize(
     "plane_angle,exp_orientation",
     [
-        (0, None),
-        (15, None),
-        (25, None),
+        (0, 0),
+        (15, 0),
+        (25, 0),
         (25, 45),
         (25, 90),
-        (40, None),
+        (40, 0),
         (40, 45),
         (40, 90),
     ],
 )
 def test_max_projection(
-    plane_angle: int,
-    exp_orientation: Optional[int],
+    plane_angle: float,
+    exp_orientation: float,
     atol: int = 15,
 ) -> None:
     """
@@ -95,12 +95,12 @@ def test_max_projection(
 
     Parameters
     ----------
-    plane_angle : int
+    plane_angle : float
         The angle of the rotation plane in degrees.
-    exp_orientation : int, optional
+    exp_orientation : float
         Expected orientation of the rotation plane, by default None.
     atol : int, optional
-        Allowed tolerance for orientation, by default 5.
+        Allowed tolerance for orientation, by default 15.
     """
     _, rotated_image_stack, *_ = rotate_image_stack(
         plane_angle=plane_angle, orientation=exp_orientation