Add Rotator (#25)

* Add basic rotator and a simple usage of it

* Add docstrings

* Add basic tests

* Make function name consistent

* Fix Rotator and join it with the derotator in an example

* Add tests

* Delete duplicate file

* Resolve naming (in many cases with "rotation" I meant derotation)

* Remove linear interpolation and recompute test images - now it's a square of different gray lines

* Fix angle indices bug in Rotator

* Use blank pixel

* Update workflows and add dependabot yml

* Simplify filter

* Apply suggestions from PR review

derotation/analysis/full_rotation_pipeline.py → derotation/analysis/full_derotation_pipeline.py

@@ -14,7 +14,7 @@
 from sklearn.mixture import GaussianMixture
 from tifffile import imsave
 
-from derotation.derotate_by_line import rotate_an_image_array_line_by_line
+from derotation.derotate_by_line import derotate_an_image_array_line_by_line
 from derotation.load_data.custom_data_loaders import (
     get_analog_signals,
     read_randomized_stim_table,
@@ -56,7 +56,7 @@ def __call__(self):
         - saving the masked image stack
         """
         self.process_analog_signals()
-        rotated_images = self.rotate_frames_line_by_line()
+        rotated_images = self.derotate_frames_line_by_line()
         masked = self.add_circle_mask(rotated_images, self.mask_diameter)
         self.save(masked)
         self.save_csv_with_derotation_data()
@@ -790,7 +790,7 @@ def plot_rotation_angles(self):
         plt.savefig(self.debug_plots_folder / "rotation_angles.png")
 
     ### ----------------- Derotation ----------------- ###
-    def rotate_frames_line_by_line(self) -> np.ndarray:
+    def derotate_frames_line_by_line(self) -> np.ndarray:
         """Rotates the image stack line by line, using the rotation angles
         by line calculated from the analog signals.
 
@@ -820,7 +820,7 @@ def rotate_frames_line_by_line(self) -> np.ndarray:
 
         offset = self.find_image_offset(self.image_stack[0])
 
-        rotated_image_stack = rotate_an_image_array_line_by_line(
+        rotated_image_stack = derotate_an_image_array_line_by_line(
             self.image_stack,
             self.rot_deg_line,
             blank_pixels_value=offset,

derotation/analysis/incremental_rotation_pipeline.py → derotation/analysis/incremental_derotation_pipeline.py

@@ -5,12 +5,11 @@
 
 import numpy as np
 import pandas as pd
+from full_derotation_pipeline import FullPipeline
 from matplotlib import pyplot as plt
 from scipy.ndimage import rotate
 from tqdm import tqdm
 
-from derotation.analysis.full_rotation_pipeline import FullPipeline
-
 
 class IncrementalPipeline(FullPipeline):
     """Derotate the image stack that was acquired using the incremental

derotation/derotate_by_line.py

@@ -5,7 +5,7 @@
 from scipy.ndimage import rotate
 
 
-def rotate_an_image_array_line_by_line(
+def derotate_an_image_array_line_by_line(
     image_stack: np.ndarray,
     rot_deg_line: np.ndarray,
     blank_pixels_value: float = 0,
@@ -18,7 +18,7 @@ def rotate_an_image_array_line_by_line(
     - creates a new image with only that line
     - rotates the line by the given angle
     - substitutes the line in the new image
-    - adds the new image to the rotated image stack
+    - adds the new image to the derotated image stack
 
     Edge cases and how they are handled:
     - the rotation starts in the middle of the image -> the previous lines
@@ -29,48 +29,48 @@ def rotate_an_image_array_line_by_line(
     Parameters
     ----------
     image_stack : np.ndarray
-        The image stack to be rotated.
+        The image stack to be derotated.
     rot_deg_line : np.ndarray
         The rotation angles by line.
 
     Returns
     -------
     np.ndarray
-        The rotated image stack.
+        The derotated image stack.
     """
 
     num_lines_per_frame = image_stack.shape[1]
-    rotated_image_stack = copy.deepcopy(image_stack)
+    derotated_image_stack = copy.deepcopy(image_stack)
     previous_image_completed = True
     rotation_completed = True
 
-    for i, rotation in tqdm.tqdm(
+    for i, angle in tqdm.tqdm(
         enumerate(rot_deg_line), total=len(rot_deg_line)
     ):
         line_counter = i % num_lines_per_frame
         image_counter = i // num_lines_per_frame
 
-        is_rotating = np.absolute(rotation) > 0.00001
+        is_rotating = np.absolute(angle) > 0.00001
         image_scanning_completed = line_counter == (num_lines_per_frame - 1)
         if i == 0:
             rotation_just_finished = False
         else:
             rotation_just_finished = not is_rotating and (
-                np.absolute(rot_deg_line[i - 1]) > np.absolute(rotation)
+                np.absolute(rot_deg_line[i - 1]) > np.absolute(angle)
             )
 
         if is_rotating:
             if rotation_completed and (line_counter != 0):
                 # when starting a new rotation in the middle of the image
-                rotated_filled_image = (
+                derotated_filled_image = (
                     np.ones_like(image_stack[image_counter])
                     * blank_pixels_value
                 )  # non sampled pixels are set to the min val of the image
-                rotated_filled_image[:line_counter] = image_stack[
+                derotated_filled_image[:line_counter] = image_stack[
                     image_counter
                 ][:line_counter]
             elif previous_image_completed:
-                rotated_filled_image = (
+                derotated_filled_image = (
                     np.ones_like(image_stack[image_counter])
                     * blank_pixels_value
                 )
@@ -83,16 +83,16 @@ def rotate_an_image_array_line_by_line(
             image_with_only_line = np.zeros_like(img_with_new_lines)
             image_with_only_line[line_counter] = line
 
-            rotated_line = rotate(
+            derotated_line = rotate(
                 image_with_only_line,
-                rotation,
+                angle,
                 reshape=False,
                 order=0,
                 mode="constant",
             )
 
-            rotated_filled_image = np.where(
-                rotated_line == 0, rotated_filled_image, rotated_line
+            derotated_filled_image = np.where(
+                derotated_line == 0, derotated_filled_image, derotated_line
             )
             previous_image_completed = False
         if (
@@ -101,11 +101,11 @@ def rotate_an_image_array_line_by_line(
             if rotation_just_finished:
                 rotation_completed = True
 
-                rotated_filled_image[line_counter + 1 :] = image_stack[
+                derotated_filled_image[line_counter + 1 :] = image_stack[
                     image_counter
                 ][line_counter + 1 :]
 
-            rotated_image_stack[image_counter] = rotated_filled_image
+            derotated_image_stack[image_counter] = derotated_filled_image
             previous_image_completed = True
 
-    return rotated_image_stack
+    return derotated_image_stack

derotation/simulate/basic_rotator.py

@@ -0,0 +1,151 @@
+import copy
+
+import numpy as np
+from scipy.ndimage import affine_transform
+
+
+class Rotator:
+    def __init__(self, angles: np.ndarray, image_stack: np.ndarray) -> None:
+        """Initializes the Rotator object.
+        The Rotator aims to imitate a 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.
+
+        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.
+
+        Parameters
+        ----------
+        angles : np.ndarray
+            An array of angles in degrees, representing the rotation of the
+            sample at the time of acquisition. The length of the array should
+            be equal to the number of lines per frame multiplied by the number
+            of frames.
+        image_stack : np.ndarray
+            The image stack represents the acquired images, as if there was no
+            rotation. The shape of the image stack should be (num_frames,
+            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.
+
+        Raises
+        ------
+        AssertionError
+            If the number of angles is not equal to the number of lines
+            per frame multiplied by the number of frames
+        """
+        #  there should be one angle per line pe frame
+        assert len(angles) == image_stack.shape[0] * image_stack.shape[1], (
+            f"Number of angles ({len(angles)}) should be equal to the number "
+            + "of lines per frame multiplied by the number of frames "
+            + 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]
+
+    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.
+
+        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.
+
+        Returns
+        -------
+        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)
+        blank_pixel = self.get_blank_pixels_value()
+
+        for i, image in enumerate(self.image_stack):
+            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)
+            )
+            if is_this_frame_rotating:
+                for j, angle in enumerate(self.angles[i]):
+                    if angle == 0:
+                        continue
+                    else:
+                        rotated_frame = self.rotate(image, angle, blank_pixel)
+                        rotated_image_stack[i][j] = rotated_frame[j]
+
+        return rotated_image_stack
+
+    @staticmethod
+    def rotate(
+        image: np.ndarray, angle: float, blank_pixel: float
+    ) -> np.ndarray:
+        """Rotate the entire image by a given angle. Uses affine transformation
+        with no interpolation.
+
+        Parameters
+        ----------
+        image : np.ndarray
+            The image to rotate.
+        angle : float
+            The angle in degrees to rotate the image in degrees. If positive,
+            rotates clockwise. If negative, rotates counterclockwise.
+        blank_pixel : float
+            The value to fill the blank pixels with.
+
+        Returns
+        -------
+        np.ndarray
+            The rotated image.
+        """
+        # Compute rotation in radians
+        angle_rad = np.deg2rad(angle)
+        cos, sin = np.cos(angle_rad), np.sin(angle_rad)
+
+        # Calculate center of the image
+        center_y, center_x = np.array(image.shape) / 2
+
+        # Rotation matrix clockwise if angle is positive
+        rotation_matrix = np.array(
+            [
+                [cos, -sin],
+                [sin, cos],
+            ]
+        )
+
+        # Compute offset so rotation is around the center
+        offset = np.array([center_y, center_x]) - rotation_matrix @ np.array(
+            [center_y, center_x]
+        )
+
+        # Apply affine transformation
+        rotated_image = affine_transform(
+            image,
+            rotation_matrix,
+            offset=offset,
+            output_shape=image.shape,  # Keep original shape
+            order=0,
+            mode="constant",  # NO interpolation
+            cval=blank_pixel,
+        )
+
+        return rotated_image
+
+    def get_blank_pixels_value(self) -> float:
+        """Get a default value to fill the edges of the rotated image.
+        This is necessary because we are using affine transformation with no
+        interpolation, so we need to fill the blank pixels with some value.
+
+        As for now, it returns the minimum value of the image stack.
+
+        Returns
+        -------
+        float
+            The minimum value of the image stack.
+        """
+        return np.min(self.image_stack)

examples/deformation_by_line_of_a_dog.py

@@ -3,7 +3,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from derotation.derotate_by_line import rotate_an_image_array_line_by_line
+from derotation.derotate_by_line import derotate_an_image_array_line_by_line
 
 path = Path(__file__).parent.parent / "images/dog.png"
 img = plt.imread(path)
@@ -14,7 +14,7 @@
 img_len = img.shape[0]
 rotation_angles = np.linspace(0, 180, img_len * 3)
 
-img_rotated = rotate_an_image_array_line_by_line(img_stack, rotation_angles)
+img_rotated = derotate_an_image_array_line_by_line(img_stack, rotation_angles)
 
 fig, ax = plt.subplots(1, 4, figsize=(10, 5))
 

examples/derotate.py

@@ -1,4 +1,4 @@
-from derotation.analysis.full_rotation_pipeline import FullPipeline
+from full_derotation_pipeline import FullPipeline
 
 derotate = FullPipeline("full_rotation")
 derotate()

examples/derotate_incremental.py

@@ -1,4 +1,4 @@
-from derotation.analysis.incremental_rotation_pipeline import (
+from incremental_derotation_pipeline import (
     IncrementalPipeline,
 )
 

examples/derotation_slurm_job.py

@@ -3,9 +3,8 @@
 from pathlib import Path
 
 import yaml
-
-from derotation.analysis.full_rotation_pipeline import FullPipeline
-from derotation.analysis.incremental_rotation_pipeline import (
+from full_derotation_pipeline import FullPipeline
+from incremental_derotation_pipeline import (
     IncrementalPipeline,
 )