Add mode doscrings and some refactoring

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)

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(

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",
+)

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
     -------