diff --git a/derotation/analysis/fit_ellipse.py b/derotation/analysis/fit_ellipse.py index 46ec9a9..f92a887 100644 --- a/derotation/analysis/fit_ellipse.py +++ b/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) diff --git a/derotation/simulate/line_scanning_microscope.py b/derotation/simulate/line_scanning_microscope.py index 00476d3..e94ce3d 100644 --- a/derotation/simulate/line_scanning_microscope.py +++ b/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( diff --git a/examples/elliptical_rotations.py b/examples/elliptical_rotations.py new file mode 100644 index 0000000..975de05 --- /dev/null +++ b/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", +) diff --git a/tests/test_integration/test_finding_center_of_rotation.py b/tests/test_integration/test_finding_center_of_rotation.py index 46c38d1..4d229ef 100644 --- a/tests/test_integration/test_finding_center_of_rotation.py +++ b/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 ------- diff --git a/tests/test_integration/test_rotation_out_of_plane.py b/tests/test_integration/test_rotation_out_of_plane.py index 78c7d35..ec40248 100644 --- a/tests/test_integration/test_rotation_out_of_plane.py +++ b/tests/test_integration/test_rotation_out_of_plane.py @@ -3,8 +3,6 @@ import numpy as np import pytest -from matplotlib import pyplot as plt -from skimage.draw import ellipse_perimeter from derotation.analysis.fit_ellipse import ( fit_ellipse_to_points, @@ -65,73 +63,6 @@ def rotate_image_stack( return image_stack, rotated_image_stack, rotator, num_frames -def make_plot( - image_stack: np.ndarray, - rotated_image_stack: np.ndarray, - rotator: Rotator, - num_frames: int, - title: str = "", -) -> None: - """ - Plots the original and rotated image stacks. - - 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") - - @pytest.mark.parametrize( "plane_angle,expected_orientation", [ @@ -245,58 +176,3 @@ def test_max_projection( f"Difference between major and minor axes should be close to " f"{expected_b}, instead got {b}" ) - - -def plot_max_projection_with_rotation_out_of_plane( - max_projection, yc, xc, a, b, o, title -): - # plot for debugging purposes - - fig, ax = plt.subplots(1, 2, figsize=(10, 5)) - ax[0].imshow(max_projection, cmap="gray") - ax[0].set_title("Max projection") - ax[0].axis("off") - - # # plot edges - # ax[1].imshow(edges, cmap="gray") - # ax[1].set_title("Edges") - # ax[1].axis("off") - - print(f"Ellipse center: ({yc}, {xc}), a: {a}, b: {b}, orientation: {o}") - - # plot the ellipse - cy, cx = ellipse_perimeter(yc, xc, a, b, orientation=o) - empty_image = np.zeros_like(max_projection) - empty_image[cy, cx] = 255 - - # show ellipse - ax[1].imshow(empty_image, cmap="gray") - ax[1].set_title("Fitted ellipse") - ax[1].axis("off") - - plt.savefig(title) - - -if __name__ == "__main__": - # To be used for debugging purposes - 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", - )