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