Skip to content

Commit

Permalink
Apply suggestions from Igor
Browse files Browse the repository at this point in the history
  • Loading branch information
@lauraporta
lauraporta committed Dec 17, 2024
1 parent 4e855e8 commit 4129dac
Show file tree
Hide file tree
Showing 3 changed files with 77 additions and 97 deletions.
142 changes: 59 additions & 83 deletions derotation/simulate/line_scanning_microscope.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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.
Expand All @@ -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)
Expand Down Expand Up @@ -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}")
Expand Down Expand Up @@ -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],
]
)
Expand All @@ -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
Expand Down Expand Up @@ -202,83 +193,68 @@ 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
-------
np.ndarray
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
Expand Down
4 changes: 4 additions & 0 deletions examples/elliptical_rotations.py
View file
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -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
)
Expand Down
28 changes: 14 additions & 14 deletions tests/test_integration/test_rotation_out_of_plane.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
from pathlib import Path
from typing import Optional, Tuple
from typing import Tuple

import numpy as np
import pytest
Expand All @@ -17,24 +17,24 @@


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
optional orientation.
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
Expand Down Expand Up @@ -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:
"""
Expand All @@ -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
Expand Down

0 comments on commit 4129dac

Please sign in to comment.