From 8587817e7217cd88423aa6556116b80e49ae4cfa Mon Sep 17 00:00:00 2001 From: lauraporta <29216006+lauraporta@users.noreply.github.com> Date: Thu, 17 Oct 2024 21:21:07 +0100 Subject: [PATCH] Apply suggestions from PR review --- derotation/simulate/basic_rotator.py | 83 +++++++++++++++++------- examples/rotate_and_derotate_a_square.py | 20 +++++- 2 files changed, 80 insertions(+), 23 deletions(-) diff --git a/derotation/simulate/basic_rotator.py b/derotation/simulate/basic_rotator.py index 6f5f207..a9d1583 100644 --- a/derotation/simulate/basic_rotator.py +++ b/derotation/simulate/basic_rotator.py @@ -7,13 +7,28 @@ 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 - The rotation angles by line per frame + 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 to be rotated + 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 ------ @@ -28,46 +43,66 @@ def __init__(self, angles: np.ndarray, image_stack: np.ndarray) -> None: + f"({image_stack.shape[0] * image_stack.shape[1]})" ) - self.angles = angles + # 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: - """Rotates the image stack line by line, using the rotation angles - provided. + """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. + 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): - start_angle_idx = i * self.num_lines_per_frame - end_angle_idx = self.num_lines_per_frame * (i + 1) - - is_this_frame_rotating = np.any( - np.abs(self.angles[start_angle_idx:end_angle_idx]) > 0.00001 + 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: - frame = copy.deepcopy(image) - for j, angle in enumerate( - self.angles[start_angle_idx:end_angle_idx] - ): + for j, angle in enumerate(self.angles[i]): if angle == 0: continue else: rotated_frame = self.rotate(image, angle, blank_pixel) - frame[j] = rotated_frame[j] - rotated_image_stack[i] = frame + rotated_image_stack[i][j] = rotated_frame[j] return rotated_image_stack + @staticmethod def rotate( - self, image: np.ndarray, angle: float, blank_pixel: float + 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) @@ -75,7 +110,7 @@ def rotate( # Calculate center of the image center_y, center_x = np.array(image.shape) / 2 - # Rotation matrix + # Rotation matrix clockwise if angle is positive rotation_matrix = np.array( [ [cos, -sin], @@ -94,15 +129,19 @@ def rotate( rotation_matrix, offset=offset, output_shape=image.shape, # Keep original shape - order=0, # NO interpolation - mode="constant", + order=0, + mode="constant", # NO interpolation cval=blank_pixel, ) return rotated_image def get_blank_pixels_value(self) -> float: - """Returns the minimum value of the image stack. + """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 ------- diff --git a/examples/rotate_and_derotate_a_square.py b/examples/rotate_and_derotate_a_square.py index ce84b48..dc7ec52 100644 --- a/examples/rotate_and_derotate_a_square.py +++ b/examples/rotate_and_derotate_a_square.py @@ -5,10 +5,12 @@ from derotation.simulate.basic_rotator import Rotator # make a simple image, a square in a black background -image = np.zeros((100, 100)) +image = np.empty((100, 100)) +# it will have 5 gray levels to visualize better the rotation gray_values = [i % 5 * 100 + 155 for i in range(100)] for i in range(100): image[i] = gray_values[i] +# make a black border image[:20] = 0 image[-20:] = 0 image[:, :20] = 0 @@ -22,6 +24,7 @@ num_angles = image_stack.shape[0] * image_stack.shape[1] angles = np.arange(num_angles) +# rotate the image stack rotator = Rotator(angles, image_stack) rotated_image_stack = rotator.rotate_by_line() @@ -30,6 +33,7 @@ rotated_image_stack, angles ) +# plot the original image, the rotated images and the derotated images fig, ax = plt.subplots(2, 4, figsize=(20, 5)) ax[0, 0].imshow(image, cmap="gray") @@ -41,6 +45,19 @@ ax[0, i + 1].set_title(f"Rotated image {i + 1}") ax[0, i + 1].axis("off") + # add information about the rotation angle + angles = rotator.angles[i] + angle_range = f"angles: {angles.min():.0f}-{angles.max():.0f}" + ax[0, i + 1].text( + 0.5, + 0.9, + angle_range, + horizontalalignment="center", + verticalalignment="center", + transform=ax[0, i + 1].transAxes, + color="white", + ) + ax[1, 0].imshow(image, cmap="gray") ax[1, 0].set_title("Original image") ax[1, 0].axis("off") @@ -50,4 +67,5 @@ ax[1, i + 1].set_title(f"Derotated image {i + 1}") ax[1, i + 1].axis("off") + plt.show()